使用（numpy）浮动时Sympy的结果不正确

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:

我试图从时间相关的旋转矩阵RE（t）（即纬度48.3°的地球自转）计算速度张量。这是通过确定偏斜对称矩阵SE（t）= dRE（t）/ dt * RE.T来实现的。使用float而不是Sympy表达式时，我得到的结果不正确，如下例所示：

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.

对于浮动纬度，尽管只有-3e-17与Sympy值的差异，但结果完全错误。我不清楚为什么会这样。从数字上看，这种计算似乎没有问题。

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.

我的问题是，如何解决这些缺陷。我应该避免混合Sympy和float / Numpy数据类型吗？对于更复杂的设置，很难检测到它们。

PS: The Sympy version is 0.7.6.

PS：Sympy版本是0.7.6。

2 个解决方案

#1

TL; DR

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.

两年前，RealField .__ init__中参数tol的默认值在commit polys中从None更改为False：在RR和CC中禁用自动减少为零。后来，tol被恢复为None以修复简化问题，在将Complex和Real字段的提交更改为无。似乎开发人员没想到这种回归会带来一些其他问题。

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

如果你在Realfield.py中的RealField .__ init__修改tol = None，到tol = False，你将得到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.

tol的变化可以解释为什么你得到了错误的结果，但我不认为它是问题的根源。恕我直言，SymPy中的多项式因子分解存在缺陷。我将说明这种不足。为方便起见，让我们做一些准备工作。将以下内容添加到您的示例中。

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.

我们不需要研究整个矩阵。让我们关注一个元素，第1行和第2列的元素。它已经显示结果不正确。

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.

当我们称之为简化时，在这种情况下，它几乎相当于TR6和因子的组合。

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.

现在，我们知道在调用factor（）期间会产生错误的结果。实际上，因子只是一个包装器，主要工作是由_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:

我们来看看_symbolic_factor_list。关键部分是：

        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.

我们使用上面的my_symbolic_factor_list来模拟这个过程。

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.

按照设计，常量多项式应该执行，除了PolificationFailed为exc：suite，而其他多项式应该执行else：suite。 expr_sy，它是expand（）之后的数字，1是常量多项式，因此抛出了PolificationFailed。 poly_fl是-7.292115e-5 * sin（7.292115e-5 * t）** 0，即-7.292115e-5，常数多项式，而expr_fl不是。它们应该是相同的多项式，只是不同的表示。现在他们不是。这是我提到的不足之处。

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.

遗失的地方1.35525271560688e-20 * sin（7.292115e-5 * t）** 2？让我们回想一下：tol被恢复为None，这意味着RR再次启用自动减少到零。 1.35525271560688e-20减少到零。因此，poly_fl变为常数多项式。如果tol为False，则不会发生这种情况。

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:

现在，我们可以解释为什么你有-2785579325.0。在else：套件中，调用Poly.factor_list。根据文件：

factor_list(f)[source]

factor_list（F）[来源]

Returns a list of irreducible factors of f.

返回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.

poly_fl应该是一个非常数多项式，但它只是一个数。因此，SymPy试图使用有理数来逼近poly_fl。分子被保留，而分母被丢弃。

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.

我不认为我们应该责怪混合Sympy和float / Numpy数据类型。这个问题不是由SymPy提到的那些陷阱造成的。即使是非常简单的因子分解也会产生违反直觉的结果。

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.

所以这是一个错误。让开发人员修复它。

#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

我认为这可能是Sympy的一个错误;当我在我的系统上运行你的脚本（Ubuntu 14.04 64位，Python 2.7，Sympy 0.7.4.1）时，我得到了

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.

我不确定该建议：您可以尝试使用Sympy的旧版本而不是0.7.6，或者是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.

[回答评论]关于为什么对角线不为零，我的第一个评论是3e-21 / 7e-5约为4e-17; IEEE754 64位（“浮点”）数值精度约为2e-16。在3e-21rad / s时，一次旋转需要60万亿年（约2e21秒）。别担心。

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.

第一个求和的系数应该总和为零，但不是。通过重新创建精度较低的系数，我设法在最后几行中对这种效果进行了排序（这只是运气，而且可能不是那么重要）。它是“排序”，因为列表中的第三个值（-3.2355883140891e-5）与summand列表（-3.23558831408908e-5）中的系数不匹配，该系数被赋予15个位置。

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:

Sympy文档在http://docs.sympy.org/dev/gotchas.html#evaluating-expressions-with-floats-and-rationals中讨论了这些问题，并就如何缓解问题提出了一些建议。这是您的代码的直接变化，将浮动的替换推迟到最后：

# 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.

我更喜欢这个，不管数字优势，因为人们可以从符号解决方案的形式中学到一些东西。