Python: Matrix exponential from equation

Question

I'm trying to implement the following equation in Python.

It's an equation for computing the matrix exponential for a given matrix A and scalar x. My code doesn't seem to work, when I'm comparing it to the Python expm from scipy.

import math
import numpy as np
from scipy.linalg import expm

# Scalar x (will later on be for user input)

x = 1

matrix = np.array([[-5, 2, 3], [2, -6, 4], [4, 5, -9]])

# Using scipy to compute the matrix exponential (for comparison)

B = expm(matrix)
print(B)
    
# Defining the equation

def equation(n):
    y = ((pow(x, n) * np.linalg.matrix_power(matrix, n)) / int(math.factorial(n)))
    return y


# Summing the equation with finite iterations

result = sum([equation(n) for n in range(0, 1000)])
print(result)

I have defined the matrix matrix = np.array([[-5, 2, 3], [2, -6, 4], [4, 5, -9]]) and with the expm function from scipy I get the output:

 [[0.3659571  0.35453832 0.27950458]
 [0.36527461 0.35510049 0.27962489]
 [0.36551524 0.35489926 0.27958549]]

But my implementation of the equations gives me:

[[282.7927229097503 439.9138578271309 2167.1107527813792]
 [548.8430305150805 -1876.4510112837029 1328.9683527937962]
 [1753.0719360816013 3838.501983853133 -5590.574633487889]]

I've been staring at my code for hours, but I just recently picked up Python so my skills are very limited.

Answer 1

This is because of lost of precision. It turns out, doing matrix exponential needs too many terms to converge (around 35 in this case) and the matrix M^35 already exploded your integer. Using the same algorithm, let's see how Julia can do it:

julia> M = [-5 2 3; 2 -6 4; 4 5 -9]
3×3 Array{Int64,2}:
 -5   2   3
  2  -6   4
  4   5  -9

julia> exp(M)
3×3 Array{Float64,2}:
 0.365957  0.354538  0.279505
 0.365275  0.3551    0.279625
 0.365515  0.354899  0.279585

julia> term = (n) -> (M^n)/factorial(big(n))
#1 (generic function with 1 method)

julia> sum(term, 0:40)
3×3 Array{BigFloat,2}:
  282.793    439.914   2167.11
  548.843  -1876.45    1328.97
 1753.07    3838.5    -5590.57

julia> M^20
3×3 Array{Int64,2}:
 8757855768227185495   5428672870161680643   4260215435320685478
 2846510725988846806  -6309877790968251876   3463367064979405070
 1252813985306285990   3038127759137839419  -4290941744444125409

julia> M = Matrix{Int128}(M)
3×3 Array{Int128,2}:
 -5   2   3
  2  -6   4
  4   5  -9

julia> M^20
3×3 Array{Int128,2}:
   691287386495480595287   1259807269882411190531  -1951094656377891785818
  1423245804401624321238   2594681036602078525980  -4017926841003702847218
 -2710418564849997801562  -4940689283995021993669   7651107848845019795231

julia> sum(term, 0:40)
3×3 Array{BigFloat,2}:
 0.365246  0.353079  0.28076
 0.363873  0.353114  0.283013
 0.367464  0.358922  0.305631

From the above, you see the big difference for M^20 when the matrix is in Int64 vs Int128. Indeed, it silently overflow the integer without throwing exceptions. That's the same reason you got the answer wrong if you sum the terms up.

Unfortunately, numpy does not have int128 type as Julia. But we do have float128. So let's modify your code to use float128 instead:

>>> from scipy.linalg import expm
>>> import numpy as np
>>> import math
>>> M = np.array([[-5, 2, 3], [2, -6, 4], [4, 5, -9]])
>>> M
array([[-5,  2,  3],
       [ 2, -6,  4],
       [ 4,  5, -9]])
>>> expm(M)
array([[0.3659571 , 0.35453832, 0.27950458],
       [0.36527461, 0.35510049, 0.27962489],
       [0.36551524, 0.35489926, 0.27958549]])
>>> np.linalg.matrix_power(M, 20)
array([[ 8757855768227185495,  5428672870161680643,  4260215435320685478],
       [ 2846510725988846806, -6309877790968251876,  3463367064979405070],
       [ 1252813985306285990,  3038127759137839419, -4290941744444125409]])
>>> term = lambda n: np.linalg.matrix_power(M, n)/float(math.factorial(n))
>>> sum([term(n) for n in range(40)])
array([[  282.79272291,   439.91385783,  2167.11075278],
       [  548.84303052, -1876.45101128,  1328.96835279],
       [ 1753.07193608,  3838.50198385, -5590.57463349]])
>>> M = M.astype('float128')
>>> M
array([[-5.,  2.,  3.],
       [ 2., -6.,  4.],
       [ 4.,  5., -9.]], dtype=float128)
>>> np.linalg.matrix_power(M, 20)
array([[ 6.91287386e+20,  1.25980727e+21, -1.95109466e+21],
       [ 1.42324580e+21,  2.59468104e+21, -4.01792684e+21],
       [-2.71041856e+21, -4.94068928e+21,  7.65110785e+21]],
      dtype=float128)
>>> sum([term(n) for n in range(40)])
array([[0.36595003, 0.35452543, 0.27952454],
       [0.36526005, 0.35507395, 0.279666  ],
       [0.36554297, 0.35494981, 0.27950722]], dtype=float128)

Same here, you see the difference of matrix M raise to 20th power when the data type are different. You lost some precision by using float, but you don't overflow to early and get your answer right, at least for this particular matrix.

Lesson learned: Do not implement your own function if scipy provide one for you. There is a reason people implement it in a library.