不能在循环中添加矩阵超过 20 次

Question

我试图在 for 循环中创建 2x2 矩阵的总和，但是当我循环总和超过 21 次时（当我的 n > 20 时，如下所示）它给我以下错误消息：

TypeError: ufunc 'add' output (typecode 'O') could not be coerced to provided output parameter (typecode 'd') according to the casting rule ''same_kind''

这是我的代码：

k = 2
n = 21

A2 = np.matrix('0.5 -0.5; 0.5 0.5')
SumA2 = np.zeros((k,k))

for i in range(0, n+1):
    SumA2 += np.linalg.matrix_power(A2, i)/np.math.factorial(i)

print(A2)
print("\n", SumA2)

我怀疑它与阶乘变得太大有关，但这真的是个问题吗？在 Matlab 中，我可以毫无问题地循环它 1000 次。

Answer 1

在 21，它将数组类型切换为对象：

In [776]: np.linalg.matrix_power(A2,20)/np.math.factorial(20)
Out[776]: 
matrix([[-4.01398205e-22,  0.00000000e+00],
        [ 0.00000000e+00, -4.01398205e-22]])
In [777]: np.linalg.matrix_power(A2,21)/np.math.factorial(21)
Out[777]: 
matrix([[-9.557100128609015e-24, 9.557100128609015e-24],
        [-9.557100128609015e-24, -9.557100128609015e-24]], dtype=object)

更具体地说，切换的是 factorial：

In [778]: np.array(np.math.factorial(20))
Out[778]: array(2432902008176640000)
In [779]: np.array(np.math.factorial(21))
Out[779]: array(51090942171709440000, dtype=object)

Python3 对 factorial 使用整数。这些可以是任何长度。但此时值变得太大而无法用 np.int64 表示。所以它切换到使用一个包含长 Python 整数的对象 dtype 数组。该开关传播到 power 计算。

尝试将此数组转换为与 SumA2 兼容的 dtype 时出现错误。

In [782]: SumA2 = np.zeros((k,k))
In [783]: SumA2 += Out[777]
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-783-53cbd27f9514> in <module>()
----> 1 SumA2 += Out[777]

TypeError: ufunc 'add' output (typecode 'O') could not be coerced to provided output parameter (typecode 'd') according to the casting rule ''same_kind''
In [784]: SumA2 = np.zeros((k,k), object)
In [785]: SumA2 += Out[777]
In [786]: SumA2
Out[786]: 
array([[-9.557100128609015e-24, 9.557100128609015e-24],
       [-9.557100128609015e-24, -9.557100128609015e-24]], dtype=object)

在 170 时开始出现将整数转换为浮点数的问题

先做一个1/factorial(...)似乎有帮助。将 A2 的 dtype 更改为更高精度的浮点数可能会有所帮助：

In [812]: np.linalg.matrix_power(A2.astype('float128'),171)*(1/np.math.factorial(171))
Out[812]: 
matrix([[-1.04145922e-335, -1.04145922e-335],
        [ 1.04145922e-335, -1.04145922e-335]], dtype=float128)

对于 2x2 矩阵，这确实没有特别利用 numpy。使用列表和 'raw' Python 数字几乎可以很容易地计算出重复次数。但即使是那些也不是为无限精确的数学而设计的。整数可以很长，但我认为 Python 浮点数没有那么灵活。

不能在循环中添加矩阵超过 20 次

Cannot add matrices more than 20 times in loop

python

numpy

matrix

factorial