Numerical instability when trying to use lambdify

[Olivier Grisel]

Hi all,

I use sympy.solve to find the roots of a degree 3 polynomial and then would like to convert the solution as a numpy expression. However lambdify does cannot compute a numerically stable evaluation (whereas evalf has no issue).

Here is a gist giving the details and a small reproduction script along with the observed output.

  https://gist.github.com/1582089

Any hint?

--
Olivier

[krastano...@gmail.com]

Hi,

I hope that I am not misleading you. Hopefully the developers will correct me if I'm wrong.

evalf uses the mpmath library for doing multiprecision math. It can not overflow under normal conditions. On the other hand lambdify will exchange by default all sympy functions with numpy equivalents. Those can overflow. Try adding "mpmath" as a third argument to lambdify. That way it will use mpmath instead of numpy. It will be slower than numpy, but it will not overflow. Check the docstring for details.

Stefan

--
Olivier

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To view this discussion on the web visit https://groups.google.com/d/msg/sympy/-/Kg9vhQMHIZQJ.
To post to this group, send email to sy...@googlegroups.com.
To unsubscribe from this group, send email to sympy+un...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/sympy?hl=en.

[Olivier Grisel]

Thank you very much it indeed works with "mpmath" in lambdify. However my final goal would be to actually use numpy to evaluate that value (as a function of n and d) in another program that has no dependency on sympy nor mpmath.

Do you think there is a way to make sympy refactor the expression so that the naive numpy translation would be more stable, especially for the range of d and n is such that the solution has an imaginary part close to or exactly zero?

I will try to further investigate.

--
Olivier

[krastano...@gmail.com]

On 10 January 2012 12:56, Olivier Grisel <olivier...@gmail.com> wrote:

Thank you very much it indeed works with "mpmath" in lambdify. However my final goal would be to actually use numpy to evaluate that value (as a function of n and d) in another program that has no dependency on sympy nor mpmath.

Do you think there is a way to make sympy refactor the expression so that the naive numpy translation would be more stable, especially for the range of d and n is such that the solution has an imaginary part close to or exactly zero?

I don't think it can be done automatically by sympy. But there are many functions that rewrite expressions. I know too little about them to tell but maybe some of them will help. Hopefully someone on the group will know a better way to do it. 

I will try to further investigate.

--
Olivier

--

You received this message because you are subscribed to the Google Groups "sympy" group.

To view this discussion on the web visit https://groups.google.com/d/msg/sympy/-/yQddw56BoSwJ.

[Olivier Grisel]

On Tuesday, January 10, 2012 2:14:50 PM UTC+1, Stefan Krastanov wrote:

On 10 January 2012 12:56, Olivier Grisel <olivier...@gmail.com> wrote:

Thank you very much it indeed works with "mpmath" in lambdify. However my final goal would be to actually use numpy to evaluate that value (as a function of n and d) in another program that has no dependency on sympy nor mpmath.

Do you think there is a way to make sympy refactor the expression so that the naive numpy translation would be more stable, especially for the range of d and n is such that the solution has an imaginary part close to or exactly zero?

I don't think it can be done automatically by sympy. But there are many functions that rewrite expressions. I know too little about them to tell but maybe some of them will help. Hopefully someone on the group will know a better way to do it. 

I might be able to come up with something that has no dependency on any external library by extracting the real and imaginary parts of the root separately using sympy.re / sympy.im and then translate them to a pair of python expressions that use http://docs.python.org/library/decimal.html for arbitrary precision floats.

Thanks for your help, I will let you know when I'll have time to carry on with that experiment.

--
Olivier

[Olivier Grisel]

Just a followup:

Apparently sympy (0.7.1) is unable to find an expansion of the real part of the root that interests me: it yields a RuntimeError:

>>> sp.expand(sol, complex=True)

=> RuntimeError: maximum recursion depth exceeded while calling a Python object

Apparently wolframlapha is able to find a something interesting in terms of phase:

http://www.wolframalpha.com/input/?i=Re[-%28%281%2Bi+sqrt%283%29%29+d%29%2F%284+%28d^3-48+d^2+log%28n%29%2B4+sqrt%286%29+sqrt%2824+d^4+log^2%28n%29-d^5+log%28n%29%29%29^%281%2F3%29%29-%28%281-i+sqrt%283%29%29+%28d^3-48+d^2+log%28n%29%2B4+sqrt%286%29+sqrt%2824+d^4+log^2%28n%29-d^5+log%28n%29%29%29^%281%2F3%29%29%2F%284+d%29%2B1%2F2]

I wonder whether the fact that sympy is crashing with a recursion depth exceeded is a bug or this a consequence of a missing algorithm only found in mathematica.

[Chris Smith]

Could you post the sympy imput because I get:

>>> -(1+I*sqrt(3))*d/(4*d)

-1/4 - sqrt(3)*I/4

>>> re(_)

-1/4

 

[Olivier Grisel]

Hi smichr,

I have updated my original gist:

https://gist.github.com/1582089

Here is a summary:

"""
import sympy as sp
import numpy as np

d, n = sp.symbols('d n', real=True, positive=True)
eps = sp.symbols('eps', real=True, positive=True)

jl_bound = 4 * sp.log(n) / (eps ** 2 / 2 - eps ** 3 / 3)

# solve the inverted expression as sympy is much faster as finding the roots of
# a polynomial and the point eps=0 is not interesting to us anyway.
solutions = sp.solve(sp.Equality(1 / d, 1 / jl_bound), eps)

sol = solutions[1]

sp.expand(sol, complex=True)

"""

[Chris Smith]

You solutions contain symbols n and d and they don't appear to cancel as nicely as in the expression that you showed as input to Mathematica. 

[Olivier Grisel]

I thought about that too and tried to replace them in the sympy program by a new symbol a = log(n) / d and tell sympy to treat it as a positive real but it does not help. I still get the RuntimeError when trying to tell the real part from the imaginary part of the root.