Krylov 方法计算指数作用

形式陈述

大型稀疏矩阵 $A$ 的指数通常稠密，但很多应用只需要向量 $y(t)=e^{tA}b$ 。能否只做矩阵向量乘，而不形成 $e^{tA}$ ？

以下取实矩阵 $A\in\mathbb R^{n\times n}$ 、 $n\ge1$ 、 $b\in\mathbb R^n$ 和 $t\ge0$ 。负时间可将 $A$ 换成 $-A$ ，再处理对应的非负时间。设 $b\ne0$ 、 $\beta=\|b\|_2$ 。在Krylov 子空间上运行 Arnoldi，得到正交基 $V_m=(v_1,\ldots,v_m)$ ，其中 $v_1=b/\beta$ ，以及关系

AV_m=V_mH_m+h_{m+1,m}v_{m+1}e_m^T.

本页复用 Arnoldi 的构造，关注它怎样用于矩阵指数。定义近似

y_m(t)=\beta V_m e^{tH_m}e_1.

小矩阵 $H_m$ 的指数可用缩放平方算法，再将结果映回原空间。 $b=0$ 时直接返回零向量；若 Arnoldi 出现精确的 $h_{m+1,m}=0$ ，当前空间已经不变，精确运算下这个公式就是准确答案。

缺陷不是终点误差

因为 $y=e^{tA}b$ 解线性微分方程组 $y'=Ay$ 、 $y(0)=b$ ，定义近似轨道的缺陷

d_m(t)=y_m'(t)-Ay_m(t) =-\beta h_{m+1,m}v_{m+1}e_m^Te^{tH_m}e_1.

真误差 $\epsilon_m=y-y_m$ 满足

\epsilon_m'=A\epsilon_m-d_m,\quad \epsilon_m(0)=0, \qquad \epsilon_m(t)=-\int_0^t e^{(t-s)A}d_m(s)\,ds.

所以若对所有 $\tau\ge0$ 已知 $\|e^{\tau A}\|_2\le M e^{\omega\tau}$ ，才可进一步得到

\|\epsilon_m(t)\|_2 \le\int_0^tM e^{\omega(t-s)}\|d_m(s)\|_2\,ds.

只计算终点的 $\|d_m(t)\|$ ，既没有累计此前缺陷，也没有计入传播放大，通常不能直接作为误差上界。

直觉

基向量记录初值在反复作用下能够探索的方向，小矩阵 $H_m$ 则记录这些方向之间的传播。投影指数在这个小世界里精确演化，遗漏的信息只从最后一个基方向向 $v_{m+1}$ 泄出，因此缺陷有一个特别简单的方向。

图中灰框给出完整指数作用的参照分量，蓝柱给出当前近似，红色箭头是遗漏的传播方向。空间扩大后，新节点开始得到非零分量。终点缺陷可能不随维数单调减小，误差界应来自整条缺陷轨道。

例子与边界

四节点扩散，逐维看输出

取 $t=1$ 、 $b=e_1$ 和

A=\begin{pmatrix}-2&1&0&0\\1&-2&1&0\\0&1&-2&1\\0&0&1&-2\end{pmatrix}.

它对称负定，所以 $\|e^{\tau A}\|_2\le1$ 。从 $e_1$ 开始的 Arnoldi 基依次就是标准基 $e_1,e_2,e_3,e_4$ ， $H_m$ 是 $A$ 的左上 $m\times m$ 块，未结束前的 $h_{m+1,m}=1$ 。

直接指数作用约为

y(1)=(0.21526562,\ 0.18644808,\ 0.08616086,\ 0.02616210)^T.

在第一维， $H_1=(-2)$ ，因此 $y_1=e^{-2}e_1$ 。加入第二维后，两个节点之间可以来回传递；加入第三维后，第三个分量也被表示出来：

$m$	$y_m(1)^T$	实际 $2$ -范数误差	终点缺陷范数
1	$(0.13533528,0,0,0)$	$0.22194570$	$0.13533528$
2	$(0.20883325,0.15904619,0,0)$	$0.09434187$	$0.15904619$
3	$(0.21506019,0.18517912,0.07972490,0)$	$0.02697276$	$0.07972490$
4	$y(1)^T$	$0$ （精确运算）	$0$

第一行已经说明缺陷不等于误差；前两行还说明，实际误差减少时终点缺陷反而可能增加。

把缺陷积分真正算成界

本例的 $H_m$ 非对角元非负，故 $e^{sH_m}$ 逐元素非负。对 $m<4$ ，收缩性给出

\|y(1)-y_m(1)\|_2 \le\int_0^1e_m^Te^{sH_m}e_1\,ds =e_m^TH_m^{-1}(e^{H_m}-I)e_1.

这里 $H_m$ 可逆；右式实现为线性求解。三个上界依次约为 $0.43233236$ 、 $0.15769146$ 、 $0.04385172$ ，均覆盖表中的实际误差。它们比终点缺陷更贵一点，但带有明确的数学保证。一般奇异 $H_m$ 可用整函数 $\varphi_1(z)=(e^z-1)/z$ 的连续延拓计算积分，无需强行求逆。

非正规传播和分段推进

若 $A$ 非正规，仅有特征值实部为负，并不能保证 $\|e^{tA}\|\le1$ 。例如上三角的大耦合可在衰减前产生瞬态放大，这时必须保留传播因子或使用可证明的替代界。

长时间问题可按时间段重建 Krylov 空间。若每段只控制局部误差，先前误差仍会被后续传播放大；全局保证要将各段误差按传播因子累积，不能把每段同一个容差直接当作最终容差。

推论与应用

$m$ 步 Arnoldi 需要 $m$ 次矩阵向量乘；完整正交化约需 $O(nm^2)$ 工作和 $O(nm)$ 基存储，小矩阵指数在 Padé 阶数和缩放次数有固定上界时需 $O(m^3)$ 工作，否则还要计入其实际缩放工作。记稀疏表示中实际存储槽位数为 $s$ （包括重复项与显式零）。若返回完整向量，每次矩阵乘法需 $O(n+s)$ 工作，包含输出初始化； $b=0$ 的提前返回也需 $O(n)$ 输出工作。输出始终是一个 $n$ 维向量，避免了 $n^2$ 个稠密指数元素。

它与求解线性系统的 Krylov 方法共享空间，却没有自动继承某种线性残差最小化性质。这里的正确性入口是投影演化和缺陷传播；停止准则也应围绕这两个对象设计。

参考资料

Yousef Saad, “Analysis of Some Krylov Subspace Approximations to the Matrix Exponential Operator,” SIAM Journal on Numerical Analysis 29(1), 1992, pp. 209–228, §§2.1、3.3、5：投影指数、精确终止与后验误差展开。作者技术报告版本
Nicholas J. Higham and Lijing Lin, Matrix Functions: A Short Course, 2013, §7.1：Krylov 矩阵函数作用及其计算规模。作者版本

形式陈述 ​

缺陷不是终点误差 ​

四节点扩散，逐维看输出 ​

把缺陷积分真正算成界 ​

非正规传播和分段推进 ​

反向链接

形式陈述

缺陷不是终点误差

四节点扩散，逐维看输出

把缺陷积分真正算成界

非正规传播和分段推进