这是一个即兴，近期看到一篇文章专门讨论 e^A（或者 e^{tA} ） 的19种算法，令我感到震惊。

（后面预计会更一大波）

原文链接：

方法一：泰勒级数

这是最直接的方法。根据定义式，

e^A=\sum_{n=0}^{\infty}\dfrac{A^n}{n!}

我们可以将之取到收敛为止。但这样做往往收敛极慢，更好的修正是对角化

A=TDT^{-1}

那么 e^A=T\left(\sum_{n=0}^{\infty}\dfrac{D^n}{n!}\right)T^{-1}=Te^{D}T^{-1} ， e^{D} 只需把对角元计算指数就可。

这种方法对低维矩阵非常合适，但对高维矩阵（如数值解微分方程）就需要解本征值本征向量，这带来的计算量也是不可忽视的。

方法二：Pade近似

这是Python中expm函数使用的方法，其基于 e^x 的帕德近似。

有关Pade近似的东西我默认大家至少道听途说过，或者可以参考

那么 e^A 的 (p,q) 型帕德近似即为

R_{pq}(A)=\dfrac{\sum_{j=0}^p\dfrac{(p+q-j)!p!}{(p+q)!j!(p-j)!}A^j}{\sum_{j=0}^q\dfrac{(p+q-j)!q!}{(p+q)!j!(q-j)!}(-A)^j}

该方法的好处是收敛要快一些，但当 p,q 较大的时候会出现分母条件数过大的问题。

方法三：放缩+乘方

以上两种方法的精度都会受限于范数 \|A\| 。可以预料，若范数足够小，那么这两种方法都可以得到相对满意的结果。

而对于\|A\|不够小的情况，我们可以作分解

e^{A}=(e^{\frac Am})^m,\quad m\gg1

那么就可以将 \left\|\dfrac Am\right\| 化的足够小，乘方也不会带来太多误差。

方法四：切比雪夫有理多项式逼近（CRAM）

这种方法利用的是用一个有理多项式 R_{mn}(x)=\dfrac{\sum_{i=0}^{m}p_ix^i}{\sum_{i=0}^nq_ix^i} 去在一个区间内逼近所需的函数 f(x) 。

该方法在区间 (a,b) 内取 m+n+1 个点 x_i ，使逼近多项式和实际函数在这些点相同：

R_{mn}(x_i)=f(x_i)

这些点选取为归一化区间的 m+n+1 阶切比雪夫多项式 T_{m+n+1} 的零点:

x_{i}=a+\dfrac{b-a}2\cos\left[\dfrac\pi n\left(k-\dfrac12\right)\right]

上面的方程利用正常的多元方程求解方式解出系数即可，如LU分解等。

然后我们将该方法利用于 e^{-x} ，我们就可以得到一些逼近函数。已有的结果如下（看一乐）：

那么对于具有负本征值的厄米矩阵 A ,我们可以利用该逼近多项式来估计 R_{qq}(-A)\approx e^{A} 。

该方法好处是可以用已知的多项式来直接进行简单矩阵运算，但对矩阵的性质也提出了一些要求，难以推广。

该方法的基本思想是利用常微分方程组

\dfrac{\mathrm dx}{\mathrm dt}=f(x,t)

的解 x=x_0e^{tA} ，取不同的初值，然后利用不同方法数值求解微分方程组，我们即可得到 t=1 的解 x ，从而确定矩阵指数 e^A 。根据使用数值方法的不同，可以分为以下三类：

方法五：通用ODE求解器

顾名思义，就是利用机器中自带的ODE求解器来求解方程。当时原文写出的时候通用求解器并不是很丰富，如RKF45、IMPSUB等，这种方法也不易使用。当然现在就很方便了，利用matlab等软件可以轻松达到，如ode45命令等。

方法六：单步ODE方法

这种方法就是用一些简单的单步ODE求解方法来计算。一种是四阶Taylor展开

x_{j+1}=\left(I+hA+\dfrac{h^2A^2}2+\dfrac{h^3A^3}{6}+\dfrac{h^4A^4}{24}\right)x_j

还有一种就是熟知的Runge-Kutta方法

x_{j+1}=x_j+\dfrac{k_1+2k_2+2k_3+k_4}6,\quad \begin{cases} k_1=hAx_i\\k_2=hA\left(x_i+\dfrac12k_1\right)\\k_3=hA\left(x_i+\dfrac12k_2\right)\\ k_2=hA\left(x_i+k_3\right) \end{cases}

方法七：多步ODE方法

这其实就是对单步方法的一个精度提升，关联了多组变量，如Adam-Bashforth方法等。五步以内的Adam-Bashforth多步法如下

这一部分是考虑用数值方法，通过ODE方式来算出矩阵指数。后面会考虑一些矩阵指数的解析表达式。

该方法将以Hamilton-Cayley定理为基础，导出矩阵函数 e^A 的多项式的形式，并关联特征值。虽然不一定好算，但不妨为一种封闭的表达方式。

方法八： Hamilton-Cayley方法

该方法是考虑将任意 A 的幂 A^m 表示为有限多个 A^k 的多项式组合，因为根据Hamilton-Cayley定理，考虑 A 的特征多项式，

f(\lambda)=|\lambda I_n-A|=\lambda^n-\sum_{k=0}^{n-1}c_k\lambda^k\implies f(A)=0

也就是说任意 A^m 可以用 I,A,\cdots,A^{n-1} 表示：

A^m=\sum_{j=0}^{n-1}\beta_{mj}A^j

再根据 f(A)=0 我们很易得到 \beta_{mj} 的递推表达式

\beta_{mj}=\begin{cases} \delta_{mj},&mn,j=0\\c_m\beta_{m-1,n-1}+\beta_{m-1,j-1},&m>n,j>0\end{cases}

那么就有

\displaystyle e^{tA}=\sum_{m=0}^{\infty}\frac{t^mA^m}{m!}=\sum_{m=0}^{\infty}\sum_{j=0}^{n-1}\frac{t^m}{m!}\beta_{mj}A^j

定义 \displaystyle\alpha_j(t)\equiv\sum_{m=0}^{\infty}\frac{t^m}{m!}\beta_{mj} ，那么

\displaystyle e^{tA}=\sum_{j=0}^{n-1}\alpha_j(t)A^j ，就表示成了A的有限多项式形式。

缺点是 \alpha_j(t) 不易确定—— \beta_{mj} 也是要递推的。

但我们也可以换一组矩阵基——即用 A_0,\cdots,A_{n-1} 来替代 I_n,A,\cdots,A^{n-1} 。

而这在单变量函数里相当于把多项式基 1,x,\cdots,x^{n-1} 转化成了其他的多项式基（如勒让德多项式）。

同时根据插值函数的知识，我们可以借鉴拉格朗日插值多项式和牛顿插值多项式：

L(x)=\sum_{j=0}^{n}y_jl_j(x),\quad l_j(x)=\prod_{i=0,i\neq j}^n\dfrac{x-x_i}{x_j-x_i}

N(x)=\sum_{j=0}^{n}[y_0,\cdots,y_j]n_j(x),\quad n_j(x)=\prod_{i=0}^{j-1}(x-x_i)

其中高阶差分的递归定义为

[y_k]=y_k,\quad [y_k,\cdots,y_{k+j}]=\dfrac{[y_{k+1},\cdots,y_{k+j}]-[y_k,\cdots,y_{k+j-1}]}{x_{k+j}-x_k},j>0

方法九：拉格朗日插值法

考虑到 A 为n阶，那么我们猜测用n阶插值得到的就是准确结果。我们猜想 x_j=\lambda_j ，y_j 则可以用 A=\lambda_jI_n 处的值来代替：

\displaystyle e^{tA}=\sum_{j=0}^{n-1}e^{\lambda_jt}\prod_{k=1,k\neq j}^n\frac{A-\lambda_kI_n}{\lambda_j-\lambda_k}

而严格证明也很容易——令右式为 F(t) ，然后证明 F'(t)=AF(t) 即可。

方法十：牛顿插值法

同样的借鉴：

\displaystyle e^{tA}=e^{\lambda_1t}I_n+\sum_{j=2}^n[\lambda_1,\cdots,\lambda_j]\prod_{k=1}^{j-1}(A-\lambda_kI_n)

同样的证明思路。

这两种插值都是用 A 的多项式而非 A^j 来作为基。但我们是否也有一种方式来用 A^j 作为基？

方法十一：范德蒙德法

我们考虑拉格朗日插值法用到的基

\displaystyle A_j=\prod_{k=1,k\neq j}^n\frac{A-\lambda_kI_n}{\lambda_j-\lambda_k}

其展开为 A 的多项式为

\displaystyle A_j=\prod_{k=1,k\neq j}^n\frac{(-1)^{n-j}\sigma_{n-j}(\lambda_1,\cdots,\lambda_{k-1},\lambda_{k+1},\cdots,\lambda_n)}{\lambda_j-\lambda_k}A^k

而系数正好是范德蒙德矩阵

V=\begin{pmatrix} 1& 1 & \cdots & 1\\ \lambda_1 &\lambda_2 &\cdots&\lambda_n\\ \vdots & \vdots &\ddots & \vdots \\ \lambda_1^{n-1} &\lambda_2^{n-1}&\cdots &\lambda_n^{n-1} \end{pmatrix}

的逆矩阵 V^{-1} 的 (j,k) 元。（这一点见线性代数教材）也就是说

A_j=\sum_{k=1}^n(V^{-1})_{jk}A^k

进而

e^{tA}=\sum_{j,k=1}^ne^{\lambda_jt}(V^{-1})_{jk}A^k

由此便表示成了 A 的封闭多项式显式。

方法十二：拉普拉斯逆变换

这一方法是注意到了 e^{\lambda t} 的拉普拉斯变换为 \dfrac{1}{s-\lambda} 。同样对于 e^{tA} 有

\mathscr L(e^{tA})=(sI_n-A)^{-1}

但这样还不容易作逆变换，需要化成 A 的多项式。这可以根据Faddeev-LeVerrier定理实现：

\displaystyle (sI_n-A)^{-1}=\sum_{k=0}^{n-1}\frac{s^{n-k-1}}{c(s)}A_k

其中 c(s)=|sI_n-A|=s^n-\sum_{k=0}^{n-1}c_ks^k ，而 A_k 满足递推关系

A_0=I,\quad A_k=A_{k-1}A-c_{n-k}I,\quad c_{n-k}=-\dfrac{\operatorname{tr}(A_{k-1}A)}k

进而根据逆变换即可得到

\displaystyle e^{tA}=\sum_{k=0}^{n-1}\mathscr L^{-1}\left[\frac{s^{n-k-1}}{c(s)}\right]A_k

方法十三：友矩阵

该方法首先是针对一类特殊的矩阵即友矩阵（companion matrix）：

C=\begin{pmatrix} 0& 1&0 &\cdots &0\\ 0& 0 & 1 & \cdots & 0\\ \vdots & \vdots &\vdots &\ddots & \vdots \\ 0 & 0 & 0 &\cdots & 1\\ c_0 & c_1 & c_2 &\cdots & c_{n-1} \end{pmatrix}

这种矩阵由于它的稀疏性，求其幂 C^k 的计算量将远少于一般的矩阵幂。从而可以用方法三简便计算。

另外，对于所有特征值几何重数为1的矩阵，其相似于友矩阵：

A=Y^{-1}CY\implies e^A=Y^{-1}e^{C}Y

这种方法也要求幂，但是化成了一类可以快速求幂的矩阵。

前一种方法把矩阵的幂化成了矩阵的有限多项式。但对于高维矩阵，算矩阵的幂也不是一件容易事。

而矩阵分解则提供了一种手段，它将矩阵的幂转换为一类易于计算幂的矩阵的幂，进而把矩阵幂这一步骤大大简化。

方法十四：对角化

这是最直接的一种方法——直接对角化：

A=VDV^{-1},\quad e^{tA}=V\operatorname{diag}\{e^{\lambda_1t},\cdots e^{\lambda_nt}\}V^{-1}

但这种方法对于不可对角化的矩阵就失效了，比如 \begin{pmatrix} 1&1\\0&1 \end{pmatrix} 。

方法十五：上三角矩阵化

这其实就是将对角化弱化了一下——Schur分解，即正交相似于上三角矩阵（假定本征值均为实数）

Q^TAQ=T ，其中 T 为上三角（注意并没有假定 A 可以对角化）

然后考虑将 T 作类似于对角化的处理（但不是正交相似）

TR=RD ，其中 R 为上三角， D 为对角

也就是说 AQR=QRD 。那么我们计算 e^{tA}x_0 的时候也就可以作变量替换

x_0=QRy_0\implies Ry_0=Q^Tx_0

给定 Q,R,x_0 后 y_0 很容易算出。那么

e^{tA}x_0=e^{tA}QRy_0=QRe^{tD}y_0

如果是计算 e^{tA} ，那么设定不同的 x_0 （如欧氏空间基）也就很容易得到 e^{tA} 。

本征值复数也类似，只是上三角矩阵 T 要变成“准上三角”——由一系列 2\times2 块构成（见高等代数教材）。

方法十六：约当标准型

其实考虑到不一定对角化后很自然就会考虑相似于约当型的一般结论：

A=PJP^{-1},\quad J=\oplus_{i=1}^kJ_k

而约当块的非对角元部分的幂零性又给予约当块指数计算以极大的方便。从而

e^{tA}=P(\oplus_{i=1}^ke^{tJ_k})P^{-1}

方法十七：Schur分解

出发点和方法十五相同，只是这时应对的是上三角矩阵的指数：

A=Qe^{tT}Q^T

但上三角矩阵的解析函数有一个良好性质：其矩阵元可以用如下递推来求出（Parlett算法）

F\equiv f(T),\quad T_{rr}F_{rs}-F_{rs}T_{ss}=\sum_{k=0}^{s-r-1}(F_{r,r+k}T_{r+k,s}-T_{r.s-k}F_{s-k,s})

（自行验证不难）

这样便很好应对 e^{tT} 。对于复本征值的情形只需将 2\times2 块看成整体即可。

方法十八：块对角化

前面几种方法都有一个共性问题： Q(P) 的计算比较麻烦且有可能具有较大条件数。因此可以取一个折中方案即块对角化。这指的是一类scheme——将 A=SBS^{-1} 中的 B 以若干个矩阵块的形式给出。

一种方案是取出每个块的本征值平均：

B_j=\gamma_jI_j+E_j,\quad \gamma_j=\overline{\lambda_{j}^{(k)}}

这样 E_j 几乎就是幂零——从而其指数收敛很快。那么

e^{tB_j}=e^{\gamma_jt}e^{tE_j}

这种方法是基于Trotter定理：

e^{A}=\lim_{m\to\infty}(e^{\frac{B}m}e^{\frac Cm})^m,\quad A=B+C

该定理也是量子计算、量子模拟的一个基础（由Lloyd提出）。一种方案是对称化，

B=\dfrac{A+A^T}{2},\quad C=\dfrac{A-A^T}{2}

那么有 \|e^A-(e^{\frac{B}m}e^{\frac Cm})^m\|\sim O(\dfrac1m) ，同时 e^{\frac{B}m},e^{\frac{C}m} 也是比较好算的。

原文献中给出了19种理论和计算上的求矩阵指数的方法，但针对不同的问题和矩阵而言又有更多特定的算法，当矩阵A满足某种形式（如友矩阵、块对角等）时各方法的优劣也大不相同。一些表达式虽然简洁漂亮但计算起来很麻烦，也有一些方法原理上极丑无比但确实算的快（类似于矩阵的双位移QR分解）。针对不同的context（高维矩阵指数计算 or 低维一般矩阵指数理论表示 etc）可供选择的方法也有很大区别，几十年来人们也在构思更有针对性的算法。