矩阵函数计算全指南：从定义到工程实现

如果你调试过基于状态空间模型的仿真系统应该对这样一个式子不陌生x(t)e^{At}x(0)。这里的核心问题正是矩阵函数值计算——把一个函数作用在矩阵上得到一个新的矩阵而不是把A的每个元素单独带进去求函数值。矩阵函数在控制理论、振动分析、马尔可夫链、谱图理论里到处都是但很多初学者卡在第一步矩阵函数的定义到底是什么为什么有时候能拆着算有时候又不能这篇文章我就把从定义、手算方法到工程实现的完整路径梳理一遍适合正在学矩阵分析的本科生也适合工作中需要手算验证数值结果的工程师。1. 矩阵函数到底是什么搞清定义才能动手算1.1 为什么工程计算绕不开矩阵函数矩阵函数不是数学家的自娱自乐。控制理论里线性时不变系统的状态转移矩阵就是e^{At}结构动力学里振动响应的解析解经常出现sin(At)、cos(At)马尔可夫链的转移概率矩阵的整数次幂本质也是一种矩阵函数谱图理论里的热核矩阵同样是矩阵指数。几乎只要出现“耦合系统随时间演化”或“对一组线性关系做非线性变换”矩阵函数就会冒出来。有个类比我一直觉得很贴切标量函数像给单个数字做变换矩阵函数则是给一整套相互耦合的变量同时做变换同时还要保持变量之间的线性关系。你可以把它理解为“换了一套坐标系之后对每个独立模式分别处理再换回来”。这个直觉后面会反复用到因为矩阵函数计算的绝大多数方法本质上都在做“找对坐标系”这件事。1.2 三种等价的定义路径矩阵函数最常见的定义路径有三条它们在理论上等价但各有各的用处。第一条是Jordan标准型定义。任意矩阵A都可以相似于Jordan标准型J写成APJP^{-1}其中P是可逆矩阵。既然有了A^kPJ^kP^{-1}那么对解析函数f自然可以定义f(A)P f(J) P^{-1}。Jordan块上的函数值由导数级数确定这个定义最严谨是很多证明的起点。第二条是多项式插值定义。如果函数f在矩阵A的谱上足够光滑我们可以找一个次数足够低的代数多项式p使p和f在所有特征值处的函数值、导数值都相等然后直接令f(A)p(A)。这条路径是手算最常用的后面的待定系数法就来自这里。第三条是无穷级数定义。对解析函数直接代入矩阵的幂级数比如e^AΣ A^k/k!。它形式简单但实际计算时收敛问题明显直接截断往往不可靠。举个能说明问题的小例子A[[0,1],[-1,0]]。观察到A^2-I于是sin(A)的级数展开中所有偶数项消失奇数项都变成(-1)^k A最后结果是sinh(1)Acos(A)则是所有奇数项消失偶数项都变成(-1)^k I最后结果是cosh(1)I。三套定义在这个例子上给出的结果完全一致理解这一点之后你就不会被“矩阵函数到底是不是唯一”这类问题卡住了。2. 计算矩阵函数的几种主流思路2.1 特征值分解法可对角化矩阵的捷径如果A可以对角化也就是存在可逆矩阵V使AVDV^{-1}其中Ddiag(λ_1,...,λ_n)那么矩阵函数有一个非常漂亮的形式f(A) V · diag(f(λ_1), ..., f(λ_n)) · V^{-1}原因不复杂A^k V D^k V^{-1}所以任何关于A的多项式都能拆成V乘以D的多项式再乘以V^{-1}。对收敛级数定义的函数同样能逐项拆开。举一个手算例子。取A[[1,2],[2,1]]。这个矩阵对称特征值分别是3和-1对应特征向量取[1,1]^T和[1,-1]^T。令V[[1,1],[1,-1]]因为V^TV2I所以V^{-1}1/2V^T1/2[[1,1],[1,-1]]。于是f(A) 1/2 · [[f(3)f(-1), f(3)-f(-1)], [f(3)-f(-1), f(3)f(-1)]]这个公式很有用。想算e^A就把f(3)e^3、f(-1)e^{-1}代进去想算sin(A)就把sin3和sin(-1)代进去。特征值分解法把所有问题都化成了标量函数在特征值上的取值非常直观。但这里有一个前提A必须可对角化。遇到重特征值或者亏损矩阵这条路就断了需要走Jordan分解或者插值法。2.2 Jordan分解与广义特征向量并不是所有矩阵都可对角化比如A[[1,1],[0,1]]。特征值只有1但它的特征子空间是一维的找不到两个线性无关的特征向量。这类矩阵只能化成Jordan标准型。Jordan标准型由若干Jordan块组成每个Jordan块形如J_k(λ)λIN其中N是上三角的移位矩阵主对角线右上方一斜线为1其余为0。N有一个关键性质N^k0也就是幂零。对单个Jordan块矩阵函数有明确的公式f(J_k(λ)) f(λ)I f(λ)N f(λ)/2! N^2 ... f^{(k-1)}(λ)/(k-1)! N^{k-1}这就是为什么矩阵函数的定义里会出现导数Jordan块内部的信息需要靠导数来补偿。用A[[1,1],[0,1]]验证。这里AINN^20。于是e^Ae^{IN}e^I·e^Ne(IN)[[e,e],[0,e]]。如果用待定系数法也会得到同一结果。这个例子里矩阵指数出现了“e和e的线性项相乘”的结构对应到微分方程解里就是x(t)会含有te^{λt}这样的项这是控制系统里非常经典的结论。实际做数值计算时我不建议直接求Jordan分解因为广义特征向量对矩阵元素的微小扰动极其敏感舍入误差会被放大到难以接受。Jordan分解是理论工具和手算工具但几乎不是数值工具。2.3 最小多项式与待定系数插值法手算矩阵函数时我最推荐的是待定系数法它的核心是“谱插值”。思路如下矩阵的最小多项式m(λ)是满足m(A)0的最低次首一多项式。如果A的最小多项式次数是m那么f(A)一定可以表示成I,A,A^2,...,A^{m-1}的线性组合。也就是说存在系数c_0,...,c_{m-1}使f(A) c_0 I c_1 A ... c_{m-1} A^{m-1}关键是怎么确定这些c_j。规则是对每个特征值λ_i设它在最小多项式中的重数为m_i则需要保证近似多项式p(λ)Σc_jλ^j满足p(λ_i)f(λ_i)p(λ_i)f(λ_i)...p^{(m_i-1)}(λ_i)f^{(m_i-1)}(λ_i)这些条件构成一个线性方程组解出c_j即可。一个细节很重要这里用的是最小多项式中的重数而不是特征多项式中的代数重数。简单说最小多项式的重数对应最大的Jordan块尺寸不是所有同特征值Jordan块尺寸之和。用更高次多项式去插值不是不可以但方程更多、更啰嗦还容易出现冗余条件所以用最小多项式最经济。拿A[[1,1],[0,1]]举例。最小多项式是(λ-1)^2所以设f(A)c_0Ic_1A需要p(1)f(1)p(1)f(1)。第一个条件给c_0c_1f(1)第二个条件给c_1f(1)。于是f(A)[f(1)-f(1)]If(1)A。如果f(z)e^ze^A(1-e)IeA[[e,e],[0,e]]和前面完全一致。2.4 无穷级数展开简单但不总是好用很多人第一反应是直接把e^AΣA^k/k!截断到前几项。这个方法不是不行但对谱半径比较大的矩阵收敛很慢算到几十项还可能不准确。工程上计算矩阵指数现在主流做法是Scaling-and-Squaring先把A除以2^s让||A/2^s||尽量小于1再用Padé逼近算矩阵的有理近似最后反复平方s次。这相当于利用e^A(e^{A/2^s})^{2^s}把大规模矩阵函数计算变成小范数矩阵函数计算。特别要提醒一个概念性错误矩阵函数不等于逐元素函数。很多人会把np.exp(A)当成矩阵指数来用这只有在A本身是对角矩阵时才成立。反例非常直观A[[0,1],[1,0]]。用逐元素exp得到[[1,e],[e,1]]但真正的矩阵指数是[[cosh1,sinh1],[sinh1,cosh1]]。两者差别很大。所以看到“对矩阵取指数/取正弦”这类需求先想清楚要的是矩阵函数还是逐元素变换再动手写代码。3. 实战四个典型矩阵函数值计算案例3.1 案例一计算e^A常微分方程的矩阵指数看一个稍复杂一点的例子A[[2,1],[0,2]]。特征值λ2是二重根但矩阵不是对角化的因为它只有一个线性无关特征向量。最小多项式是(λ-2)^2。用待定系数法设e^Ac_0Ic_1A。条件为c_02c_1e^2 c_1e^2解得c_0e^2-2e^2-e^2c_1e^2。于是e^A -e^2Ie^2A e^2(A-I) [[e^2, e^2], [0, e^2]]用Jordan块验证更直接A2INN^20e^Ae^{2I}e^Ne^2(IN)[[e^2,e^2],[0,e^2]]。这个结果里为什么会出现“元素乘以e^2”因为Jordan块的存在。放到微分方程组xAx里解x(t)e^{At}x_0就会含有t e^{2t}这样的项。工程上遇到重根、临界阻尼、共振现象矩阵指数里出现“多项式×指数”的结构是常态。3.2 案例二计算sin(A)与cos(A)现在算A[[0,θ],[-θ,0]]的矩阵正弦函数。这个矩阵在二维旋转理论里很常见它是旋转生成元。特征值是±iθ最小多项式是λ^2θ^2。设sin(A)c_0Ic_1A条件为c_0iθc_1sin(iθ)i sinhθ c_0-iθc_1sin(-iθ)-i sinhθ两个式子相加得c_00相减得c_1sinhθ/θ。因此sin(A) (sinhθ/θ) · A [[0, sinhθ], [-sinhθ, 0]]再看cos(A)。设cos(A)d_0Id_1A条件为d_0iθd_1cos(iθ)coshθ d_0-iθd_1coshθ解得d_0coshθd_10。所以cos(A)coshθ·I是一个数量矩阵的倍数。这个结果第一次见都会有点意外对这样一个反对称矩阵取cos得到的居然是单位阵的倍数。但它完全符合谱映射规律A的特征值是±iθcos在特征值上的取值都是coshθ所以对应矩阵函数的特征值全相等矩阵函数就成了单位阵倍数。3.3 案例三矩阵平方根A^{1/2}矩阵平方根的工程用途很广协方差矩阵的白化处理、几何变换的开方操作都会用到。取A[[5,4],[4,5]]求它的主平方根。A的特征值是9和1特征向量矩阵V[[1,1],[1,-1]]V^{-1}1/2V^T。取正平方根分支则A^{1/2} V diag(3,1) V^{-1} 1/2 [[1,1],[1,-1]] [[3,0],[0,1]] [[1,1],[1,-1]]算出结果为[[2,1],[1,2]]。验证一下[[2,1],[1,2]]^2[[5,4],[4,5]]正确。注意这里有个隐藏问题矩阵平方根不是唯一的因为特征值取正分支还是负分支可以组合。比如这道题如果取diag(-3,1)会得到另一对平方根。工程上默认取主平方根也就是所有特征值取主支并要求矩阵没有非正实特征值否则结果可能进入复数域。数值计算时直接调用sqrtm会默认处理这些细节但手算或者验证结果时一定要检查分支是否选对。3.4 案例四同一矩阵多个函数的谱映射对比最后做一个综合性实验。取A[[0,1],[-1,0]]它满足A^2-I。对这个矩阵几乎所有解析函数都能快速算出来。因为A的偶数次幂交替等于I或-I奇数次幂交替等于A或-A所以f(A) [(f(i)f(-i))/2]·I [(f(i)-f(-i))/(2i)]·A把常见函数代进去得到一张很直观的对照表函数 f(z)f(A) 的最终结果e^z[[cos1, sin1], [-sin1, cos1]]sin(z)[[0, sinh1], [-sinh1, 0]]cos(z)[[cosh1, 0], [0, cosh1]]z^2[[-1, 0], [0, -1]]这个表非常值得多看两眼。它说明一件事矩阵函数f(A)的特征值确实是f(λ_i)但矩阵的结构不是只看特征值就够的——特征向量、Jordan块和函数的导数信息都会影响最终的每个元素。4. 工程实现中的数值注意事项与常见坑4.1 库函数别自己重复造轮子手算方法再好大规模场景下也别手写矩阵函数。MATLAB里用expm、logm、sqrtm、funmPython里用scipy.linalg下的expm、logm、sqrtm、funm。这些函数背后是成熟的高精度算法不是简单截断级数。一段简单的参考代码import numpy as np from scipy.linalg import expm, funm, sqrtm A np.array([[2., 1.], [0., 2.]]) # 矩阵指数结果应该是 [[exp(2), exp(2)], [0, exp(2)]] print(expm(A)) # 矩阵正弦 B np.array([[0., 1.], [-1., 0.]]) print(funm(B, np.sin)) # 矩阵平方根并验证 C np.array([[5., 4.], [4., 5.]]) R sqrtm(C) print(R) print(R R) # 应该接近 C有一点要反复强调numpy里的np.exp是逐元素指数scipy.linalg.expm才是矩阵指数。二者不能混用。对于funm这个函数它接收一个可调用对象内部会把输入当成标量值去调用但它对函数的光滑性有要求如果函数定义域包含矩阵的负特征值且不能解析延拓结果可能会带复数或直接报错。4.2 数值稳定性问题为什么数值库普遍不用特征值分解法因为特征向量矩阵V可能病态。如果A的特征值比较接近V的条件数会变得很大即使特征值计算非常精确f(A)Vf(D)V^{-1}的结果也会被放大误差。Jordan分解就更不用说了广义特征向量本身对矩阵扰动极度敏感数值上几乎不可用。因此成熟的矩阵函数库走的是另外的路对矩阵指数用Scaling-and-Squaring加Padé逼近对一般矩阵函数用Schur分解再回代。Schur分解AQTQ^H是正交变换不放大条件数算完上三角矩阵上的函数值后再通过回代过程得到f(A)的完整结果。这段逻辑理解到“为什么不能用特征分解”就够了细节交给库。4.3 我踩过的几个坑第一个坑是早期把np.exp当成矩阵指数用。当时在做一个线性系统响应仿真出来的结果总是不对最后逐行排查才发现是这里出了问题。从那以后我养成了习惯凡是结果里出现指数函数先看一眼自己调的是exp还是expm。第二个坑是平方根不验证。用sqrtm算完矩阵平方根后直接拿去做下一步计算结果发现后续方程不满足回头检查才发现算出来的R乘回去误差很大。现在我的默认流程是RR和原矩阵的二范数差值小于1e-10才继续用。第三个坑是logm的复数问题。矩阵本身是实矩阵但logm的结果可能是复数因为对数函数在负实轴上有分支。工程场景一般不想要复数输出所以算之前最好先检查特征值实部或者对矩阵做偏移处理。第四个坑是手算时忽略最小多项式。很多教材习题会给一个重特征值矩阵比如[[2,1],[0,2]]有人直接套特征多项式去做三次插值方程多、算得慢还容易把系数代错。用最小多项式或直接看出Jordan块结构会快很多。4.4 快速排查清单症状可能原因处理方式结果出现NaN或Inf级数截断、矩阵谱半径过大检查特征值范围改用库函数逐元素结果看起来“合理”但就是不对把逐元素函数当成了矩阵函数区分np.exp和expm、逐元素sin和funmsqrtm之后平方不等于原矩阵分支选择错误或精度不足用RR与原矩阵对比范数funm结果精度低或报错函数在矩阵谱上有奇点判断特征值是否在解析域内手算结果与库函数不一致忽略了重根导数条件重新确认最小多项式和插值条件做完这些排查矩阵函数计算这个环节基本就不会再拖后腿了。做了这么多年矩阵函数计算我最大的体会是先理解谱再动手算。凡是f(A)的特征值不等于f(λ_i)的情况结果基本可以判定为错。对2×2、3×3矩阵我很推荐先手算一遍不是为了替代库函数而是为了建立直觉知道函数作用在矩阵上到底会产生什么样的结构和相位。希望这篇能帮你把矩阵函数这块从“知道定义”变成“能算、会验、敢用”。