时间戳

首页 > 解决方案 > Python：来自方程的矩阵指数

python - Python：来自方程的矩阵指数

问题描述

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

在此处输入图像描述

这是一个计算给定矩阵 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, 2, 3], [2, -6, 4], [4, 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，所以我的技能非常有限。

标签： pythonmatrixnumerical-analysis

解决方案

这是因为失去了精度。事实证明，做矩阵指数需要太多项来收敛（在这种情况下大约是 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, 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

从上面，您可以看到M^20矩阵在 Int64 与 Int128 中的巨大差异。事实上，它默默地溢出整数而不抛出异常。如果你总结这些术语，这就是你得到错误答案的原因。

不幸的是，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,  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)

M同样在这里，当数据类型不同时，您会看到矩阵的差异提高到 20 次方。使用浮点数会丢失一些精度，但至少对于这个特定的矩阵，您不会过早溢出并得到正确的答案。

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

推荐阅读