企业官网建设流程全解析

如果你调试过基于状态空间模型的仿真系统，应该对这样一个式子不陌生：x(t)=e^{At}x(0)。这里的核心问题正是矩阵函数值计算——把一个函数作用在矩阵上，得到一个新的矩阵，而不是把A的每个元素单独带进去求函数值。矩阵函数在控制理论、振动分析、马尔可夫链、谱图理论里到处都是，但很多初学者卡在第一步：矩阵函数的定义到底是什么？为什么有时候能拆着算，有时候又不能？这篇文章我就把从定义、手算方法到工程实现的完整路径梳理一遍，适合正在学矩阵分析的本科生，也适合工作中需要手算验证数值结果的工程师。

1. 矩阵函数到底是什么：搞清定义才能动手算

1.1 为什么工程计算绕不开矩阵函数

矩阵函数不是数学家的自娱自乐。控制理论里，线性时不变系统的状态转移矩阵就是e^{At}；结构动力学里，振动响应的解析解经常出现sin(At)、cos(At)；马尔可夫链的转移概率矩阵的整数次幂本质也是一种矩阵函数；谱图理论里的热核矩阵同样是矩阵指数。几乎只要出现“耦合系统随时间演化”或“对一组线性关系做非线性变换”，矩阵函数就会冒出来。

有个类比我一直觉得很贴切：标量函数像给单个数字做变换，矩阵函数则是给一整套相互耦合的变量同时做变换，同时还要保持变量之间的线性关系。你可以把它理解为“换了一套坐标系之后，对每个独立模式分别处理，再换回来”。这个直觉后面会反复用到，因为矩阵函数计算的绝大多数方法，本质上都在做“找对坐标系”这件事。

1.2 三种等价的定义路径

矩阵函数最常见的定义路径有三条，它们在理论上等价，但各有各的用处。

第一条是Jordan标准型定义。任意矩阵A都可以相似于Jordan标准型J，写成A=PJP^{-1}，其中P是可逆矩阵。既然有了A^k=PJ^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)A；cos(A)则是所有奇数项消失，偶数项都变成(-1)^k I，最后结果是cosh(1)I。三套定义在这个例子上给出的结果完全一致，理解这一点之后，你就不会被“矩阵函数到底是不是唯一”这类问题卡住了。

2. 计算矩阵函数的几种主流思路

2.1 特征值分解法：可对角化矩阵的捷径

如果A可以对角化，也就是存在可逆矩阵V使A=VDV^{-1}，其中D=diag(λ_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^TV=2I，所以V^{-1}=1/2V^T=1/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(λ)=λI+N，其中N是上三角的移位矩阵，主对角线右上方一斜线为1，其余为0。N有一个关键性质：N^k=0，也就是幂零。

对单个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]]验证。这里A=I+N，N^2=0。于是e^A=e^{I+N}=e^I·e^N=e(I+N)=[[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_0I+c_1A，需要p(1)=f(1)，p'(1)=f'(1)。第一个条件给c_0+c_1=f(1)，第二个条件给c_1=f'(1)。于是f(A)=[f(1)-f'(1)]I+f'(1)A。如果f(z)=e^z，e^A=(1-e)I+eA=[[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^A=c_0I+c_1A。条件为

c_0+2c_1=e^2 c_1=e^2

解得c_0=e^2-2e^2=-e^2，c_1=e^2。于是

e^A = -e^2I+e^2A = e^2(A-I) = [[e^2, e^2], [0, e^2]]

用Jordan块验证更直接：A=2I+N，N^2=0，e^A=e^{2I}e^N=e^2(I+N)=[[e^2,e^2],[0,e^2]]。

这个结果里为什么会出现“元素乘以e^2”？因为Jordan块的存在。放到微分方程组x'=Ax里，解x(t)=e^{At}x_0就会含有t e^{2t}这样的项。工程上遇到重根、临界阻尼、共振现象，矩阵指数里出现“多项式×指数”的结构是常态。

3.2 案例二：计算sin(A)与cos(A)

现在算A=[[0,θ],[-θ,0]]的矩阵正弦函数。这个矩阵在二维旋转理论里很常见，它是旋转生成元。

特征值是±iθ，最小多项式是λ^2+θ^2。设sin(A)=c_0I+c_1A，条件为

c_0+iθc_1=sin(iθ)=i sinhθ c_0-iθc_1=sin(-iθ)=-i sinhθ

两个式子相加得c_0=0，相减得c_1=sinhθ/θ。因此

sin(A) = (sinhθ/θ) · A = [[0, sinhθ], [-sinhθ, 0]]

再看cos(A)。设cos(A)=d_0I+d_1A，条件为

d_0+iθd_1=cos(iθ)=coshθ d_0-iθd_1=coshθ

解得d_0=coshθ，d_1=0。所以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(A)的特征值确实是f(λ_i)，但矩阵的结构不是只看特征值就够的——特征向量、Jordan块和函数的导数信息都会影响最终的每个元素。

4. 工程实现中的数值注意事项与常见坑

4.1 库函数：别自己重复造轮子

手算方法再好，大规模场景下也别手写矩阵函数。MATLAB里用expm、logm、sqrtm、funm；Python里用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分解A=QTQ^H是正交变换，不放大条件数，算完上三角矩阵上的函数值后，再通过回代过程得到f(A)的完整结果。这段逻辑理解到“为什么不能用特征分解”就够了，细节交给库。

4.3 我踩过的几个坑

第一个坑是早期把np.exp当成矩阵指数用。当时在做一个线性系统响应仿真，出来的结果总是不对，最后逐行排查才发现是这里出了问题。从那以后我养成了习惯：凡是结果里出现指数函数，先看一眼自己调的是exp还是expm。

第二个坑是平方根不验证。用sqrtm算完矩阵平方根后，直接拿去做下一步计算，结果发现后续方程不满足，回头检查才发现算出来的R乘回去误差很大。现在我的默认流程是R@R和原矩阵的二范数差值小于1e-10才继续用。

第三个坑是logm的复数问题。矩阵本身是实矩阵，但logm的结果可能是复数，因为对数函数在负实轴上有分支。工程场景一般不想要复数输出，所以算之前最好先检查特征值实部，或者对矩阵做偏移处理。

第四个坑是手算时忽略最小多项式。很多教材习题会给一个重特征值矩阵，比如[[2,1],[0,2]]，有人直接套特征多项式去做三次插值，方程多、算得慢，还容易把系数代错。用最小多项式或直接看出Jordan块结构，会快很多。

4.4 快速排查清单

做完这些排查，矩阵函数计算这个环节基本就不会再拖后腿了。

做了这么多年矩阵函数计算，我最大的体会是：先理解谱，再动手算。凡是f(A)的特征值不等于f(λ_i)的情况，结果基本可以判定为错。对2×2、3×3矩阵，我很推荐先手算一遍，不是为了替代库函数，而是为了建立直觉，知道函数作用在矩阵上到底会产生什么样的结构和相位。希望这篇能帮你把矩阵函数这块从“知道定义”变成“能算、会验、敢用”。