Bartels–Stewart 矩阵方程算法

形式陈述

求解 Sylvester 方程 $AX-XB=C$ 时，若把 $X$ 的所有元素排成向量， $n\times n$ 的未知矩阵就会变成 $n^2$ 阶线性系统。Bartels–Stewart 算法保留两侧乘法的结构，把任务变成两个较小的分解和一次有序回代。

设 $A\in\mathbb C^{m\times m}$ 、 $B\in\mathbb C^{n\times n}$ 且两谱不交。算法如下：

复 Schur 形的 $T,S$ 为上三角矩阵。设 $Y$ 的第 $j$ 列为 $y_j$ ，则

(T-s_{jj}I)y_j=d_j+\sum_{k<j}y_k s_{kj}.

因此列号 $j$ 从小到大推进；每列内部使用上三角回代，行号 $i$ 从大到小推进。标量更新为

y_{ij}=\frac{d_{ij}+\sum_{k<j}y_{ik}s_{kj} -\sum_{\ell>i}t_{i\ell}y_{\ell j}} {t_{ii}-s_{jj}}.

右端来自 $-YS$ 移项，已知列的和应取加号。

开始计算第 $j$ 列时，所有更早的列都已满足方程。计算本列第 $i$ 行时，本列更低的各行也已满足方程。上三角结构保证当前方程不再依赖其他未知量，谱不交保证分母非零。更新后不变量向上扩展；整列完成后再向右扩展。

实数据可使用实 Schur 形。这时对角线上有 $1\times1$ 或 $2\times2$ 块，后者表示复共轭特征值。仍按块列从左到右、块行从下到上处理，但当前未知量是至多 $2\times2$ 的矩阵 $Y_{IJ}$ ，要解

T_{II}Y_{IJ}-Y_{IJ}S_{JJ}=\text{已知右端}.

这是至多四个未知数的小系统，不能拆成互不相关的标量除法。

直觉

Schur 变换把两侧的耦合都整理成有方向的依赖。右乘 $S$ 只把以前的列送进当前列；左乘 $T$ 只把当前列下面的元素送进上面。于是计算从左下角开始，沿每列向上，再转向下一列。

图中蓝色是当前未知量，绿色是已经完成的量。单向箭头表示公式真正读取的已知数据。酉变换还有一项数值好处：它不放大 Frobenius 范数下的误差；但原方程的 sep 很小时，解本身仍可能敏感。

例子与边界

给定

T=\begin{pmatrix}1&2&-1\\0&2&3\\0&0&4\end{pmatrix},\quad S=\begin{pmatrix}-1&1&2\\0&-2&-1\\0&0&-3\end{pmatrix},\quad D=\begin{pmatrix}-2&7&5\\3&5&10\\10&-2&-11\end{pmatrix}.

第一列解 $(T+I)y_1=(-2,3,10)^T$ 。从底部开始： $5y_{31}=10$ 得 $y_{31}=2$ ； $3y_{21}+3\cdot2=3$ 得 $y_{21}=-1$ ；最后 $2y_{11}+2(-1)-2=-2$ 得 $y_{11}=1$ 。

第二列右端要先加上 $s_{12}y_1$ ：

d_2+y_1=(8,4,0)^T.

解 $(T+2I)y_2=(8,4,0)^T$ ，依次得到 $y_{32}=0$ 、 $y_{22}=1$ 、 $y_{12}=2$ 。

第三列右端为

d_3+2y_1-y_2=(5,7,-7)^T.

解 $(T+3I)y_3=(5,7,-7)^T$ ，依次得到 $y_{33}=-1$ 、 $y_{23}=2$ 、 $y_{13}=0$ 。最终

Y=\begin{pmatrix}1&2&0\\-1&1&2\\2&0&-1\end{pmatrix},\qquad TY-YS=D.

这九个值不仅给出答案，也逐项展示了“已经完成的列怎样修改下一列右端”。若原问题还有非平凡的 $Q,U$ ，必须进一步计算 $X=QYU^*$ ；把 $Y$ 直接当作 $X$ 会混淆坐标。

取

A=\begin{pmatrix}0&-1\\1&0\end{pmatrix},\quad B=(2),\quad C=\begin{pmatrix}1\\0\end{pmatrix}.

方程为 $(A-2I)X=C$ ，即 $-2x_1-x_2=1$ 、 $x_1-2x_2=0$ ，所以 $X=(-2/5,-1/5)^T$ 。只看两个零对角元并逐项除以 $-2$ ，会错误地得到 $(-1/2,0)^T$ 。对角二阶块内部的旋转耦合不能丢掉。

即使原始数据有限，解也可能大到溢出。LAPACK 的 DTRSYL 允许返回 $0<\mathrm{scale}\le1$ ，使计算结果满足

TY-YS=\mathrm{scale}\,D.

此时应按缩放后的方程核验残差，或在确认不会溢出后再除以 scale。把缩放后的 $Y$ 与未缩放的 $D$ 直接比较，会把正常的防溢出处理误判为算法失败。例程也会报告相同或过近特征值导致的小系统近奇异情形；这种状态应传递给调用者。

推论与应用

两个 Schur 分解需 $O(m^3+n^3)$ 工作，右端变换和恢复需 $O(m^2n+mn^2)$ ，三角回代也需 $O(m^2n+mn^2)$ 。总成本为

O\bigl(m^3+n^3+mn(m+n)\bigr),

存储为 $O(m^2+n^2+mn)$ 。当 $m=n$ 时是立方级工作；对显式 $n^2$ 阶稠密向量化系统做一般消元则是六次方级工作。这里的收益来自保留矩阵方程结构。

原方程残差可以使用归一化指标

\eta=\frac{\|C-(A\widehat X-\widehat XB)\|_F} {\|A\|_2\|\widehat X\|_F+\|B\|_2\|\widehat X\|_F+\|C\|_F},

零分母的全零情形单独处理。这个指标适合检查执行，但要把它转成解误差仍需 sep。矩阵函数的块递推、平方根的方向导数以及控制中的 Lyapunov 方程，都可以复用这个求解核心。

参考资料

Richard H. Bartels and G. W. Stewart, “Solution of the Matrix Equation $AX+XB=C$ ,” Communications of the ACM 15(9), 1972, pp. 820–826：算法的原始出处。论文记录
LAPACK, DTRSYL，参数说明及非转置分支：实准三角块回代、scale 输出和近奇异标志。官方文档与源代码
Nicholas J. Higham, “What Is the Sylvester Equation?”, 2020：Schur 变换与矩阵方程可解性。作者文章