Sympy 和 Numpy 无法对使用 Sympy 创建的矩阵求逆

Question

我有一个方程，它有一个使用 Sympy 包命名的符号变量 a。使用这个变量我创建了一个 6 * 6 矩阵：

array([[1, exp(-0.04*a**2), exp(-0.16*a**2), exp(-0.36*a**2),
        exp(-0.64*a**2), exp(-1.0*a**2)],
       [exp(-0.04*a**2), 1, exp(-0.04*a**2), exp(-0.16*a**2),
        exp(-0.36*a**2), exp(-0.64*a**2)],
       [exp(-0.16*a**2), exp(-0.04*a**2), 1, exp(-0.04*a**2),
        exp(-0.16*a**2), exp(-0.36*a**2)],
       [exp(-0.36*a**2), exp(-0.16*a**2), exp(-0.04*a**2), 1,
        exp(-0.04*a**2), exp(-0.16*a**2)],
       [exp(-0.64*a**2), exp(-0.36*a**2), exp(-0.16*a**2),
        exp(-0.04*a**2), 1, exp(-0.04*a**2)],
       [exp(-1.0*a**2), exp(-0.64*a**2), exp(-0.36*a**2),
        exp(-0.16*a**2), exp(-0.04*a**2), 1]], dtype=object)

上面的 a 变量是一个 sympy 符号。现在我想反转这个矩阵。使用 numpy 我得到以下错误：

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-117-3c958037db81> in <module>()
----> 1 inv(np.array(final_result))

D:\datascience\lib\site-packages\numpy\linalg\linalg.py in inv(a)
    526     signature = 'D->D' if isComplexType(t) else 'd->d'
    527     extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
--> 528     ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
    529     return wrap(ainv.astype(result_t, copy=False))
    530 

TypeError: No loop matching the specified signature and casting
was found for ufunc inv

并使用 Sympy 自己的语法来反转这个矩阵，即 M ** -1，需要太长时间并且不会计算结果。 有没有办法使用numpy来反转这个矩阵？还是纯python语法？

Answer 1

如果您正在寻找符号逆，则需要使用 SymPy。 NumPy 只计算数值逆，它使用的算法不适用于符号矩阵。

看起来 SymPy 使用的默认算法对这个矩阵很慢，但是您可以通过将 method 关键字传递给 Matrix 的 inv 方法来尝试替代算法。可用选项为'GE'（默认）、'LU' 和'ADJ'。

M.subs(exp(-a**2/25), b) 也可以帮助计算。这样得到的矩阵只有b 的幂（你需要使用有理数而不是浮点数才能工作）。我能够得到一个逆

a, b = symbols('a b')
M = Matrix([[1, exp(-a**2/25), exp(-4*a**2/25), exp(-9*a**2/25), exp(-16*a**2/25), exp(-a**2)], [exp(-a**2/25), 1, exp(-a**2/25), exp(-4*a**2/25), exp(-9*a**2/25), exp(-16*a**2/25)], [exp(-4*a**2/25), exp(-a**2/25), 1, exp(-a**2/25), exp(-4*a**2/25), exp(-9*a**2/25)], [exp(-9*a**2/25), exp(-4*a**2/25), exp(-a**2/25), 1, exp(-a**2/25), exp(-4*a**2/25)], [exp(-16*a**2/25), exp(-9*a**2/25), exp(-4*a**2/25), exp(-a**2/25), 1, exp(-a**2/25)], [exp(-a**2), exp(-16*a**2/25), exp(-9*a**2/25), exp(-4*a**2/25), exp(-a**2/25), 1]])
M2 = M.subs(exp(-a**2/25), b)
M2inv = simplify(M2.inv('ADJ')) # takes a little time to compute, about 1:30 on my machine
Minv = M2inv.subs(b, exp(-a**2/25))

在一天结束时，您可能需要考虑是否需要矩阵的符号逆。用数值替换a 然后用 NumPy 计算数值逆是否足够？你真的需要逆吗？通常计算矩阵的显式逆对于数值计算来说不是一个好主意。通常最好使用numpy.linalg.solve 之类的东西。

符号逆通常是非常复杂的表达式，本质上难以计算。数值逆在数值上非常不稳定，并且几乎从来都不是一个好主意来计算来解决您要解决的任何问题。

Answer 2

NumPy 是一个数值包，不能与 SymPy 符号一起使用。 如果我理解正确，您使用 a 作为 SymPy 符号，而没有为其分配具体值。这使得计算逆非常困难。对于a 的某些分配，矩阵甚至可能是不可逆的。恐怕 SymPy 的矩阵求逆不能胜任这样的任务。可能会有所帮助的方法是将数字表示为分数，因为实数很难用符号表示。

如果您尝试对较小的矩阵进行同样的操作，您会注意到结果很快就会变成一堆复杂的符号。对于您需要的a 的每个值，对矩阵进行数值反转（使用 NumPy）可能会更容易。

另一个尝试的技巧是引入一个辅助符号来获得更简单的表达式，然后再替换回来。例如。首先是M.subs({a**2:-log(b)})。得到的逆矩阵可能更容易简化。然后M.subs({b: exp(-a**2)}) 调查原始公式会发生什么。仅使用此替换计算原始矩阵的行列式会导致表达式长度超过 1000 个字符，这看起来不容易简化。完整的逆矩阵是一个易于管理的表达式的希望仍然很小。