Невозможно добавить матрицы более 20 раз в цикле

Question

Я пытаюсь создать сумму матриц 2x2 в цикле for, но когда я повторяю суммирование более 21 раза (при моем n> 20, как показано ниже), выдается следующее сообщение об ошибке:

Ошибка типа: ошибка вывода ufunc 'add' (typecode 'O') не может быть принудительно приведена к предоставленному выходному параметру (typecode 'd') в соответствии с правилом приведения типов '' 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.Поэтому он переключается на использование массива dtype объекта, который содержит длинное целое число Python.Этот переключатель распространяется на вычисление power.

Ошибка возникает, когда он пытается преобразовать этот массив в d-тип, совместимый с SumA2.

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(...) первым кажется полезным.И может помочь изменение типа d A2 на более высокую точность с плавающей точкой:

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.Повторную мощность можно почти так же легко вычислить с помощью списков и «сырых» чисел Python.Но даже они не предназначены для математики с бесконечной точностью.Целые числа могут быть длинными, но я не думаю, что Python плавает настолько гибко.