Schur–Parlett 矩阵函数算法

形式陈述 ​

已知一个上三角矩阵的对角元素，如何恢复整个 矩阵函数公理库主矩阵函数与 Hermite 插值Primary matrix function · Matrix function by Hermite interpolation用极小多项式规定所需的函数值与导数，经 Hermite 插值定义一般方阵的函数，并区分 primary 与 principal 两种不同的分支要求。？Schur–Parlett 方法先作酉 Schur 分解公理库Schur 分解Schur decomposition · Schur triangularization · Real Schur form用酉或正交相似变换把一般矩阵化为上三角或实准上三角形，作为浮点全谱计算的稳定结构目标。 $A=QT{Q}^{\ast }$，计算 $F=f(T)$，最后返回 $QF{Q}^{\ast }$。下面假设 $f$ 在谱的邻域解析。

因为 $f(T)$ 是 $T$ 的多项式，$TF=FT$。对 $i<j$ 比较元素，得到标量 Parlett 递推：

$$(t_{ii}-t_{jj})f_{ij}=t_{ij}(f_{ii}-f_{jj})+\sum _{i<k<j}(f_{ik}{t}_{kj}-t_{ik}{f}_{kj}),f_{ii}=f(t_{ii}).$$

按超对角距离 $j-i=1,2,\dots$ 推进，右侧所需的函数元素均已完成。这个标量公式要求 $t_{ii}\ne t_{jj}$；接近的特征值也可能使分子相消后再除以小数，造成精度损失。

用块吸收接近的特征值 ​

将 Schur 形酉重排为块上三角矩阵，使需要共同处理的特征值落入同一对角块。先计算 $F_{ii}=f(T_{ii})$，然后按块超对角距离解

$$T_{ii}{F}_{ij}-F_{ij}{T}_{jj}=F_{ii}{T}_{ij}-T_{ij}{F}_{jj}+\sum _{i<k<j}(F_{ik}{T}_{kj}-T_{ik}{F}_{kj}).$$

每一步都是 Sylvester 方程公理库Sylvester 方程与谱分离Sylvester equation刻画 AX−XB=C 对所有右端唯一可解的充要条件，并用非正规矩阵说明特征值间距不能代替 sep 敏感性。。不同对角块的谱必须不交；要获得准确结果，还需关注 $\mathrm{sep}(T_{ii},T_{jj})$，不能只检查两个最近特征值的距离。块乘法的左右顺序必须保持。

完整流程是：Schur 分解、酉重排与分块、对角块函数求值、块间 Sylvester 回代、酉变换还原。对角块较小时，可用局部展开；函数有专用算法时也可直接调用。三角 Sylvester 核可复用 Bartels–Stewart 算法公理库Bartels–Stewart 矩阵方程算法Bartels–Stewart algorithm用 Schur 变换、从左下开始的块回代及坐标恢复求解矩阵方程，明确实二阶块、缩放输出、执行不变量和立方级成本。中的回代阶段。

对角块的局部展开 ​

对一个 $b\times b$ 块，取 $\sigma =\mathrm{tr}(T_{ii})/b$、$M=T_{ii}-\sigma I$。若 $f$ 在包含 $\sigma +\sigma (M)$ 的合适圆盘内有 Taylor 展开，则

$$f(T_{ii})=\sum _{k\ge 0}\frac{f^{(k)}(\sigma )}{k!}{M}^k.$$

重复特征值时，$M$ 可能幂零，展开变成有限和。接近特征值时，对角部分小，但严格上三角部分仍可能很大，因此需要实际的截断界。不能看到一项很小或为零就停止：矩阵幂的范数可能交替变大变小，某些函数的若干导数也会周期性为零。

直觉

标量递推试图用两个几乎相同的函数值恢复导数，数值上容易丢失信息。把它们合成一块，就直接在该块内部计算导数效应；只有跨越足够分离的块时，才用差值关系恢复耦合。

图中前两个特征值被放在蓝色块内，红色小分母不再出现在块间求解中。绿色耦合列由已经完成的两个对角块决定。分块不是简单交换矩阵行列：任意置换可能破坏上三角结构，真正的 Schur 重排要同时更新酉变换。

例子与边界

一个连续走向 Jordan 块的矩阵族 ​

取

$$T_{\epsilon }=\left(\begin{array}{ccc}1& 2& 1\\ 0& 1+\epsilon & 3\\ 0& 0& 4\end{array}\right),0\le \epsilon \le 1,$$

计算指数。标量递推的第一项是

$$f_{12}=2\frac{e^{1+\epsilon }-e}{\epsilon }.$$

当 $\epsilon \to 0$，它趋于 $2e$，不存在数学奇点；但直接相减会丢失有效数字。对存入双精度的 $\epsilon ={10}^{-12}$ 例，脚本的直接差商约为 $5.43694494$，高精度参照约为 $5.436563656920809$。单个二阶块可以用 expm1 型标量公式保护；分块方法则提供可推广到更大块的路线。

将前两个指标合成一块。对指数，局部展开取

$$e^{T_{11}}=e^{\sigma }\sum _{k=0}^r\frac{M^k}{k!}+\text{余项},\Vert \text{余项}{\Vert }_1\le e^{\sigma +\Vert M{\Vert }_1}\frac{\Vert M{\Vert }_1^{r+1}}{(r+1)!},$$

这里 $\sigma$ 为实数。脚本使用这个保守但明确的界，而不靠单项停止。对于本节脚本采用的 $\epsilon ={10}^{-12}$，它形成二十三次块幂更新后停止。

在极限 $\epsilon =0$，二阶块的 $M^2=0$，手算直接得到

$$F_{11}=e\left(\begin{array}{cc}1& 2\\ 0& 1\end{array}\right),F_{22}=e^4.$$

记右上耦合列为 $x$，则

$$(\left(\begin{array}{cc}1& 2\\ 0& 1\end{array}\right)-4I)x=e\left(\begin{array}{cc}1& 2\\ 0& 1\end{array}\right)\left(\begin{array}{c}1\\ 3\end{array}\right)-e^4\left(\begin{array}{c}1\\ 3\end{array}\right).$$

从第二行得 $x_2=e^4-e$，再由第一行得 $x_1=e^4-3e$。于是

$$e^{T_0}=\left(\begin{array}{ccc}e& 2e& e^4-3e\\ 0& e& e^4-e\\ 0& 0& e^4\end{array}\right).$$

重复特征值没有造成除零；其影响已经包含在对角块的导数项里。

分块也有代价和限制 ​

块太细，跨块 sep 可能很小；块太大，局部函数求值的代价和舍入风险可能增加。以特征值距离聚类可以提供启发，却不能为非正规块自动保证良好 sep。需要把聚类、局部求值误差与 Sylvester 敏感性一起考虑。

对数等非整函数还要检查展开圆盘是否避开分支切线。仅知道各特征值处有函数值，不足以证明所选 Taylor 展开可用。

推论与应用

Schur 分解、块重排及全部块间 Sylvester 求解通常是 $O(n^3)$ 级工作。若第 $i$ 个对角块大小为 $b_i$，局部 Taylor 使用 $r_i$ 次稠密块乘法，另需 $O(\sum _i{r}_i{b}_i^3)$。$r_i$ 不一定是与维数无关的常数；整块幂零展开也可能需要随块长增长的项数，所以不能无条件把一般 Taylor 版本称为立方成本。

这个框架适合没有专用实现的解析函数，也能复用一次 Schur 分解计算多个函数。对指数、对数和平方根，专用算法往往能进一步利用函数结构；通用框架的价值在于清楚划分对角块求值与跨块耦合这两个任务。

参考资料

• Philip I. Davies and Nicholas J. Higham, “A Schur–Parlett Algorithm for Computing Matrix Functions,” SIAM Journal on Matrix Analysis and Applications 25(2), 2003, pp. 464–485, §§1.1–2：块递推、局部 Taylor、停止准则和非单调项序列。作者版本
• Nicholas J. Higham, “Functions of Matrices,” Handbook of Linear Algebra, 2nd ed., 2014, §17.8：Schur–Parlett 方法的结构与限制。作者版本