JuliaMath
Meijer G function
Hi there, it would be possible to add the meijer-G function to the Special Functions package? Thanks a lot! Crist.
I doubt it — it seems like one of those incredibly general functions (look at how many input parameters it has!) that is virtually impossible to implement well.
Is there any efficient implementation of this function (i.e. not just summing a slowly converging series or doing a brute-force numerical integration) in any language or any publication anywhere?
Hi! Thank you for your response. That is true, there is a lot of input parameters. To be honest I don't know how internally it operates, but this function is implemented in Mathematica and matlab and it seems to work well. Best, Crist.
Symbolic-algebra packages like Mathematica use slow brute-force methods to compute special functions. It looks like Matlab may use its symbolic-algebra backend for this, so it is probably similar.
(A simple test: in Matlab, try computing this function on rand(10^6,1) and timing it, and compare it to evaluating something like a Bessel function on the same inputs.)
Here is a test on Matlab R2020b:
>> x = rand(10^5,1);
>> tic; meijerG(3, [], [], 2, x); toc
Elapsed time is 39.323138 seconds.
>> tic; expint(x); toc
Elapsed time is 0.092412 seconds.
>> tic; besselj(2,x); toc
Elapsed time is 0.065117 seconds.
i.e. it is 500–1000× slower than other special functions, which is indicative of a slow brute-force calculation (e.g. by numerical integration). And if you add more parameters, it slows down even more — e.g. meijerG(3, [4.6, 8.7], [3.2, 4.1], 2, x) is even slower, by an another factor of 25×.
Nearly the whole point of using a special function is to avoid slow brute-force numerical integration, so if that's the only way to compute it then I don't think it's appropriate for this package.
Dear Steven, Indeed this functions is so slow, even in matlab, and that was my reason to try it in Julia. But now, with your explanation, I see that run time depends of the internal implementation, not the language. I have been using this function to calculate the fractional discrete laplacian matrix for a two dimensional field. Apparently, it would be possible to obtain the same result by integrating a function written in terms of special functions. But when I do that the answer is not the same provided by meijerG, which is the correct one. In spite of the distribution of the result is correct something is wrong with the "intensity", probably it comes from the lower limit of integral which is "0". It is worth to mention that by using the meijerG function the elapsed time is around 1 hour, however, by using the another definition in Julia it takes around 20 secs. Best, Cristian.
The Matlab special functions are mostly not implemented in Matlab, but rather are calling C and Fortran libraries, so their performance is usually reasonably good on a vector of inputs.
@stevengj I'm very interested in using these and I can work with you to add them. I have a use-case for them in SymbolicRegression.jl. Because Meijer G-functions let you represent a very general class of special functions, they are a useful basis for evolving new operators.
I think that in Julia we could exploit multiple dispatch, specializing on the parameters (which are only integers - i.e., via Val{p}), to perform some initial symbolic simplification before computing the integral.
Also, here's a note I found on the Mathematica implementation from stackexchange:
According to this old draft by Folkmar Bornemann, the implementation of Meijer G-function is a bit complicated and it happens in various stages. The article mentions a comment by Daniel Lichtblau on July 27, 2003:
I'll be discussing aspects of MeijerG issues at ACA next week. It is basically a lookup that converts various functions to MeijerG, then figures out the integral of a product of 2 MeijerG's via Slater convolution.
Furthermore, the Notes on Internal Implementation say the following:
Many other definite integrals are done using Marichev-- Adamchik Mellin transform methods. The results are often initially expressed in terms of Meijer G functions, which are converted into hypergeometric functions using Slater's Theorem and then simplified.
So in a nutshell, when one asks for the value of a Meijer G-function, this is converted to hypergeometric functions and then evaluated using the latter. The attached article also contains some mathematical details, that an interested reader may enjoy.
So perhaps there are faster ways to compute them?
What do you think? Best, Miles
