Sylvester and Lyapunov matrix equations: the Bartels-Stewart algorithm

Intuition Beginner

You know how to solve a single linear equation like $3x=12$: divide and read off $x=4$. There is a version of this where the unknown is a whole matrix $X$, not a single number. The Sylvester equation asks you to find a matrix $X$ that satisfies $AX+XB=C$, where $A$, $B$, and $C$ are given matrices. It looks exotic, but it is still a linear problem in the entries of $X$, so there is exactly one answer in the good case, and the task is to compute it without doing anything reckless.

Why care? This equation is the workhorse of control engineering. When you ask whether a designed system settles down to rest rather than blowing up, the answer comes from a close relative called the Lyapunov equation, where $B$ is a mirror image of $A$. Solving it tells you, with certainty, whether every starting state decays over time. Aircraft autopilots, robot balancers, and power-grid regulators all lean on this test.

The naive way to solve $AX+XB=C$ is to flatten $X$ into one long column of unknowns and solve a giant ordinary linear system. That works, but it is wasteful and clumsy for sizable matrices. The clever way, due to Bartels and Stewart, first reshapes $A$ and $B$ into a staircase (triangular) form, which makes the flattened system fall apart into tiny pieces you can solve one after another. Then it reshapes back. The whole trick is to make the matrix structure do the work.

Visual Beginner

The picture shows the Bartels-Stewart idea. You cannot read off the answer to $AX+XB=C$ when $A$ and $B$ are full, tangled matrices. So you first rotate $A$ and $B$ into triangular staircases, which carries $C$ along to a new right-hand side. In the staircase form the equation splits column by column: the first column of $X$ depends only on itself, the next depends on the first, and so on, so you sweep across solving small pieces.

   FULL (hard)                  TRIANGULAR (easy, solve in order)        UNDO

   A X + X B = C                A' Y + Y B' = C'                          X = Q Y Z^T
   A,B full                     A' , B' staircases
                                                                         rotate back
   x x x   x x x                a * *      b * *        col 1: 1 small solve
   x x x   x x x      ===>      . a *      . b *        col 2: uses col 1
   x x x   x x x                . . a      . . b        col 3: uses cols 1,2
   (tangled)                    (read off, one block at a time)

Two things are worth seeing. First, putting $A$ and $B$ into staircase form is the same Schur reduction used to find eigenvalues, so the diagonal entries you end up dividing by are exactly the eigenvalues of $A$ and $B$. Second, the answer exists and is unique exactly when no eigenvalue of $A$ cancels an eigenvalue of $B$ with the opposite sign, because that cancellation is the only thing that can make a divisor land on zero.

Worked example Beginner

Solve the small Sylvester equation $AX+XB=C$ with

$$A=\left(\begin{array}{cc}2& 0\\ 0& 3\end{array}\right),B=\left(\begin{array}{cc}1& 0\\ 0& 4\end{array}\right),C=\left(\begin{array}{cc}6& 12\\ 8& 35\end{array}\right).$$

Both $A$ and $B$ are already diagonal, so this is the staircase case and we can read the answer off entry by entry. Write $X$ with unknown entries $x_{ij}$.

Step 1. Find the equation for one entry. The $(i,j)$ entry of $AX$ is $a_{ii}{x}_{ij}$ because $A$ is diagonal, and the $(i,j)$ entry of $XB$ is $x_{ij}{b}_{jj}$ because $B$ is diagonal. So each entry obeys $a_{ii}{x}_{ij}+x_{ij}{b}_{jj}=c_{ij}$, that is $(a_{ii}+b_{jj})x_{ij}=c_{ij}$.

Step 2. Solve each entry. Divide by the sum $a_{ii}+b_{jj}$:

$$x_{11}=\frac{6}{2+1}=2,x_{12}=\frac{12}{2+4}=2,x_{21}=\frac{8}{3+1}=2,x_{22}=\frac{35}{3+4}=5.$$

Step 3. Assemble the answer: $X=\left(\begin{array}{cc}2& 2\\ 2& 5\end{array}\right)$. You can check directly that $AX+XB$ rebuilds $C$.

What this tells us: in the diagonal case every entry of $X$ is just the matching entry of $C$ divided by a sum of one eigenvalue from $A$ and one from $B$. The solution exists because none of those four sums $a_{ii}+b_{jj}$ is zero. If one were zero, that entry could not be solved, and the equation would have no unique answer.

Check your understanding Beginner

Formal definition Intermediate+

Let $A\in M_m(\mathbb{C})$, $B\in M_n(\mathbb{C})$, and $C\in M_{m\times n}(\mathbb{C})$. The Sylvester equation is the linear matrix equation

$$AX+XB=C,X\in M_{m\times n}(\mathbb{C}),$$

for the unknown $X$. The left-hand side defines the Sylvester operator $\mathcal{S}:M_{m\times n}(\mathbb{C})\to M_{m\times n}(\mathbb{C})$, $\mathcal{S}(X)=AX+XB$, a linear map on the matrix space; the equation reads $\mathcal{S}(X)=C$.

Definition (vec and Kronecker form). The vec operator stacks the columns of a matrix into one tall vector: $\mathrm{vec}(X)\in {\mathbb{C}}^{mn}$ lists the first column of $X$, then the second, and so on. The Kronecker product $P\otimes Q$ is the block matrix whose $(i,j)$ block is $p_{ij}Q$. With the standard identity $\mathrm{vec}(PXQ)=(Q^{\mathsf{T}}\otimes P)\mathrm{vec}(X)$, applied to $AX=AX{I}_n$ and $XB=I_{m}XB$, the Sylvester equation becomes the ordinary linear system

$$(I_n\otimes A+B^{\mathsf{T}}\otimes I_m)\mathrm{vec}(X)=\mathrm{vec}(C),$$

whose coefficient matrix is the Kronecker matrix of the Sylvester operator. Here $I_n\otimes A$ realises $X\mapsto AX$ and $B^{\mathsf{T}}\otimes I_m$ realises $X\mapsto XB$ after column stacking.

Definition (the Lyapunov equation). Taking $m=n$, $B=A^{\mathsf{T}}$, and $C=-Q$ with $Q=Q^{\mathsf{T}}$ gives the Lyapunov equation

$$A^{\mathsf{T}}X+XA=-Q,$$

the symmetric specialisation of the Sylvester equation; when $Q$ is symmetric the solution $X$ is symmetric. A matrix $A$ is Hurwitz (stable) when every eigenvalue has strictly negative real part, the condition under which the linear system $\dot{x}=Ax$ of 02.06.03 decays to the origin.

Definition (separation). The separation of $A$ and $-B$ is

$$\mathrm{sep}(A,-B)=\underset{\Vert X{\Vert }_F=1}{\min }\Vert AX+XB{\Vert }_F={\sigma }_{\min }(I_n\otimes A+B^{\mathsf{T}}\otimes I_m),$$

the smallest singular value of the Sylvester operator's Kronecker matrix. It is positive exactly when $\mathcal{S}$ is invertible, and its reciprocal controls how rounding errors and perturbations in $A$, $B$, $C$ propagate into $X$.

Notation: $\mathcal{S}$ is the Sylvester operator; $\mathrm{vec}(X)$ is the column-stacking map; $P\otimes Q$ is the Kronecker product; $\sigma (A)$ is the spectrum; $\mathrm{sep}(A,-B)$ is the separation; ${\sigma }_{\min }$ is the smallest singular value; $\Vert \cdot {\Vert }_F$ is the Frobenius norm. The spectrum of $\mathcal{S}$, computed in the Key theorem, is $\{\lambda +\mu :\lambda \in \sigma (A),\mu \in \sigma (B)\}$.

Counterexamples to common slips

• Solvability is not about $A$ and $B$ sharing eigenvalues; it is about $\sigma (A)$ meeting $\sigma (-B)$. For A = \begin{psmallmatrix} 1 & 0 \\ 0 & 2 \end{psmallmatrix} and B = \begin{psmallmatrix} -1 & 0 \\ 0 & 5 \end{psmallmatrix}, the matrices share no eigenvalue, yet $\lambda =1$ and $\mu =-1$ give $\lambda +\mu =0$, so the Sylvester operator is singular and the equation may have no solution or infinitely many.

• The Sylvester operator's spectrum is the set of sums $\lambda +\mu$, not products or ratios. A reader who imports the generalised-eigenvalue ratio $\lambda /\mu$ of 43.06.10 here will mis-locate the singular cases; the additive structure comes from $AX+XB$ being additive in its two terms.

• A symmetric right-hand side $Q$ does not by itself force a positive-definite Lyapunov solution. With A = \begin{psmallmatrix} 1 & 0 \\ 0 & 1 \end{psmallmatrix} (not Hurwitz) and $Q=I$, the equation $A^{\mathsf{T}}X+XA=-Q$ gives $2X=-I$, so $X=-\frac{1}{2}I$ is negative-definite; definiteness of $X$ tracks the stability of $A$, the content of the Lyapunov theorem of 02.12.08, not merely the sign of $Q$.

Key theorem with proof Intermediate+

Theorem (solvability via the Bartels-Stewart reduction; Bartels-Stewart 1972 ^{[source pending]}; Golub-Van Loan §7.6 ^{[source pending]}). The Sylvester equation $AX+XB=C$ has a unique solution for every $C$ if and only if $\sigma (A)\cap \sigma (-B)=\varnothing$. When this holds, the solution is computed by the Bartels-Stewart algorithm: reduce $A$ to real Schur form $A=U{R}_A{U}^{\mathsf{T}}$ and $B$ to real Schur form $B=V{R}_B{V}^{\mathsf{T}}$ by the QR algorithm of 43.06.03, with $R_A$ lower quasi-triangular and $R_B$ upper quasi-triangular; transform the right-hand side to $\tilde{C}=U^{\mathsf{T}}CV$; solve the quasi-triangular equation $R_{A}Y+Y{R}_B=\tilde{C}$ for $Y$ by back-substitution over $1\times 1$ and $2\times 2$ diagonal blocks; and recover $X=UY{V}^{\mathsf{T}}$. The method costs $\mathcal{O}(m^3+n^3)$ and is backward stable.

Proof. The eigenvalues of the Sylvester operator $\mathcal{S}$ are the sums $\lambda +\mu$ over $\lambda \in \sigma (A)$, $\mu \in \sigma (B)$; this spectral fact, proven in the Full proof set, makes $\mathcal{S}$ invertible exactly when no sum vanishes, that is when no $\lambda$ equals $-\mu$, that is when $\sigma (A)\cap \sigma (-B)=\varnothing$. Invertibility of $\mathcal{S}$ is precisely the existence of a unique solution for every $C$, establishing the equivalence.

For the algorithm, substitute the Schur factorisations. Replacing $A$, $B$ and left-right multiplying by the orthogonal factors:

$$AX+XB=C⟺U{R}_A{U}^{\mathsf{T}}X+XV{R}_B{V}^{\mathsf{T}}=C.$$

Multiply on the left by $U^{\mathsf{T}}$ and on the right by $V$, and insert $I=V{V}^{\mathsf{T}}$ and $I=U^{\mathsf{T}}U$ to expose $Y=U^{\mathsf{T}}XV$:

$$R_A(U^{\mathsf{T}}XV)+(U^{\mathsf{T}}XV)R_B=U^{\mathsf{T}}CV,\text{i.e.}{R}_{A}Y+Y{R}_B=\tilde{C}.$$

This transformed equation has the same operator spectrum because $R_A$, $R_B$ are orthogonally similar to $A$, $B$ and so share their eigenvalues; the reduction is therefore lossless and preserves solvability. Now read the transformed equation column by column. Let $y_j$ and ${\tilde{c}}_j$ be the $j$-th columns of $Y$ and $\tilde{C}$, and let $R_B=(r_{kj})$ be upper quasi-triangular. The $j$-th column of $Y{R}_B$ is ${\sum }_{k\le j}{r}_{kj}{y}_k$, so the $j$-th column of the equation is

$$R_A{y}_j+\sum _{k\le j}{r}_{kj}{y}_k={\tilde{c}}_j,\text{hence}(R_A+r_{jj}I)y_j={\tilde{c}}_j-\sum _{k<j}{r}_{kj}{y}_k.$$

The right-hand side uses only earlier columns $y_1,\dots ,y_{j-1}$, so the columns are found in order $j=1,2,\dots ,n$. Each step solves a quasi-triangular system with coefficient matrix $R_A+r_{jj}I$, itself solved by back-substitution up the rows of $R_A$, the $2\times 2$ diagonal blocks of $R_A$ (and of $R_B$, when a $2\times 2$ block couples two columns) handled as small dense solves. Each such small solve is non-singular because its eigenvalues are ${\lambda }_i+{\mu }_j\ne 0$ under the spectral hypothesis. Reversing the transformation, $X=UY{V}^{\mathsf{T}}$, returns the solution to the original equation. The Schur reductions cost $\mathcal{O}(m^3)$ and $\mathcal{O}(n^3)$, the transformations of $C$ and back are $\mathcal{O}(mn(m+n))$, and the back-substitution is $\mathcal{O}(m^{2}n+m{n}^2)$; backward stability follows because every transformation is orthogonal and the small block solves are backward stable. $\square$

Bridge. This solvability-plus-algorithm result builds toward the Lyapunov and control-theory applications of the Advanced results, where the same Schur back-substitution computes the stability Gramians of a linear system. The Bartels-Stewart reduction is exactly the same orthogonal Schur machinery that the QR algorithm of 43.06.03 produces, repurposed so that a matrix-equation solve inherits its backward stability; this is the foundational reason the method never forms the giant $mn\times mn$ Kronecker system in full but instead lets the triangular structure split it into $n$ small column solves. The operator-spectrum criterion $\sigma (A)\cap \sigma (-B)=\varnothing$ appears again in the matrix-function unit 43.06.11, where the block Parlett recurrence is a Sylvester solve across separated eigenvalue clusters, and this is exactly the well-posedness proposition proven there, which the present unit cross-references rather than reproving. The central insight is that an equation linear in a matrix unknown is an ordinary linear system in disguise, and putting these together, the bridge is the recognition that triangularising the coefficient matrices triangularises that disguised system, turning an $\mathcal{O}(m^3{n}^3)$ dense solve into an $\mathcal{O}(m^3+n^3)$ structured one.

Exercises Intermediate+

Advanced results Master

Theorem (conditioning via separation; Golub-Van Loan §7.6 ^{[source pending]}). Let $\mathcal{S}(X)=AX+XB$ be invertible, with $\mathrm{sep}(A,-B)={\min }_{\Vert X{\Vert }_F=1}\Vert AX+XB{\Vert }_F>0$. If $\hat{X}$ solves the perturbed equation $(A+\Delta A)\hat{X}+\hat{X}(B+\Delta B)=C+\Delta C$ with the perturbations small relative to $\mathrm{sep}(A,-B)$, then

$$\frac{\Vert \hat{X}-X{\Vert }_F}{\Vert X{\Vert }_F}\le \frac{\Vert \Delta C{\Vert }_F+(\Vert \Delta A{\Vert }_F+\Vert \Delta B{\Vert }_F)\Vert X{\Vert }_F}{\mathrm{sep}(A,-B)\Vert X{\Vert }_F}+O({\delta }^2).$$

Thus $1/\mathrm{sep}(A,-B)$ is the condition number of the Sylvester solve, and the separation can be far smaller than the smallest gap ${\min }_{i,j}|{\lambda }_i+{\mu }_j|$ between operator eigenvalues when $A$ or $B$ is non-normal.

The separation, not the bare eigenvalue gap, is the right conditioning measure because $\mathrm{sep}(A,-B)$ is the smallest singular value of the Kronecker matrix $I_n\otimes A+B^{\mathsf{T}}\otimes I_m$, whereas the eigenvalue gap is its smallest eigenvalue modulus; for a non-normal operator these diverge, and a Sylvester equation with well-separated operator eigenvalues can still be severely ill-conditioned. Estimating $\mathrm{sep}(A,-B)$ without forming the $mn\times mn$ Kronecker matrix is done by a power-iteration on ${\mathcal{S}}^{-1}{\mathcal{S}}^{-\ast }$, each application a Bartels-Stewart solve on the Schur factors already in hand, so the condition estimate reuses the same triangular machinery at the cost of a few extra solves.

Theorem (Hurwitz stability and the Lyapunov equation; Golub-Van Loan §7.6.3 ^{[source pending]}; Lancaster-Tismenetsky Ch. 12 ^{[source pending]}). For $A\in M_n(\mathbb{R})$ the following are equivalent: (i) $A$ is Hurwitz; (ii) for some symmetric positive-definite $Q$ the Lyapunov equation $A^{\mathsf{T}}X+XA=-Q$ has a symmetric positive-definite solution $X$; (iii) for every symmetric positive-definite $Q$ that equation has a unique symmetric positive-definite solution, given by the convergent Gramian $X={\int }_0^{\infty }{e}^{A^{\mathsf{T}}t}Q{e}^{At}dt$. The solution $X$ is the observability Gramian when $Q=C^{\mathsf{T}}C$, and the controllability Gramian solves the dual equation $AP+P{A}^{\mathsf{T}}=-B{B}^{\mathsf{T}}$.

This is the algebraic face of the Lyapunov direct method of 02.12.08: a positive-definite $X$ solving the Lyapunov equation is exactly a quadratic Lyapunov function $V(x)=x^{\mathsf{T}}Xx$ whose derivative along trajectories of $\dot{x}=Ax$ is $\dot{V}=-x^{\mathsf{T}}Qx<0$, certifying decay to the origin. The equivalence converts the analytic stability question — do all trajectories decay? — into a single matrix-equation solve, and the Bartels-Stewart algorithm makes that solve a backward-stable $\mathcal{O}(n^3)$ computation; the eigenvalue-free certificate this provides is the basis of how control software verifies stability without ever computing the spectrum of $A$. The controllability and observability Gramians it produces feed balanced truncation, the standard model-reduction method, where the Hankel singular values are the square roots of the eigenvalues of the product of the two Gramians.

Theorem (Hammarling's Cholesky-factor variant; Hammarling 1982 ^{[source pending]}). Let $A$ be Hurwitz and let the right-hand side be $Q=G^{\mathsf{T}}G$ in factored form. Then the unique positive-definite solution of $A^{\mathsf{T}}X+XA=-G^{\mathsf{T}}G$ admits a Cholesky factorisation $X=U^{\mathsf{T}}U$ with $U$ upper-triangular, and $U$ can be computed directly from the real Schur form of $A$ and the factor $G$, column block by column block, without ever forming $X$ or $Q$.

Computing $U$ rather than $X$ halves the storage of the symmetric solution and, more importantly, guarantees the computed solution is positive-semidefinite to working precision — a property that forming $X=U^{\mathsf{T}}U$ explicitly and then refactoring can lose to rounding. The Hammarling recursion threads the Cholesky factorisation through the same Schur back-substitution as Bartels-Stewart: at each diagonal block it solves a small triangular Lyapunov equation for a block of $U$, and the off-diagonal blocks fall out of a Sylvester solve against the already-computed blocks. This factored form is the numerically preferred route in model-reduction codes, where the Gramian factors, not the Gramians, are what balanced truncation consumes.

Synthesis. The Sylvester equation is a linear system wearing matrix clothing, and the equivalence of its three faces — the matrix equation $AX+XB=C$, the operator equation $\mathcal{S}(X)=C$, and the Kronecker system $(I_n\otimes A+B^{\mathsf{T}}\otimes I_m)\mathrm{vec}(X)=\mathrm{vec}(C)$ — is the central insight that organises every algorithm and every conditioning bound for it. This is exactly why the operator spectrum is the additive set $\{\lambda +\mu \}$ and why solvability is the separation condition $\sigma (A)\cap \sigma (-B)=\varnothing$, the well-posedness proposition that 43.06.11 proves and the present unit cross-references; the Bartels-Stewart algorithm generalises the scalar divide-and-read-off recipe by triangularising $A$ and $B$ through the Schur reduction of 43.06.03, so the disguised Kronecker matrix triangularises with them and the $\mathcal{O}(m^3{n}^3)$ dense solve collapses into $n$ structured column solves of cost $\mathcal{O}(m^3+n^3)$ total. The Lyapunov specialisation is dual to the Sylvester general case: setting $B=A^{\mathsf{T}}$ and $Q=-C$ turns the solve into a stability certificate, and the Gramian integral ${\int }_0^{\infty }{e}^{A^{\mathsf{T}}t}Q{e}^{At}dt$ ties the algebraic equation back to the matrix exponential of 02.06.03 and the analytic Lyapunov direct method of 02.12.08, so the computed solution $X$ is simultaneously a number a computer returns and a quadratic Lyapunov function a theorem guarantees.

Putting these together, the bridge from the static spectral theory of $A$ and $B$ to a production solver is the recognition that the separation $\mathrm{sep}(A,-B)$ — the smallest singular value of the operator, not the smallest gap of its spectrum — is what conditions the answer, the foundational reason a problem with well-separated operator eigenvalues can still be ill-conditioned. Hammarling's factored variant preserves the definiteness the stability application demands, and the entire computation rides on the orthogonal Schur form, never on the eigenvector basis that non-normality would corrupt.

Full proof set Master

Proposition (spectrum of the Sylvester operator). The eigenvalues of $\mathcal{S}:X\mapsto AX+XB$ on $M_{m\times n}(\mathbb{C})$ are exactly the sums $\lambda +\mu$ with $\lambda \in \sigma (A)$ and $\mu \in \sigma (B)$, counted with the products of the corresponding generalised-eigenspace dimensions. Hence $\mathcal{S}$ is invertible iff $\sigma (A)\cap \sigma (-B)=\varnothing$.

Proof. Pass to the Kronecker matrix $K=I_n\otimes A+B^{\mathsf{T}}\otimes I_m$, which represents $\mathcal{S}$ in the basis ordered by $\mathrm{vec}$, so $\sigma (\mathcal{S})=\sigma (K)$. Triangularise both coefficient matrices by Schur: $A=U{R}_A{U}^{\mathsf{H}}$ with $R_A$ upper-triangular carrying $\sigma (A)$ on its diagonal, and $B^{\mathsf{T}}=W{R}_B{W}^{\mathsf{H}}$ with $R_B$ upper-triangular carrying $\sigma (B^{\mathsf{T}})=\sigma (B)$ on its diagonal. Conjugate $K$ by the unitary $W\otimes U$ and apply the mixed-product law $(P\otimes Q)(R\otimes S)=(PR)\otimes (QS)$: $$ (W \otimes U)^{\mathsf H} K (W \otimes U) = I_n \otimes R_A + R_B \otimes I_m, $$ since $(W\otimes U{)}^{\mathsf{H}}(I_n\otimes A)(W\otimes U)=(W^{\mathsf{H}}W)\otimes (U^{\mathsf{H}}AU)=I_n\otimes R_A$ and similarly the second term gives $R_B\otimes I_m$. The matrix $I_n\otimes R_A+R_B\otimes I_m$ is upper-triangular: $I_n\otimes R_A$ is block-diagonal with the upper-triangular block $R_A$ repeated, and $R_B\otimes I_m$ is upper-block-triangular with diagonal blocks the scalar multiples $({\mu }_j)I_m$. Its diagonal entries are therefore the sums of the $i$-th diagonal of $R_A$ and the $j$-th diagonal of $R_B$, namely ${\lambda }_i+{\mu }_j$ ranging over all $i,j$. These are the eigenvalues of $K$, hence of $\mathcal{S}$. Invertibility fails exactly when some ${\lambda }_i+{\mu }_j=0$, that is ${\lambda }_i=-{\mu }_j$, that is $\sigma (A)\cap \sigma (-B)\ne \varnothing$. $\square$

Proposition (Bartels-Stewart column recurrence is well-posed). Under $\sigma (A)\cap \sigma (-B)=\varnothing$, the quasi-triangular equation $R_{A}Y+Y{R}_B=\tilde{C}$ (with $R_A$, $R_B$ in real Schur form) is solved column-block by column-block, and each block solve $(R_A+r_{jj}I)y_j={\tilde{c}}_j-{\sum }_{k<j}{r}_{kj}{y}_k$ has a non-singular coefficient matrix.

Proof. Order the columns $j=1,\dots ,n$. Expanding the $j$-th column of $Y{R}_B$ with $R_B$ upper quasi-triangular gives ${\sum }_{k\le j}{r}_{kj}{y}_k$, so the $j$-th column of $R_{A}Y+Y{R}_B=\tilde{C}$ reads $R_A{y}_j+r_{jj}{y}_j+{\sum }_{k<j}{r}_{kj}{y}_k={\tilde{c}}_j$, that is $(R_A+r_{jj}I)y_j={\tilde{c}}_j-{\sum }_{k<j}{r}_{kj}{y}_k$, whose right-hand side involves only previously computed columns. The coefficient matrix $R_A+r_{jj}I$ has eigenvalues ${\lambda }_i+r_{jj}$ where ${\lambda }_i\in \sigma (A)$ and $r_{jj}$ is a diagonal entry of $R_B$, hence an eigenvalue ${\mu }_j\in \sigma (B)$. Under the hypothesis no ${\lambda }_i+{\mu }_j$ vanishes, so $R_A+r_{jj}I$ is non-singular and $y_j$ is determined. When a $2\times 2$ block of $R_B$ couples columns $j,j+1$, the pair is solved jointly as a $2m\times 2m$ system whose coefficient matrix is $I_2\otimes R_A+R_B^{(jj)}\otimes I_m$ with $R_B^{(jj)}$ the $2\times 2$ block; its eigenvalues are ${\lambda }_i+\mu$ for the two (complex-conjugate) eigenvalues $\mu$ of the block, again non-zero, so the joint solve is non-singular. $\square$

Proposition (uniqueness and symmetry of the Lyapunov solution). For Hurwitz $A$ and symmetric $Q$, the Lyapunov equation $A^{\mathsf{T}}X+XA=-Q$ has a unique solution, and that solution is symmetric.

Proof. Uniqueness is the Sylvester criterion with $B=A$: the operator $X\mapsto A^{\mathsf{T}}X+XA$ has spectrum $\{\bar{\nu }+{\nu }^{\prime }:\nu ,{\nu }^{\prime }\in \sigma (A)\}$ where the first summand ranges over $\sigma (A^{\mathsf{T}})=\sigma (A)$. (Writing the eigenvalues of $A$ as $\nu =\xi +i\eta$, a sum $\nu +{\nu }^{\prime }=0$ would force $\mathrm{Re}\nu +\mathrm{Re}{\nu }^{\prime }=0$; but $A$ Hurwitz makes every real part strictly negative, so every such sum has strictly negative real part and cannot vanish.) Thus the operator is invertible and the solution is unique. For symmetry, transpose the equation: $(A^{\mathsf{T}}X+XA{)}^{\mathsf{T}}=A^{\mathsf{T}}{X}^{\mathsf{T}}+X^{\mathsf{T}}A$, and $(-Q{)}^{\mathsf{T}}=-Q$, so $X^{\mathsf{T}}$ solves the same Lyapunov equation; uniqueness forces $X^{\mathsf{T}}=X$. $\square$

Proposition (Gramian formula and positive-definiteness). For Hurwitz $A$ and symmetric positive-definite $Q$, $X={\int }_0^{\infty }{e}^{A^{\mathsf{T}}t}Q{e}^{At}dt$ converges, solves $A^{\mathsf{T}}X+XA=-Q$, and is symmetric positive-definite.

Proof. Convergence: since $A$ is Hurwitz, $\Vert e^{At}\Vert \le c{e}^{-\alpha t}$ for some $c,\alpha >0$, so the integrand is bounded in norm by $c^2\Vert Q\Vert e^{-2\alpha t}$ and the integral converges absolutely. That $X$ solves the equation follows by integrating the derivative identity $\frac{d}{dt}(e^{A^{\mathsf{T}}t}Q{e}^{At})=A^{\mathsf{T}}(e^{A^{\mathsf{T}}t}Q{e}^{At})+(e^{A^{\mathsf{T}}t}Q{e}^{At})A$ from $0$ to $\infty$: the left side telescopes to $[e^{A^{\mathsf{T}}t}Q{e}^{At}{]}_0^{\infty }=-Q$, the right to $A^{\mathsf{T}}X+XA$. Symmetry: each integrand satisfies $(e^{A^{\mathsf{T}}t}Q{e}^{At}{)}^{\mathsf{T}}=e^{A^{\mathsf{T}}t}Q{e}^{At}$ because $Q=Q^{\mathsf{T}}$, so $X=X^{\mathsf{T}}$. Positive-definiteness: for $v\ne 0$, $v^{\mathsf{T}}Xv={\int }_0^{\infty }(e^{At}v{)}^{\mathsf{T}}Q(e^{At}v)dt$; the integrand is $\ge {\lambda }_{\min }(Q)\Vert e^{At}v{\Vert }^2\ge 0$ and strictly positive near $t=0$ since $e^{A\cdot 0}v=v\ne 0$ and $e^{At}$ is invertible for all $t$, so the integral over the positive-measure neighbourhood of $0$ is strictly positive. Hence $v^{\mathsf{T}}Xv>0$. $\square$

Connections Master

The Bartels-Stewart algorithm is built directly on the real Schur form produced by the QR algorithm of 43.06.03: the entire method is a Schur reduction of $A$ and $B$ followed by quasi-triangular back-substitution, so the backward stability and the $\mathcal{O}(n^3)$ cost are inherited from the QR reduction, and the $2\times 2$ diagonal blocks that the real Schur form carries for complex-conjugate eigenvalue pairs are exactly the blocks the column recurrence handles as small dense solves rather than scalar divisions.

The well-posedness of the Sylvester equation is the operator-spectrum proposition $\sigma (\mathcal{S})=\{\lambda +\mu \}$ that already appears in the matrix-function unit 43.06.11, where the block Parlett recurrence $T_{ii}{F}_{ij}-F_{ij}{T}_{jj}=C_{ij}$ is a Sylvester solve across well-separated eigenvalue clusters; that unit owns and proves the solvability proposition for the block-triangular subroutine, and the present unit specialises it to the general dense $AX+XB=C$ and supplies the algorithm for two arbitrary matrices, so the two units are the well-posedness theorem and its general solver.

The Lyapunov equation is the algebraic counterpart of the analytic Lyapunov direct method of 02.12.08: a positive-definite solution $X$ of $A^{\mathsf{T}}X+XA=-Q$ is precisely a quadratic Lyapunov function $V(x)=x^{\mathsf{T}}Xx$ certifying asymptotic stability of $\dot{x}=Ax$, and the Gramian formula $X={\int }_0^{\infty }{e}^{A^{\mathsf{T}}t}Q{e}^{At}dt$ expresses that solution through the matrix exponential of 02.06.03, tying the matrix-equation solve to the flow whose decay it certifies.

The eigenvalue-and-spectrum theory of 01.01.08 underlies the solvability criterion: the Sylvester operator's eigenvalues are the pairwise sums of the eigenvalues of $A$ and $B$, so the characteristic-polynomial and spectrum machinery of that foundational unit is what the separation condition $\sigma (A)\cap \sigma (-B)=\varnothing$ rests on, and the degenerate case where a sum vanishes is the matrix-equation analogue of a vanishing characteristic value.

The separation $\mathrm{sep}(A,-B)$ and the generalised-eigenvalue conditioning of 43.06.10 are siblings in invariant-subspace sensitivity: the same separation that conditions the Sylvester solve governs the sensitivity of the invariant subspaces spanned by clustered eigenvalues, so a small $\mathrm{sep}$ signals both an ill-conditioned matrix-equation solve and a poorly determined invariant subspace, linking matrix equations to the perturbation theory of the eigenproblem.

Historical & philosophical context Master

The equation $AX+XB=C$ carries the name of James Joseph Sylvester, who studied the operator $X\mapsto AX+XB$ and its spectrum in the nineteenth-century theory of bilinear forms and showed that its eigenvalues are the pairwise sums of those of $A$ and $B$, the result that fixes when the equation is uniquely solvable ^{[Sylvester 1884]}. The symmetric specialisation $A^{\mathsf{T}}X+XA=-Q$ is named for Aleksandr Lyapunov, whose 1892 dissertation on the stability of motion introduced the direct method in which a positive-definite quadratic form decreasing along trajectories certifies stability, the analytic source of the matrix-equation stability criterion ^{[Lyapunov 1892]}.

The numerical breakthrough was the 1972 algorithm of Richard Bartels and G. W. Stewart, which recognised that reducing $A$ and $B$ to Schur form by the then-recent QR algorithm turns the dense Sylvester equation into a quasi-triangular system solvable by back-substitution, giving the first backward-stable $\mathcal{O}(n^3)$ dense solver and displacing the $\mathcal{O}(n^6)$ Kronecker approach ^{[Bartels 1972]}. Sven Hammarling's 1982 variant computed the Cholesky factor of the Lyapunov solution directly from the Schur form, preserving positive-semidefiniteness to working precision and supplying the factored Gramians that balanced-truncation model reduction consumes ^{[Hammarling 1982]}. The Bartels-Stewart method, with its Hessenberg-Schur refinement due to Golub, Nash, and Van Loan, remains the dense solver behind the Sylvester and Lyapunov routines of LAPACK and SLICOT.

Bibliography Master

@article{BartelsStewart1972,
  author  = {Bartels, Richard H. and Stewart, G. W.},
  title   = {Solution of the Matrix Equation {$AX+XB=C$}},
  journal = {Communications of the ACM},
  volume  = {15},
  number  = {9},
  year    = {1972},
  pages   = {820--826}
}

@article{Hammarling1982,
  author  = {Hammarling, Sven J.},
  title   = {Numerical Solution of the Stable Non-negative Definite Lyapunov Equation},
  journal = {IMA Journal of Numerical Analysis},
  volume  = {2},
  number  = {3},
  year    = {1982},
  pages   = {303--323}
}

@article{GolubNashVanLoan1979,
  author  = {Golub, Gene H. and Nash, Stephen and Van Loan, Charles F.},
  title   = {A {H}essenberg-{S}chur Method for the Problem {$AX+XB=C$}},
  journal = {IEEE Transactions on Automatic Control},
  volume  = {24},
  number  = {6},
  year    = {1979},
  pages   = {909--913}
}

@article{Sylvester1884,
  author  = {Sylvester, James Joseph},
  title   = {Sur l'\'equation en matrices $px = xq$},
  journal = {Comptes Rendus de l'Acad\'emie des Sciences},
  volume  = {99},
  year    = {1884},
  pages   = {67--71, 115--116}
}

@book{Lyapunov1892,
  author    = {Lyapunov, Aleksandr M.},
  title     = {The General Problem of the Stability of Motion},
  publisher = {Kharkov Mathematical Society (transl. Taylor \& Francis, 1992)},
  year      = {1892}
}

@book{LancasterTismenetsky1985,
  author    = {Lancaster, Peter and Tismenetsky, Miron},
  title     = {The Theory of Matrices},
  edition   = {2nd},
  publisher = {Academic Press},
  address   = {Orlando},
  year      = {1985}
}

@book{HornJohnson1991,
  author    = {Horn, Roger A. and Johnson, Charles R.},
  title     = {Topics in Matrix Analysis},
  publisher = {Cambridge University Press},
  address   = {Cambridge},
  year      = {1991}
}

@book{GolubVanLoan2013,
  author    = {Golub, Gene H. and Van Loan, Charles F.},
  title     = {Matrix Computations},
  edition   = {4th},
  publisher = {Johns Hopkins University Press},
  address   = {Baltimore},
  year      = {2013}
}