矩阵指数

---

此示例说明 19 种矩阵指数计算方法中的 3 种。

有关矩阵指数计算的背景信息，请参阅：

Moler, Cleve, and Charles Van Loan.“Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later.”SIAM Review 45, no. 1 (January 2003):3–49. https://doi.org/10.1137/S00361445024180 (opens new window).

首先创建矩阵 A。

using TyMath
using TyBase
A = [0 1 2; 0.5 0 1; 2 1 0]

A =
3×3 Matrix{Float64}:
 0.0  1.0  2.0
 0.5  0.0  1.0
 2.0  1.0  0.0

导入 expmdemo1.jl、expmdemo2.jl、expmdemo3.jl 文件。

include("expmdemo1.jl")
include("expmdemo2.jl")
include("expmdemo3.jl")

Asave = A;

方法 1：加权平方

expmdemo1 是以下著作中算法 11.3.1 的实现：

Golub, Gene H. and Charles Van Loan.Matrix Computations, 3rd edition.Baltimore, MD:Johns Hopkins University Press, 1996.

f, e = log2(opnorm(A, Inf), out = true)
s = max(0, e + 1)
A = A / 2^s
X = A
c = 1 / 2
E = eye(size(A)) + c * A
D = eye(size(A)) - c * A
q = 6
p = 1
for k = 2:q
  global c = c * (q - k + 1) / (k * (2 * q - k + 1))
  global X = A * X
  cX = c * X
  global E = E + cX
  if Bool(p)
    global D = D + cX
  else
    global D = D - cX
  end
  global p = !Bool(p)
end
E = D \ E

for k = 1:s
  global E = E * E
end

E1 = E

E1 =
3×3 Matrix{Float64}:
 5.30908  4.0012   5.57784
 2.80879  2.88452  3.19301
 5.17375  4.0012   5.71318

方法 2：泰勒级数

expmdemo2 使用矩阵指数的经典定义，表示为幂级数

.

是与 A 具有相同维度的单位矩阵。作为一种实用的数值方法，如果 opnorm(A) 太大，此方法将很慢且不准确。

A = Asave

E = zeros(size(A))
F = eye(size(A))
k = 1

while opnorm(E + F - E, 1) > 0
  global E = E + F
  global F = A * F / k
  global k = k + 1
end

E2 = E

E2 =
3×3 Matrix{Float64}:
 5.30908  4.0012   5.57784
 2.80879  2.88452  3.19301
 5.17375  4.0012   5.71318

方法 3：特征值和特征向量

expmdemo3 假定矩阵包含一组完整的特征向量 V，使得 A=VD。矩阵指数可以通过对特征值的对角矩阵求幂来计算：

.

作为一种实际的数值方法，准确性由特征向量矩阵的条件确定。

A = Asave

D, V = eigen(A)
E = V * diagm(exp.(D)) / V

E3 = E

E3 =
3×3 Matrix{Float64}:
 5.30908  4.0012   5.57784
 2.80879  2.88452  3.19301
 5.17375  4.0012   5.71318

比较结果

对于此示例中的矩阵，所有三种方法都同样有效。

E = exp(Asave)
err1 = E - E1

err1 =
3×3 Matrix{Float64}:
 -1.77636e-15  -3.55271e-15  -6.21725e-15
 -1.77636e-15  -1.33227e-15  -3.10862e-15
 -6.21725e-15  -5.32907e-15  -9.76996e-15

err2 = E - E2

err2 =
3×3 Matrix{Float64}:
 -5.32907e-15  -5.32907e-15  -8.88178e-15
 -3.10862e-15  -2.66454e-15  -3.55271e-15
 -4.44089e-15  -5.32907e-15  -6.21725e-15

err3 = E - E3

err3 =
3×3 Matrix{Float64}:
 -1.24345e-14  -8.88178e-15  -1.5099e-14
 -9.76996e-15  -7.10543e-15  -1.15463e-14
 -1.68754e-14  -1.06581e-14  -1.68754e-14

泰勒级数失败

对于某些矩阵，泰勒级数中的项在变为零之前变得非常大。因此，expmdemo2 失败。

A = [-147 72; -192 93]
E1 = expmdemo1(A)

E1 =
2×2 Matrix{Float64}:
 -0.0995741  0.0746806
 -0.199148   0.149361

E2 = expmdemo2(A)

E2 =
2×2 Matrix{Float64}:
 -1.19853e6  -5.90782e5
 -2.74382e6  -2.04418e6

E3 = expmdemo3(A)

E3 =
2×2 Matrix{Float64}:
 -0.0995741  0.0746806
 -0.199148   0.149361

特征值和特征向量失败

以下是不包含一组完整的特征向量的矩阵。因此，expmdemo3 失败。

A = [-1 1; 0 -1]
E1 = expmdemo1(A)

E1 =
2×2 Matrix{Float64}:
 0.367879  0.367879
 0.0       0.367879

E2 = expmdemo2(A)

E2 =
2×2 Matrix{Float64}:
 0.367879  0.367879
 0.0       0.367879

E3 = expmdemo3(A)

E3 =
2×2 Matrix{Float64}:
 0.367879  0.0
 0.0       0.367879