矩阵指数的缩放平方算法

形式陈述

计算矩阵函数 $e^A$ 时，直接累加幂级数可能要很多项，也可能在大项相消中损失精度。缩放平方算法先把矩阵缩小，在较容易逼近的范围内计算指数，再利用

e^A=\left(e^{A/2^s}\right)^{2^s}

恢复原来的尺度。这里 $s$ 是非负整数，最后的 $2^s$ 次幂只需 $s$ 次平方。

令 $X=2^{-s}A$ 。使用对角 Padé 有理逼近

r_m(X)=p_m(X)p_m(-X)^{-1},\qquad p_m(x)=\sum_{j=0}^m b_jx^j, \quad b_j=\frac{(2m-j)!\,m!}{(2m)!\,(m-j)!\,j!}.

分子、分母都是 $X$ 的多项式，在精确运算中交换。实现时不形成逆矩阵，而解多右端系统 $p_m(-X)R=p_m(X)$ ，然后重复 $R\leftarrow R^2$ 共 $s$ 次。

下面固定 $m=13$ ，采用 Higham 2005 的双精度范数阈值 $\theta_{13}=5.371920351148152$ ，选

s=\max\left(0,\left\lceil\log_2\frac{\|A\|_1}{\theta_{13}}\right\rceil\right),

$A=0$ 单独返回单位矩阵。先形成 $X_2=X^2$ 、 $X_4=X_2^2$ 、 $X_6=X_2X_4$ ，再计算

\begin{aligned} U&=X\bigl[X_6(b_{13}X_6+b_{11}X_4+b_9X_2) +b_7X_6+b_5X_4+b_3X_2+b_1I\bigr],\\ V&=X_6(b_{12}X_6+b_{10}X_4+b_8X_2) +b_6X_6+b_4X_4+b_2X_2+b_0I. \end{aligned}

解 $(V-U)R=V+U$ 。所有 $b_j$ 同乘一个非零常数不改变结果；脚本使用等比例的整数系数。构造 $U,V$ 共需六次矩阵乘法，随后是一次多右端求解和 $s$ 次平方。这是固定阶、未平移的教学核心；自适应实现还会按规模选更低阶，并用矩阵幂的范数改进缩放量。

设精确 Padé 值满足 $e^{-X}r_m(X)=I+G$ 且 $\|G\|<1$ 。定义 $H=\operatorname{Log}(I+G)$ ；它与 $X$ 交换，因此

r_m(X)^{2^s}=e^{A+\Delta A},\qquad \Delta A=2^sH.

\frac{\|\Delta A\|}{\|A\|} \le\frac{-\log(1-\|G\|)}{\|X\|}.

阈值来自对这种 Padé 截断后向误差的控制。它还没有计入多项式计算、线性求解和反复平方的舍入误差。后向误差与前向误差也不同：输入扰动最终还要经过指数函数的条件数放大。

直觉

缩放把远处的目标拉到逼近式擅长的邻域；Padé 用一次线性求解综合很多幂级数信息；平方则把小时间步的传播重复拼接起来。三个阶段各有不同的任务，不能只看最后的数值是否“像一个指数”。

图中同一个上三角矩阵始终保留方向耦合。超对角元素不是两个对角指数的简单平均，它记录从第二个状态向第一个状态的累计传递。

例子与边界

取

A=\begin{pmatrix}-1&20\\0&-2\end{pmatrix},\qquad e^A=\begin{pmatrix}e^{-1}&20(e^{-1}-e^{-2})\\0&e^{-2}\end{pmatrix}.

精确超对角约为 $4.65088316$ 。先故意用低阶 $m=1,s=1$ 演示机制，而不把它当作双精度推荐参数。此时 $X=A/2$ ，

R=(I-X/2)^{-1}(I+X/2) =\begin{pmatrix}3/5&16/3\\0&1/3\end{pmatrix}.

平方一次得到

R^2=\begin{pmatrix}9/25&224/45\\0&1/9\end{pmatrix}.

超对角约为 $4.97777778$ ，偏差仍很明显。缩放、求解、平方的流程正确，不代表选取的阶数和缩放量已经足够。

改用十三阶核心， $\|A\|_1=22$ ，所以 $s=3$ 。脚本得到超对角 $4.650883158696587$ ，使用九次矩阵乘法和一次多右端求解；高精度参照值约为 $4.650883158696592$ 。

令

N=\begin{pmatrix}0&10^8\\0&0\end{pmatrix}.

虽然 $\|N\|_1=10^8$ ，但 $N^2=0$ ，所以 $e^N=I+N$ 。任意 $m\ge1$ 的对角 Padé 在这个矩阵上也精确等于 $I+N$ ，无需缩放。

只看范数的十三阶规则却选出 $s=25$ ，增加二十五次平方。附带脚本在所记录的 NumPy/SciPy 环境中，未校正的范数缩放核产生约 $3.73\times10^{-9}$ 的相对误差；强制 $s=0$ 的结果则接近精确。这个具体舍入数值会随实现变化，代数上的“无需缩放”结论不变。

误差可以从对角看出：若初始求解把本应为 $1$ 的对角算成 $1-\delta$ ，反复平方会把偏差近似放大为 $2^s\delta$ 。现代算法因此利用 $\|A^k\|^{1/k}$ 而非仅用 $\|A\|$ ，并在三角情形对某些条目直接用标量公式修正。

推论与应用

稠密 $n\times n$ 输入的固定十三阶核心成本为 $O((s+1)n^3)$ ，存储为 $O(n^2)$ 。若只需要 $e^{tA}b$ ，形成完整稠密指数可能远比输出本身昂贵，Krylov 指数作用把高成本转成矩阵向量乘和一个小矩阵指数。

验算时可结合已知解析结构、不同精度和 Fréchet 条件估计。仅检查 $e^Ae^{-A}\approx I$ 并不足以证明两次计算都准确，因为它们的误差可能相关。

参考资料

Nicholas J. Higham, “The Scaling and Squaring Method for the Matrix Exponential Revisited,” SIAM Journal on Matrix Analysis and Applications 26(4), 2005, pp. 1179–1193, §2、Algorithm 2.3、Table 2.3：Padé 后向误差、阈值与高效多项式组织。作者版本
Awad H. Al-Mohy and Nicholas J. Higham, “A New Scaling and Squaring Algorithm for the Matrix Exponential,” SIAM Journal on Matrix Analysis and Applications 31(3), 2009, pp. 970–989, §1：过度缩放与矩阵幂范数的作用。作者版本