Incorrect results with Sympy when utilizing (numpy's) floats

Question

I am trying to calculate a velocity tensor from a time dependent rotationmatrix RE(t) (Namely the earth rotation at latitude 48.3°). This is achieved by determining the skew symmetric matrix SE(t) = dRE(t)/dt * RE.T. I am obtaining incorrect results when utilizing a float instead of a Sympy expression, as shown in the following example:

from IPython.display import display
import sympy as sy

sy.init_printing()  # LaTeX like pretty printing for IPython


def mk_rotmatrix(alpha, coord_ax="x"):
    """ Rotation matrix around coordinate axis """
    ca, sa = sy.cos(alpha), sy.sin(alpha)
    if coord_ax == "x":
        return sy.Matrix([[1,  0,   0],
                          [0, ca, -sa],
                          [0, sa, +ca]])
    elif coord_ax == 'y':
        return sy.Matrix([[+ca, 0, sa],
                          [0,   1,  0],
                          [-sa, 0, ca]])
    elif coord_ax == 'z':
        return sy.Matrix([[ca, -sa, 0],
                          [sa, +ca, 0],
                          [0,    0, 1]])
    else:
        raise ValueError("Parameter coord_ax='" + coord_ax +
                         "' is not in ['x', 'y', 'z']!")


t, lat = sy.symbols("t, lat", real=True)  # time and latitude
omE = 7.292115e-5  # rad/s -- earth rotation rate (15.04107 °/h)
lat_sy = 48.232*sy.pi/180  # latitude in rad
lat_fl = float(lat_sy)  # latitude as float
print("\nlat_sy - lat_fl = {}".format((lat_sy - lat_fl).evalf()))

# earth rotation matrix at latitiude 48.232°:
RE = (mk_rotmatrix(omE*t, "z") * mk_rotmatrix(lat - sy.pi/2, "y"))
# substitute latitude with sympy and float value:
RE_sy, RE_fl = RE.subs(lat, lat_sy), RE.subs(lat, lat_fl)

# Angular velocity in world coordinates as skew symmetric matrix:
SE_sy = sy.simplify(RE_sy.diff(t) * RE_sy.T)
SE_fl = sy.simplify(RE_fl.diff(t) * RE_fl.T)

print("\nAngular velocity with Sympy latitude ({}):".format(lat_sy))
display(SE_sy)  # correct result
print("\nAngular velocity with float latitude ({}):".format(lat_fl))
display(SE_fl)  # incorrect result

The result is:

For the float latitude the result is totally wrong in spite of the difference of only -3e-17 to the Sympy value. It is not clear to me, why this happens. Numerically, this calculation does not seem to be problematic.

My question is, how to work around such deficits. Should I avoid mixing Sympy and float/Numpy data types? They are quite difficult to detect for more complex settings.

PS: The Sympy version is 0.7.6.

Answer 1

TL; DR

It is a bug. If you don't believe it, try this:

In [1]: from sympy import factor, Symbol

In [2]: factor(1e-20*Symbol('t')-7.292115e-5)
Out[2]: -2785579325.00000

---

Two years ago, the default value for the parameter tol in RealField.__init__ was changed from None to False in commit polys: Disabled automatic reduction to zero in RR and CC.
Later, tol was reverted back to None to fix a simplification issue, in commit Changed tol on Complex and Real field to None.
It seems the developers didn't expect this reversion would bring some other issue.

If you modify tol=None at RealField.__init__ in realfield.py, to tol=False, you will get the correct result for SE_fl.

Matrix([
[3.3881317890172e-21*sin(0.0001458423*t),                     -7.29211495242194e-5, 0],
[                    7.29211495242194e-5, -3.3881317890172e-21*sin(0.0001458423*t), 0],
[                                      0,                                        0, 0]])

---

The change of tol can explain why you've got a wrong result, but I don't thint it is the root of the issue.
IMHO, there is a deficiency in the polynomial factorization in SymPy. I'll illustrate this deficiency.
For convenience, let us do some preparation work.
Add the followings to your example.

from sympy import simplify, expand, S
from sympy.polys import factor
from sympy.polys.domains import QQ, RR, RealField
from sympy.polys.factortools import dup_convert
from sympy.polys.polytools import Poly
from sympy.polys.polytools import _symbolic_factor_list, _poly_from_expr
from sympy.polys.polyerrors import PolificationFailed
from sympy.polys import polyoptions as options
from sympy.simplify.fu import TR6

def new_opt():
    args = dict()
    options.allowed_flags(args, [])
    opt = options.build_options((), args)
    return opt
    
def my_symbolic_factor_list(base):
    opt = new_opt()
    try:
        poly, _ = _poly_from_expr(base, opt)
    except PolificationFailed as exc:
        print(exc)
        print(exc.expr)
    else:
        _coeff, _factors = poly.factor_list()
        print(poly)
        print(_coeff, _factors)
        return poly

We don't need to study the whole matrices. Let us focus on one element, element at row 1 and column 2. It has already shown the result is incorrect.

In [8]: elm_sy = (RE_sy.diff(t) * RE_sy.T)[1]

In [9]: elm_fl = (RE_fl.diff(t) * RE_fl.T)[1]

In [10]: elm_sy
Out[10]: -7.292115e-5*sin(0.267955555555556*pi)**2*sin(7.292115e-5*t)**2 - 7.292115e-5*sin(7.292115e
-5*t)**2*cos(0.267955555555556*pi)**2 - 7.292115e-5*cos(7.292115e-5*t)**2

In [11]: elm_fl
Out[11]: -7.292115e-5*sin(7.292115e-5*t)**2 - 7.292115e-5*cos(7.292115e-5*t)**2

In [12]: simplify(elm_sy)
Out[12]: -7.29211500000000e-5

In [13]: simplify(elm_fl)
Out[13]: -2785579325.00000

When we call simplify, in this case, it's almost equivalent to a combination of TR6 and factor.

In [15]: expr_sy = TR6(elm_sy)

In [16]: expr_fl = TR6(elm_fl)

In [17]: expr_fl
Out[17]: 1.35525271560688e-20*sin(7.292115e-5*t)**2 - 7.292115e-5

In [18]: factor(expr_fl)
Out[18]: -2785579325.00000

Now, we know wrong results would be produced during the invocation of factor().
Actually, factor is just a wrapper, the major work is done by _symbolic_factor_list.

In [20]: _symbolic_factor_list(expr_fl, opt, 'factor')
Out[20]: (-2785579325.00000, [])

Let us take a look at _symbolic_factor_list. The key part is:

        try:
            poly, _ = _poly_from_expr(base, opt)
        except PolificationFailed as exc:
            factors.append((exc.expr, exp))
        else:
            func = getattr(poly, method + '_list')

            _coeff, _factors = func()

We use the above my_symbolic_factor_list to simulate this procedure.

In [22]: expand(expr_sy)
Out[22]: -7.29211500000000e-5

In [23]: my_symbolic_factor_list(expr_sy)
can't construct a polynomial from -7.292115e-5*sin(0.267955555555556*pi)**2*sin(7.292115e-5*t)**2 -
7.292115e-5*(-sin(0.267955555555556*pi)**2 + 1)*sin(7.292115e-5*t)**2 + 7.292115e-5*sin(7.292115e-5*
t)**2 - 7.292115e-5
-7.29211500000000e-5

In [24]: my_symbolic_factor_list(S(1))
can't construct a polynomial from 1
1

In [25]: expr_fl
Out[25]: 1.35525271560688e-20*sin(7.292115e-5*t)**2 - 7.292115e-5    

In [26]: poly_fl = my_symbolic_factor_list(expr_fl)
Poly(-7.292115e-5, sin(7.292115e-5*t), domain='RR')
(-2785579325.00000, [])

By design, the constant polynomial should execute except PolificationFailed as exc: suite, while the other polynomials should execute else: suite.
expr_sy, which is a number after expand(), and 1 are both constant polynomials, thus PolificationFaileds were thrown.
poly_fl is -7.292115e-5 * sin(7.292115e-5*t) ** 0, namely, -7.292115e-5, a constant polynomial, whereas expr_fl is not. They were supposed to be the same polynomial, just different representation. Now they are not.
This is the deficiency I mentioned.

Where is the missing 1.35525271560688e-20*sin(7.292115e-5*t)**2?
Let us recall: tol was reverted back to None, which means automatic reduction to zero in RR is enabled again.
1.35525271560688e-20 was reduced to zero. Thus, poly_fl became a constant polynomial.
If tol is False, this won't happen.

In [31]: arg2 = expr_fl.args[1].args[0]

In [32]: arg2
Out[32]: 1.35525271560688e-20

In [33]: RR.from_sympy(arg2)
Out[33]: 0.0

In [34]: R = RealField(tol=False)

In [35]: R.from_sympy(arg2)
Out[35]: 1.35525271560688e-20

Now, we can explain why you've got -2785579325.0. In the else: suite, Poly.factor_list is called.
According to docs:

factor_list(f)[source]

Returns a list of irreducible factors of f.

poly_fl is supposed to be a non constant polynomial, but it is just a number. Thus, SymPy was tring to use a rational number to approximate poly_fl. The numerator is kept, while the denominator is discarded.

In [42]: poly_fl.factor_list()
Out[42]: (-2785579325.00000, [])

In [43]: dup_convert(poly_fl.coeffs(), RR, QQ)
Out[43]: [-2785579325/38199881995827]

In [44]: Poly([S(1.25)], t, domain='RR').factor_list()
Out[44]: (5.00000000000000, [])

In [45]: dup_convert(Poly([S(1.25)], t, domain='RR').coeffs(), RR, QQ)
Out[45]: [5/4]

In [46]: Poly((RE_fl.diff(t) * RE_fl.T)[3].args[0].args[0], t).factor_list()
Out[46]: (1767051195.00000, [])

I don't think we should blame mixing Sympy and float/Numpy data types. This problem is not caused by those pitfalls SymPy mentioned.
Even a very simple factorization can produce a counterintuitive result.

In [47]: factor(1e-20*t-1.2345e-5)
Out[47]: -539023891.000000

In [48]: factor(S(1e-20)*t-S(1.2345e-5))
Out[48]: -539023891.000000

So it is a bug. Just let the developers fix it.

Answer 2

I think this might be a bug in Sympy; when I run your script on my system (Ubuntu 14.04 64-bit, Python 2.7, Sympy 0.7.4.1), I get

lat_sy - lat_fl = -2.61291277482447e-17

Angular velocity with Sympy latitude (0.267955555555556*pi):
Matrix([
[          0, -7.292115e-5, 0],
[7.292115e-5,            0, 0],
[          0,            0, 0]])

Angular velocity with float latitude (0.841807204822):
Matrix([
[3.3881317890172e-21*sin(0.0001458423*t),                     -7.29211495242194e-5, 0],
[                    7.29211495242194e-5, -3.3881317890172e-21*sin(0.0001458423*t), 0],
[                                      0,                                        0, 0]])

which looks OK.

I'm not sure what to suggest: you could try an older version of Sympy than 0.7.6, or the latest revision from Github.

[In answer to comment] As to why the diagonals are non-zero, my first comment is that 3e-21/7e-5 is about 4e-17; IEEE754 64-bit ("float") numerical precision is around 2e-16. At 3e-21 rad/s one revolution will take 60 trillion years (about 2e21 s). Don't worry about it.

I'm not entirely sure what is happening here, but after adding this to your script

def matrix_product_element(a, b, i, j):
    v = a[3*i:3*i+3]
    w = b[j::3]
    summand_list = [v[k]*w[k]
                    for k in range(3)]

    print('element ({},{})'.format(i, j))
    print('  summand_list: {}'.format(summand_list))
    print('  sum(summand_list): {}'.format(sum(summand_list)))
    print('  sum(summand_list).simplify(): {}'.format(sum(summand_list)))

matrix_product_element(RE_fl.diff(t), RE_fl.T, 0, 0)
matrix_product_element(RE_fl.diff(t), RE_fl.T, 1, 0)
matrix_product_element(RE_fl.diff(t), RE_fl.T, 2, 0)

sumlist=[sy.Float(-4.05652668591092e-5,15), sy.Float(7.292115e-5,15), sy.Float(-3.23558831408908e-5,14)]
display(sumlist)
display(sum(sumlist))

I get

element (0,0)
  summand_list: [-4.05652668591092e-5*sin(7.292115e-5*t)*cos(7.292115e-5*t), 7.292115e-5*sin(7.292115e-5*t)*cos(7.292115e-5*t), -3.23558831408908e-5*sin(7.292115e-5*t)*cos(7.292115e-5*t)]
  sum(summand_list): 6.7762635780344e-21*sin(7.292115e-5*t)*cos(7.292115e-5*t)
  sum(summand_list).simplify(): 6.7762635780344e-21*sin(7.292115e-5*t)*cos(7.292115e-5*t)
element (1,0)
  summand_list: [4.05652668591092e-5*cos(7.292115e-5*t)**2, 7.292115e-5*sin(7.292115e-5*t)**2, 3.23558831408908e-5*cos(7.292115e-5*t)**2]
  sum(summand_list): 7.292115e-5*sin(7.292115e-5*t)**2 + 7.292115e-5*cos(7.292115e-5*t)**2
  sum(summand_list).simplify(): 7.292115e-5*sin(7.292115e-5*t)**2 + 7.292115e-5*cos(7.292115e-5*t)**2
element (2,0)
  summand_list: [0, 0, 0]
  sum(summand_list): 0
  sum(summand_list).simplify(): 0
[-4.05652668591092e-5, 7.29211500000000e-5, -3.2355883140891e-5]
6.77626357803440e-21

The coefficients of the first summation should sum to zero, but don't. I've managed to sort-of fake this effect in the last few lines by recreating the coefficients with lower precision (this was just luck, and probably not that signicant). It's "sort-of" since the third value in the list (-3.2355883140891e-5) doesn't match the coefficient in the summand list (-3.23558831408908e-5), which is given to 15 places.

The Sympy docs discuss these sorts of issue here http://docs.sympy.org/dev/gotchas.html#evaluating-expressions-with-floats-and-rationals , with some suggestions on how to mitigate the problem. Here's a straightforward variation on your code, deferring substitution of floats right to the end:

# encoding:utf-8
from IPython.display import display
import sympy as sy

sy.init_printing()  # LaTeX like pretty printing for IPython


def mk_rotmatrix(alpha, coord_ax="x"):
    """ Rotation matrix around coordinate axis """
    ca, sa = sy.cos(alpha), sy.sin(alpha)
    if coord_ax == "x":
        return sy.Matrix([[1,  0,   0],
                          [0, ca, -sa],
                          [0, sa, +ca]])
    elif coord_ax == 'y':
        return sy.Matrix([[+ca, 0, sa],
                          [0,   1,  0],
                          [-sa, 0, ca]])
    elif coord_ax == 'z':
        return sy.Matrix([[ca, -sa, 0],
                          [sa, +ca, 0],
                          [0,    0, 1]])
    else:
        raise ValueError("Parameter coord_ax='" + coord_ax +
                         "' is not in ['x', 'y', 'z']!")


# time [s], latitude [rad], earth rate [rad/s]
t, lat, omE = sy.symbols("t, lat, omE", real=True)

RE = (mk_rotmatrix(omE*t, "z") * mk_rotmatrix(lat - sy.pi/2, "y"))

SE = sy.simplify(RE.diff(t) * RE.T)

display(SE)
display(SE.subs({lat: 48.232*sy.pi/180, omE: 7.292115e-5}))

This gives:

Matrix([
[  0, -omE, 0],
[omE,    0, 0],
[  0,    0, 0]])
Matrix([
[          0, -7.292115e-5, 0],
[7.292115e-5,            0, 0],
[          0,            0, 0]])

I prefer this regardless of numerical advantages, since one may learn something from the form of the symbolic solution.