使用Numpy计算矩阵指数时两种方法结果不一致的问题

问题背景 #

对于2x2矩阵A，若T是由其特征向量组成的可逆矩阵，D是对角元素为A特征值的对角矩阵，则有(A = TDT^{{-1})。根据矩阵幂运算规则和指数函数的幂级数展开，可推导出矩阵指数公式：(e}A = Te^DT{-1})。但用Numpy实现时，两种计算方法得到的结果却不一致。

测试代码 #

import numpy as np
import numpy.linalg as LA

A = np.array([[3, -1],[-1, 3]])
e_values, e_vectors = LA.eig(A)
D = np.diag(e_values)
T = e_vectors.T

r1 = np.matmul(T, np.matmul(np.exp(D), LA.inv(T)))
r2 = np.exp(np.matmul(T, np.matmul(D, LA.inv(T))))

print("r1: ", r1,"\nr2: ", r2)

运行输出 #

r1:  [[29.99360307 23.60454697]
 [23.60454697 31.99360307]] 
r2:  [[20.08553692  2.71828183]
 [ 2.71828183 20.08553692]]

原因分析 #

核心问题是Numpy的np.exp()函数对矩阵执行的是逐元素指数运算，而非线性代数中定义的矩阵指数运算：

• r1的计算符合矩阵指数的推导逻辑：通过特征分解，先对对角矩阵D的元素取指数得到(e^{D)，再通过(Te}DT^{{-1})得到正确的矩阵指数(e}A)。
• r2中的np.exp(np.matmul(T, np.matmul(D, LA.inv(T))))本质是对矩阵的每个元素单独计算指数，这和通过幂级数定义的矩阵指数完全不是同一概念，自然结果不符。

解决方法 #

方法1：使用专门的矩阵指数函数 #

可以借助scipy.linalg.expm()直接计算矩阵指数，该函数实现了符合线性代数定义的矩阵指数运算：

import numpy as np
from scipy.linalg import expm

A = np.array([[3, -1],[-1, 3]])
# 直接计算矩阵指数
r3 = expm(A)
print("r3: ", r3)

运行结果会和r1完全一致。

方法2：修正特征分解的计算逻辑 #

LA.eig()返回的e_vectors是列向量为特征向量的矩阵，因此正确的特征分解公式应为(A = e_vectors @ D @ LA.inv(e_vectors))，对应的矩阵指数计算可调整为：

import numpy as np
import numpy.linalg as LA

A = np.array([[3, -1],[-1, 3]])
e_values, e_vectors = LA.eig(A)
D = np.diag(e_values)

# 正确的特征分解计算矩阵指数
r1_correct = e_vectors @ np.exp(D) @ LA.inv(e_vectors)
print("r1_correct: ", r1_correct)

该写法避免了原代码中对特征向量矩阵转置的不必要操作，结果同样与expm(A)一致。

内容的提问来源于stack exchange，提问作者zelpha_343