Laguerre expansion of multivariate function?

[mvsoom]

I came across your paper when I was looking for a Laguerre expansion of a multivariate function, which led me to this Julia package.

My question is, given a function k(x,y): ℝ₊² → ℝ, i.e. a function that may be positive and/or negative on the first quadrant, could I use this package to expand k(x,y) ≈ ∑ₘ ∑ₙ aₘₙ Lₘ(x) Lₙ(y) exp(-x/2) exp(-y/2), i.e., calculate the aₘₙ coefficients? Here Lₘ(x) is the Laguerre polynomial of order m.

I know such 2D expansions are possible with Chebyshev and Legendre polynomials, and from the paper it looks like it might possible with Laguerre polynomials too, but I wanted to doublecheck before diving into the math. Note the k(x,y) is not a PDF; it could be negative at some values of (x,y) and need not integrate to 1.

Thank you!

[lrnv]

Hey,

Thanks for reading the paper ! Yes you can, of course, this is the basis of the approximations in the paper. This package is very unstable yet, but there is code in it to do exactly what you ask for with two lines of Julia.

The standard in this package is to use k(x,y) ≈ ∑ₘ ∑ₙ aₘₙ Lₘ(2x) Lₙ(2y) exp(-x) exp(-y), as it is in my paper, but the translation to what you want is going to be straightforward.

Some functions should be reworked a little since I assume everywhere than $f$ is a density, but indeed the basis is more general. Do you want me to give you a code example ?

[mvsoom]

Thanks for the reply. That's some good news :)

Do you want me to give you a code example ?

That would be great!

[lrnv]

To compute theoretical coefficients, I use DoubleExponentielFormulas.jl, an unregistered package because it implements takashi-mori cubatures in pure Julia. It goes as follows:

#using Pkg
#Pkg.add(url="https://github.com/machakann/DoubleExponentialFormulas.jl.git")
#Pkg.add(url="https://github.com/lrnv/ThorinDistributions.jl.git")
import ThorinDistributions as TD
import DoubleExponentialFormulas as DEF

F = Float64 # You could use Double64 from DoubleFloats or even BigFloat f you do not trust the results. 

# Building the integrator and loading laguerre things : 
qde = DEF.QuadDE(F; h0=one(F)/8, maxlevel=10)
ϕ(x,p) = TD.laguerre_phi(x,p)
coef(fun, p) = qde(y -> qde(x -> fun(x,y)*ϕ(x,p[1])*ϕ(y,p[2]),0,Inf)[1], 0, Inf)

# Applicaiton: 
f(x,y) = exp(-x-y) # the function 

coef(f,[0,0])
coef(f,[1,0])
coef(f,[0,1])
coef(f,[1,1])

Although this is somewhat slow, it allows to change the type of the integration to be sure to have correct results if your function is 'hard' to integrate. it becomes much faster if you only want univariate coefficients.

If the function $f$ is a density, sampling it is usually a much better solution.

[mvsoom]

OK, that's great, thank you for the answer!