Inverse polylogarithm

[Ondřej Čertík]

Hi,

I would like to evaluate to higher accuracy (quadruple precision is
enough) an inverse Fermi-Dirac integral of order 1/2. The direct
function is the following integral (and it can also be written as a
polylogarithm):

I_{1/2}(x) = Integral(sqrt(t) / (1+exp(t-x)), (t, 0, oo))
= -gamma(S(3)/2) * polylog(S(3)/2, -exp(x))

and I need its inverse. SymPy can evaluate either representation:

In [19]: (-gamma(S(3)/2) * polylog(S(3)/2, -exp(x))).subs(x, 3).n(35, chop=True)
Out[19]: 3.9769853540479774178558497377805661

In [20]: Integral(sqrt(t) / (1+exp(t-x)), (t, 0, oo)).subs(x, 3).n(35)
Out[20]: 3.9769853540479774178558497377805661

I need to calculate an inverse, i.e. for 3.9769... it would return
3.00..., to arbitrary accuracy. Once I have that, I want to use
rational function approximation to obtain very fast and accurate
double precision implementation.

What is the best way to do that?

There is a recent publication [1], which does that using Mathematica,
and then provides Fortran implementation. I tested it and it works
great. What I actually need is to find a rational approximation to an
expression that contains both the inverse and direct Fermi-Dirac
integrals and I want to just have a rational approximation for the
final expression. The inverse Fermi-Dirac is the hardest part, so
that's why I am asking.

In [1], they calculate a series expansion of the Fermi-Dirac integral
and then they reverse (invert) the series, see eq. (9), (10), (14),
(15) and the Mathematica code afterwards. Unfortunately, SymPy raises
an exception for a series of the polylog function:

In [23]: polylog(S(3)/2, -z).series(z, 0, 5)

https://github.com/sympy/sympy/issues/9497

But let's say we fix it. Here i s a slightly modified expression that works:

In [6]: polylog(S(3)/2, 1+z).series(z, 0, 5)
zeta(3/2) + z*zeta(1/2) + z**2*(zeta(-1/2)/2 - zeta(1/2)/2) +
z**3*(zeta(1/2)/3 + zeta(-3/2)/6 - zeta(-1/2)/2) +
z**4*(11*zeta(-1/2)/24 + zeta(-5/2)/24 - zeta(-3/2)/4 - zeta(1/2)/4) +
O(z**5)

How can we invert the series? Mathematica has a function InverseSeries:

http://reference.wolfram.com/language/ref/InverseSeries.html

It seems SymPy can't do it yet.

Shivam, I think that would be a very useful addition to the series
module that you are implementing. We have to figure out an algorithm
for it.

Ondrej


[1] Fukushima, T. (2015). Precise and fast computation of inverse
Fermi–Dirac integral of order 1/2 by minimax rational function
approximation. Applied Mathematics and Computation, 259, 698–707.
doi:10.1016/j.amc.2015.03.015

[Aaron Meurer]

I think there's a closed form for the inversion of a formal power
series. See https://en.wikipedia.org/wiki/Formal_power_series#The_Lagrange_inversion_formula.

I'm not sure what the algorithm is for series like x*sin(x), which
that Mathematica example shows is a series in sqrt(x). I suppose
that's because the linear term is 0 (the Wikipedia article claims that
the constant needs to be 0 and the linear term needs to be nonzero).
Strictly speaking x*sin(x) is not invertible near x=0 (it's an even
function). I guess this picks one branch and inverts that.

Aaron Meurer

>
> Ondrej
>
>
> [1] Fukushima, T. (2015). Precise and fast computation of inverse
> Fermi–Dirac integral of order 1/2 by minimax rational function
> approximation. Applied Mathematics and Computation, 259, 698–707.
> doi:10.1016/j.amc.2015.03.015
>

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CADDwiVCp7rfTUk0TrooNaWqpVrADFbGPb_Ptup2Jn49twTapEg%40mail.gmail.com.
> For more options, visit https://groups.google.com/d/optout.

[Ondřej Čertík]

I don't quite get it with the constant term, but it seems to be true,
e.g. the inversion of exp(x) = 1 + x + .... doesn't exist, because
that would be log(x), which doesn't have an expansion around x=0. One
can only expand e.g. log(1+x), as the inverse is exp(x)-1.

I found a pretty good intro here:

http://mathworld.wolfram.com/SeriesReversion.html

with further references to literature.

>
> Aaron Meurer
>
>>
>> Ondrej
>>
>>
>> [1] Fukushima, T. (2015). Precise and fast computation of inverse
>> Fermi–Dirac integral of order 1/2 by minimax rational function
>> approximation. Applied Mathematics and Computation, 259, 698–707.
>> doi:10.1016/j.amc.2015.03.015
>>
>> --
>> You received this message because you are subscribed to the Google Groups "sympy" group.
>> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
>> To post to this group, send email to sy...@googlegroups.com.
>> Visit this group at http://groups.google.com/group/sympy.
>> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CADDwiVCp7rfTUk0TrooNaWqpVrADFbGPb_Ptup2Jn49twTapEg%40mail.gmail.com.
>> For more options, visit https://groups.google.com/d/optout.
>
> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAKgW%3D6K%2BXXF6-Bb7tjxXHMbt74wEG_ufEu_%3D48-2UnQ4YbJjFw%40mail.gmail.com.

[mario]

There is  `series_reversion` in PR609.

[Ondřej Čertík]

On Wed, Jun 10, 2015 at 9:30 AM, mario <mario....@gmail.com> wrote:
>
> There is `series_reversion` in PR609.

Indeed, it's there:

https://github.com/sympy/sympy/pull/609/files#diff-c36c26d9ea5b82b88777bcbc0281eefdR823

Thanks Mario. Shivam, let's port this into the ring_series as well.

Ondrej

[Shivam Vats]

Mario, the algorithm looks quite similar to the Newton

method you have used for functional inverses for

expanding tan and tanh. We could also expand `sin/cos`.

Which one do think should be faster?

Also, can the convergence of the method be proven?

I could not find any reference to it (except the algorithm
in Computer Arithmetic by Paul Zimmerman).

Shivam Vats

[Shivam Vats]

Sure, Ondrej. I will send a PR for it.

[mario]

As long as the series to be inverted has finite precision `n`, has no constant term and finite linear
term, one can revert the series iteratively; I do not think there are convergence problems,
since these series are formal,
or anyway for sufficienly small region around 0, the function is monotonic, so it can be reversed.

[Fredrik Johansson]

On Wednesday, June 10, 2015 at 5:53:16 PM UTC+2, Shivam Vats wrote:

Mario, the algorithm looks quite similar to the Newton

method you have used for functional inverses for

expanding tan and tanh. We could also expand `sin/cos`.

Which one do think should be faster?

Also, can the convergence of the method be proven?

I could not find any reference to it (except the algorithm
in Computer Arithmetic by Paul Zimmerman).

Series reversion is a purely formal operation, so it always converges.

For a faster algorithm to compute the reversion of a generic power series, see my paper http://fredrikj.net/math/reversion.html (a "generic" power series being given by a list of coefficients -- there are better methods for computing inverse series of functions composed from elementary functions like sin(sin(x))).

Fredrik

[Shivam Vats]

Thanks a lot Fredrik! This looks very interesting! I will
definitely read your paper. Could you point me to references
for inverting compositions of elementary functions?

Regards
Shivam Vats

[Ondřej Čertík]

On Tue, Jun 16, 2015 at 9:02 AM, Shivam Vats <shiv...@gmail.com> wrote:
>
> Thanks a lot Fredrik! This looks very interesting! I will
> definitely read your paper. Could you point me to references
> for inverting compositions of elementary functions?

Is your question related to series somehow? What kind of compositions
do you have in mind?

Ondrej

[Fredrik Johansson]

On Tuesday, June 16, 2015 at 5:02:28 PM UTC+2, Shivam Vats wrote:

Thanks a lot Fredrik! This looks very interesting! I will
definitely read your paper. Could you point me to references
for inverting compositions of elementary functions?

The general idea is that computing logarithmic and inverse trigonometric functions of formal power series is just algebraic operations on power series followed by formal (term by term) integration, e.g. log(f(x)) = int f'(x) / f(x) dx. From there, Newton iteration allows you to compute exponential and forward trigonometric functions.

The Brent-Kung paper is probably the best place to start. As further reading, see the other references in my reversion paper (there's also a little more content in section 4 of my PhD thesis).

Fredrik

[Shivam Vats]

Is your question related to series somehow? What kind of compositions
do you have in mind?

Yes, Ondrej. Since, we have already implemented many elementary
functions in ring_series, it will be nice if we have faster methods
to invert them and their compositions (eg sin(sin(x)))

The general idea is that computing logarithmic and inverse trigonometric functions of formal power series is just algebraic operations on power series followed by formal (term by term) integration, e.g. log(f(x)) = int f'(x) / f(x) dx. From there, Newton iteration allows you to compute exponential and forward trigonometric functions.

Right! We have attempted something similar in `ring_series`.
Formula based expansion as `_atan_series` and this method
as `rs_atan` are implemented here. `_atan_series` seems to
be much faster than `rs_atan`. As of now, we expand `tan`
from `atan` using Newton iterations and sin/cos from `tan`
using half angle formula. What do you suggest?

The Brent-Kung paper is probably the best place to start. As further reading, see the other references in my reversion paper (there's also a little more content in section 4 of my PhD thesis).

Thanks a lot Fredrik! These references are extremely
helpful :)

Regards
Shivam Vats

 

[Fredrik Johansson]

On Wednesday, June 17, 2015 at 9:04:54 PM UTC+2, Shivam Vats wrote:

Is your question related to series somehow? What kind of compositions
do you have in mind?

Yes, Ondrej. Since, we have already implemented many elementary
functions in ring_series, it will be nice if we have faster methods
to invert them and their compositions (eg sin(sin(x)))

The general idea is that computing logarithmic and inverse trigonometric functions of formal power series is just algebraic operations on power series followed by formal (term by term) integration, e.g. log(f(x)) = int f'(x) / f(x) dx. From there, Newton iteration allows you to compute exponential and forward trigonometric functions.

Right! We have attempted something similar in `ring_series`.
Formula based expansion as `_atan_series` and this method
as `rs_atan` are implemented here. `_atan_series` seems to
be much faster than `rs_atan`. As of now, we expand `tan`
from `atan` using Newton iterations and sin/cos from `tan`
using half angle formula. What do you suggest?

Yes, for atan this is the right approach.

For sin/cos, exp, tan, etc. Newton iteration is optimal for long power series if you have asymptotically fast polynomial multiplication. When your multiplication is O(n^2), the following O(mn) algorithm for exponentials (where m is the length of the input and n is the length of the output) is better:

    def exp_series(A, n):
        B = [exp(A[0])]
        for k in range(1, n):
            B.append(sum(j*A[j]*B[k-j]/k for j in range(1,min(len(A),k+1))))
        return B

It's easy to turn this into an algorithm for sin/cos (and hence also tan), and there is an analog for raising a series to a constant power (if I remember correctly, it's in Knuth vol 2).

Fredrik