09 Sylvester 方程、Schur-Parlett 与矩阵指数

1. Sylvester 方程

标准形式为

$$AX-XB=C.$$

向量化得到

$$(I\otimes A-B^T\otimes I)\mathrm{vec}(X)=\mathrm{vec}(C),$$

但系数矩阵尺寸变为 $mn\times mn$，直接求解浪费结构。唯一可解的条件是

$$\sigma (A)\cap \sigma (B)=\varnothing .$$

两组谱越接近，问题越病态。

2. Bartels-Stewart 算法

先做 Schur 分解

$$A=Q_1{T}_1{Q}_1^{\ast },B^{\ast }=Q_2{T}_2{Q}_2^{\ast }.$$

令 $Y=Q_1^{\ast }X{Q}_2$、$\tilde{C}=Q_1^{\ast }C{Q}_2$，问题化为三角 Sylvester 方程

$$T_{1}Y-Y{T}_2^{\ast }=\tilde{C}.$$

再按列或分块回代求解。核心思想与普通线性方程相同：先通过稳定正交变换把系数化到易解结构。

3. Schur-Parlett 矩阵函数

若 $A=QT{Q}^{\ast }$，则

$$f(A)=Qf(T)Q^{\ast }.$$

把 $T$ 分块为

$$T=\left(\begin{array}{cc}{T}_{11}& T_{12}\\ 0& T_{22}\end{array}\right),f(T)=\left(\begin{array}{cc}f(T_{11})& X\\ 0& f(T_{22})\end{array}\right).$$

由 $Tf(T)=f(T)T$ 得到 Sylvester 方程

$$T_{11}X-X{T}_{22}=f(T_{11})T_{12}-T_{12}f(T_{22}).$$

若两块谱太近，该方程病态；因此应先重排 Schur 形，把相近特征值聚在同一对角块中。

4. Scaling and Squaring

矩阵指数满足

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

选择 $s$ 使 $\parallel A/2^s\parallel$ 足够小，先用 Padé 逼近或截断 Taylor 计算 $e^{A/2^s}$，再平方 $s$ 次。Padé 形式通常比同阶 Taylor 在较大区域内更准确。

缩放不能无限增大：每次平方都会传播误差，且对非正规矩阵，$\parallel e^A\parallel$ 可能远大于 $e^{max\mathrm{\Re }\lambda (A)}$。

5. Lyapunov 方程

连续 Lyapunov 方程常写为

$$AX+X{A}^{\ast }=-Q.$$

它是 Sylvester 方程的特殊情形。若 $A$ 稳定（特征值实部均负）且 $Q⪰0$，则

$$X={\int }_0^{\mathrm{\infty }}{e}^{tA}Q{e}^{t{A}^{\ast }}dt⪰0.$$

这连接了矩阵指数、控制系统稳定性和正定矩阵。

6. 只计算 $f(A)b$

大规模问题通常不需要完整 $f(A)$。若只关心 $f(A)b$，形成稠密矩阵既昂贵又破坏稀疏性；第 13 章会用 Krylov 投影把它近似为一个小矩阵函数。

7. Sylvester 方程的逐列递推

若 $T_2$ 上三角，将 $Y=[y_1,\dots ,y_n]$ 按列展开。第 $j$ 列满足

$$(T_1-t_{jj}^{(2)}I)y_j={\tilde{c}}_j+\sum _{k<j}{y}_k{t}_{kj}^{(2)}.$$

右端只依赖已经求出的列，因此可依次三角求解。实 Schur 形出现 $2\times 2$ 块时，应以两列为一个块同时求解，避免转入复数运算。

定义分离度

$$\mathrm{sep}(A,B)=\underset{X\ne 0}{min}\frac{\parallel AX-XB{\parallel }_F}{\parallel X{\parallel }_F}.$$

它是 Sylvester 算子的最小奇异值。$\mathrm{sep}$ 越小，解对 $C,A,B$ 的扰动越敏感；“两组谱不相交”只保证唯一性，并不保证良态。

8. 矩阵函数的定义

若 $f$ 在 $A$ 的谱邻域解析，可用 Cauchy 积分定义

$$f(A)=\frac{1}{2\pi i}{\oint }_{\mathrm{\Gamma }}f(z)(zI-A{)}^{-1}dz.$$

若 $A=X\mathrm{\Lambda }{X}^{-1}$ 可对角化，则 $f(A)=Xf(\mathrm{\Lambda })X^{-1}$；但数值算法通常不应依赖可能病态的特征向量矩阵。Jordan 形式说明了重复特征值处需要 $f$ 的导数：一个大小为 $m$ 的 Jordan 块会用到 $f,f^{\prime },\dots ,f^{(m-1)}$。

矩阵函数满足相似不变性

$$f(S^{-1}AS)=S^{-1}f(A)S$$

以及交换关系 $Af(A)=f(A)A$。Schur-Parlett 正是利用这两条性质，在正交相似变换后递推求解。

9. Parlett 递推的标量形式

当 $T$ 上三角且对角元分离时，$F=f(T)$ 的对角元为 $f_{ii}=f(t_{ii})$，非对角元可由

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

递推得到。若 $t_{ii}\approx t_{jj}$，除以小差会放大误差，因此要把相近特征值聚成块，在块内使用更稳健的 Taylor、Padé 或专用算法。

10. Padé 逼近与缩放平方

有理逼近写成

$$r_{p,q}(A)=D_q(A{)}^{-1}{N}_p(A).$$

计算时通过线性方程组 $D_q(A)X=N_p(A)$ 得到 $X$，而不是显式求逆。缩放参数 $s$ 与 Padé 阶数共同选择：$s$ 太小会让逼近误差大，太大则增加平方次数和舍入误差。成熟算法用 $\parallel A^k{\parallel }^{1/k}$ 的估计而非只看 $\parallel A\parallel$，以减少非正规矩阵上的过度缩放。

11. 特殊矩阵函数需要专门算法

• 平方根：可用 Schur 法或稳定迭代，主平方根要求谱避开非正实轴。
• 对数：常结合逆缩放平方，把矩阵反复开方后再近似 $\log (I+E)$。
• 矩阵符号函数：可用 Newton 迭代 $X_{k+1}=\frac{1}{2}(X_k+X_k^{-1})$，并用于谱投影。
• 三角函数：可通过指数关系或专用块公式计算。

把标量公式直接替换成矩阵并不总是稳定或高效；应利用函数恒等式、谱区域和矩阵结构选择算法。

12. Fréchet 导数与条件数

矩阵函数的一阶扰动写成

$$f(A+E)=f(A)+L_f(A,E)+O(\parallel E{\parallel }^2),$$

其中 $L_f(A,\cdot )$ 是 Fréchet 导数。它可由块矩阵恒等式读取：

$$f\left(\left[\begin{array}{cc}A& E\\ 0& A\end{array}\right]\right)=\left[\begin{array}{cc}f(A)& L_f(A,E)\\ 0& f(A)\end{array}\right].$$

$\parallel L_f(A)\parallel$ 给出绝对条件数。非正规矩阵即使特征值位置温和，矩阵函数仍可能对扰动高度敏感。

13. 复杂度与验证

稠密 Schur 分解和后续递推通常都是 $O(n^3)$。验证矩阵指数可检查交换残差 $\parallel A{e}^A-e^{A}A\parallel$、缩放恒等式以及小规模高精度基准；验证 Sylvester 解则直接报告

$$\frac{\parallel AX-XB-C{\parallel }_F}{\parallel A{\parallel }_F\parallel X{\parallel }_F+\parallel B{\parallel }_F\parallel X{\parallel }_F+\parallel C{\parallel }_F}.$$

14. 自检

• [ ] 能写出 Sylvester 方程唯一可解的谱条件。
• [ ] 能描述 Bartels-Stewart 的“Schur 化—三角回代”流程。
• [ ] 能从交换关系推导 Schur-Parlett 的块方程。
• [ ] 能解释 Scaling and Squaring 为什么需要在逼近与平方误差间平衡。
• [ ] 能用 $\mathrm{sep}(A,B)$ 解释 Sylvester 方程的敏感性。
• [ ] 知道为什么相近 Schur 对角元应聚成同一块。

下一章：稀疏矩阵与定常迭代。