Python：方程式的矩阵指数

编程问答 2022-06-10

问题描述

我正在尝试在Python中实现以下等式。

enter image description here

这是一个方程，用于计算给定矩阵A和标量x的矩阵指数。 当我将其与scipy中的Python expm进行比较时，我的代码似乎不起作用。

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)

我已经定义了矩阵matrix = np.array([[-5,-9]])，并使用scipy的expm函数获得了输出：

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

但是我对等式的实现给了我

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

我已经盯着我的代码好几个小时了，但是我最近才开始学习Python，所以我的技能非常有限。

解决方法

这是因为精度下降。事实证明，执行矩阵指数运算需要太多项才能收敛（在这种情况下，大约为35），并且矩阵M ^ 35已经分解了整数。使用相同的算法，让我们看看Julia如何做到：

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,2}:
 0.365246  0.353079  0.28076
 0.363873  0.353114  0.283013
 0.367464  0.358922  0.305631

从上面可以看到，当矩阵位于Int64与Int128中时，M^20的区别很大。确实，它在不引发异常的情况下静默地溢出了整数。如果对这些条件加总，那也是您得到错误答案的相同原因。

不幸的是，numpy没有像Julia那样的int128类型。但是我们确实有float128。因此，让我们修改代码以改为使用float128：

>>> 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,[ 4,-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)

同样，在这里，您会看到当数据类型不同时，矩阵M的差将提高到20的幂。使用float会损失一些精度，但至少不会针对此特定矩阵，您不会过早地得到答案并获得正确的答案。

经验教训：如果scipy为您提供一个功能，请不要实现您自己的功能。人们在图书馆中实施它是有原因的。

matrixnumerical-analysispython