Holonomic_functions.md

+John Cook I like hypergeometric functions, but when you have an expression (say in SymPy) and want to convert it to a hypergeometric function (or a linear combination of those), it turns out it's really difficult, because when you multiply two hypergeometric functions, the result is generally not hypergeometric, except few simple cases. So you are totally screwed. However, just a few days ago, +Fredrik Johansson introduced me to holonomic functions (https://en.wikipedia.org/wiki/Holonomic_function), which are generalizations of hypergeometric functions (and still closed under differentiation and integration), but you can multiply them. I am more and more convinced that this is the ultimate way to go, and I would like to implement it into SymPy.

You mention that the definition seems arbitrary. That is not so, but it took me a long time to connect the dots. The overall idea is simple and it really makes sense. I have numbered the individual steps, as you generalize things in a very natural, I would say obvious once you see this (not arbitrary) way:

1. You start with a series sum t_k, and you say that this is geometric if t_{k+1}/t_k is a constant (it can depend on "x", but not on "k").

2. So the most simple generalization is to go from a constant to a rational function (ratio of two arbitrary polynomials). So you define a hypergeometric series if t_{k+1}/t_k is a rational function.

3. The problem is that the series doesn't always converge (you need to analytically continue it, and sometimes the series is just asymptotic), so an equivalent and more robust way is to realize that the hypergeometric function pFq(a, b, x) satisfies a differential equation:

[(xd/dx+a1)(xd/dx+a2)... - d/dx(xd/dx+b1-1)(xd/dx+b2-1)...] f(x) = 0

   With the initial conditions at x=0, which are given by:

pFq(a, b, 0) = 1 pFq'(a, b, 0) = a1a2... / (b1b2*...) pFq''(a, b, 0) = a1*(a1+1) * a2*(a2+1) ... / (b1*(b1+1) * b2*(b2+1) ...) ...

   So things are uniquely defined, and you just need to integrate the differential equation, using any usual numerical method. Things are then well defined for all complex "x", no issues with the possibly finite or zero radius of convergence of the hypergeometric series at x=0. This step 3. is equivalent to 2., but you realize that the differential equation formulation is way better, no more special cases, you can use the usual methods for ODEs.

   When you multiply the factors out, let's say for the Gauss hypergeometric function 2F1, and write the differential equation in the usual form:

x*(1-x)f' '(x) + (c-(a+b+1)x)f'(x) - abf(x) = 0

   You notice that all these equations have quite a special form. It's almost constant coefficients, with a very restricted powers of "x" attached to each derivative.

4. As we have done in step 2., the most simple generalization is to go from a constant (or almost a constant in this case) to a rational function, i.e. a ratio of two polynomials. Since you can multiply the differential equation by the denominators, this is equivalent to generalizing to a linear homogeneous differential equations with polynomial coefficients, e.g.:

a2(x)*f' '(x) + a1(x)*f'(x) + a0(x)*f(x) = 0

   And these are called holonomic functions. The function is uniquely defined if you specify the polynomials a0(x), a1(x), a2(x), ... and a set of initial conditions (values and derivatives), say at x=0. Since the differential equation for a holonomic function can be written as:

[ Sum_{i,j } A(i, j) * x^j * d^i/dx^i ] * f(x) = 0

   in order to represent any holonomic function on a computer, you just need to store the matrix A(i, j) of the coefficients in the polynomials, typically integers, and a vector of initial conditions. For example, to represent cos(x), you get A = [[1], [0], [1]] and the initial conditions at x=0 are X0=[1, 0]. So this is really not much different from the hypergeometric function, where you also need to store the vectors of coefficients "a" and "b", for example cos(x) =  0F1(3/2; x^2/4), so you get a=[], b=[3/2] and z=x^2/4. So when you store on a computer, for hypergeometric function you need to be able to represent and store z=x^2/4, but for a holonomic function you just store the matrix A + vector X0 of initial conditions. So it's definitely not more complicated, I would say it's actually simpler. Then you have algorithms for adding, multiplying, differentiating, integrating (and other operations) that just transform A and X0 to another A and X0, these algorithms can be implemented in a computer (some are straightforward, some not so much, but all can be done).

----------------

Since hypergeometric functions are obviously a subset, the holonomic functions include all the usual elementary and special functions, you can trivially integrate and differentiate, but also you can multiply them (and add them). Not all functions are holonomic, the wikipedia has some examples, but holonomic functions cover a lot, and yet they can still be manipulated algorithmically on a computer. You can convert any hypergeometric function to a holonomic one trivially, e.g. here it is:

     def hyper(x, Dx, a, b, z):
        r1 = 1
        for i in range(len(a)):
            r1 = r1*(x*Dx+a[i])
        r2 = Dx
        for i in range(len(b)):
            r2 = r2*(x*Dx+b[i]-1)
        return (r1-r2).annihilator_of_composition(z)

this is using the ore_algebra module for Sage (http://www.risc.jku.at/research/combinat/software/ore_algebra/) written by +Fredrik Johansson- and others.

To convert a holonomic function back to a hypergeometric one is also quite easy, as one can expand the holonomic function into a (formal) series, and then compute the ratio of terms and get the hypergeometric function out of it.

I am hoping there is a way to also get a linear combination of hypergeometric functions, that would cover a lot of cases.

For the rest we'll probably have to have some kind of a table or pattern matching. But ultimately it's a similar problem as converting a hypergeometric function back to elementary ones. Even if you don't know how to do that, the fact that you can evaluate the function numerically from the differential equation makes the function essentially "known". Similar to if you have an algebraic number as a root of a polynomial, even if you can't or don't know how to convert to radicals, you can still evaluate numerically, or calculate symbolically as a root of a polynomial (RootOf in SymPy or Mathematica).

Anyway, I just learned about this two days ago, and I would really like to have a good support for this in SymPy.