python - 使用 sympy 和 numpy 获取倒置浮点矩阵

Question

我正在尝试在 python 上实现一些算法。为了文档和对流程细节的清晰理解，我使用 sympy。事实证明，它在计算倒置浮点矩阵时失败了。

所以我明白了

TypeError                                 Traceback (most recent call last)
<ipython-input-20-c2193b2ae217> in <module>()
     10 np.linalg.inv(xx)
     11 symInv = lambdify(X0,X0Inv)
---> 12 symInv(xx)

/opt/anaconda3/lib/python3.6/site-packages/numpy/__init__.py in <lambda>(X0)

TypeError: ufunc 'bitwise_xor' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

如果矩阵是整数，它可以正常工作：

import numpy as np
from sympy import *
init_printing()


X0 = MatrixSymbol('X0',2,2)
xx = np.random.rand(4,4) 
#xx = np.random.randint(10,size=(4,4)) # this line makes it workable
X0Inv = X0**-1
np.linalg.inv(xx)
symInv = lambdify(X0,X0Inv)
symInv(xx)

Link to a live version of the code

如果有人知道任何解决方法，如果您能分享，那就太好了。提前致谢。

更新。正如@hpaulj 和@tel 所指出的，问题是lambdify 如何将** 转换为矩阵符号的numpy 代码：由于某种原因，它试图对元素进行异或。我将尝试找到一种简单的方法来改变这种行为。任何帮助/提示表示赞赏。

Answer 1

正如 hpaulj 指出的那样，由于某种原因，该错误似乎源于 lambdify 中发生的 ** 到 ^ 的转换。

您可以通过使用np.power 而不是** 来修复您遇到的错误：

import numpy as np
from sympy import MatrixSymbol, lambdify

X0 = MatrixSymbol('X0',2,2)
xx = np.random.rand(4,4)
X0Inv = np.power(X0, -1)
symInv = lambdify(X0,X0Inv)

print('matrix xx')
print(xx, end='\n\n')
print('result of symInv(xx)')
print(symInv(xx), end='\n\n')

输出：

matrix xx
[[0.4514882  0.84588859 0.02431252 0.25468078]
 [0.46767727 0.85748153 0.51207567 0.59636962]
 [0.84557537 0.38459205 0.76814414 0.96624407]
 [0.0933803  0.43467119 0.77823338 0.58770188]]

result of symInv(xx)
[[2.214897321138516, 1.1821887747951494], [2.1382266426713077, 1.1662058776397513]]

但是，正如您设置的那样，symInv 不会产生矩阵求逆，而只会对xx 中的每个值进行逐元素取幂。换句话说，symInv(xx)[i,j] == xx[i,j]**-1。这段代码显示了逐元素取幂和真正的逆运算之间的区别。

print('result of xx**-1')
print(xx**-1, end='\n\n')
print('result of np.linalg.inv(xx)')
print(np.linalg.inv(xx))

输出：

result of xx**-1
[[ 2.21489732  1.18218877 41.13107402  3.92648394]
 [ 2.13822664  1.16620588  1.95283638  1.67681243]
 [ 1.18262669  2.60015778  1.301839    1.0349352 ]
 [10.7088969   2.30058954  1.28496159  1.70154295]]

result of np.linalg.inv(xx)
[[-118.7558445   171.37619558  -20.37188041  -88.94733652]
 [  -0.56274492    2.49107626   -1.00812489   -0.62648633]
 [-160.35674704  230.3266324   -28.87548299 -116.75862026]
 [ 231.62940572 -334.07044947   42.21936405  170.90926978]]

编辑：解决方法

我 95% 确定您遇到的是 Sympy 代码中的错误。似乎 X0^-1 在某些时候是 Sympy Matrix 对象的有效语法，但不再是。然而，似乎有人忘了告诉维护lambdify 代码的人，因为它仍然将每个矩阵求幂转换为胡萝卜^ 语法。

因此，您应该在Sympy github 上提交问题。只需发布您的代码及其产生的错误，并询问这是否是预期的行为。与此同时，这里有一个肮脏的黑客来解决这个问题：

import numpy as np
from sympy import MatrixSymbol, lambdify

class XormulArray(np.ndarray):
    def __new__(cls, input_array):
        return np.asarray(input_array).view(cls)

    def __xor__(self, other):
        return np.linalg.matrix_power(self, other)

X0 = MatrixSymbol('X0',2,2)
xx = np.random.rand(4,4)
X0Inv = X0.inv()
symInv = lambdify(X0,X0Inv,'numpy')

print('result of symInv(XormulArray(xx))')
print(symInv(XormulArray(xx)), end='\n\n')

print('result of np.linalg.inv(xx)')
print(np.linalg.inv(xx))

输出：

result of symInv(XormulArray(xx))
[[ 3.50382881 -3.84573344  3.29173896 -2.01224981]
 [-1.88719742  1.86688465  0.3277883   0.0319487 ]
 [-3.77627792  4.30823019 -5.53247103  5.53412775]
 [ 3.89620805 -3.30073088  4.27921307 -4.68944191]]

result of np.linalg.inv(xx)
[[ 3.50382881 -3.84573344  3.29173896 -2.01224981]
 [-1.88719742  1.86688465  0.3277883   0.0319487 ]
 [-3.77627792  4.30823019 -5.53247103  5.53412775]
 [ 3.89620805 -3.30073088  4.27921307 -4.68944191]]

基本上，您必须在将所有数组传递给symInv 之前将它们转换为瘦包装类型XormulArray。由于多种原因（包括它显然打破了您对 X0 施加的 (2,2) 形状限制这一事实），这种 hack 并不是最佳实践，但在修复 Sympy 代码库之前，它可能是您能做的最好的事情.