Bartels–Stewart algorithm

In numerical linear algebra, the Bartels–Stewart algorithm is used to numerically solve the Sylvester matrix equation ${\displaystyle AX-XB=C}$. Developed by R.H. Bartels and G.W. Stewart in 1971,^[1] it was the first numerically stable method that could be systematically applied to solve such equations. The algorithm works by using the real Schur decompositions of ${\displaystyle A}$ and ${\displaystyle B}$ to transform ${\displaystyle AX-XB=C}$ into a triangular system that can then be solved using forward or backward substitution. In 1979, G. Golub, C. Van Loan and S. Nash introduced an improved version of the algorithm,^[2] known as the Hessenberg–Schur algorithm. It remains a standard approach for solving Sylvester equations when ${\displaystyle X}$ is of small to moderate size.

The algorithm

Let ${\displaystyle X,C\in \mathbb {R} ^{m\times n}}$, and assume that the eigenvalues of ${\displaystyle A}$ are distinct from the eigenvalues of ${\displaystyle B}$. Then, the matrix equation ${\displaystyle AX-XB=C}$ has a unique solution. The Bartels–Stewart algorithm computes ${\displaystyle X}$ by applying the following steps:^[2]

1.Compute the real Schur decompositions

${\displaystyle R=U^{T}AU,}$

${\displaystyle S=V^{T}B^{T}V.}$

The matrices ${\displaystyle R}$ and ${\displaystyle S}$ are block-upper triangular matrices, with diagonal blocks of size ${\displaystyle 1\times 1}$ or ${\displaystyle 2\times 2}$.

2. Set ${\displaystyle F=U^{T}CV.}$

3. Solve the simplified system ${\displaystyle RY-YS^{T}=F}$, where ${\displaystyle Y=U^{T}XV}$. This can be done using forward substitution on the blocks. Specifically, if ${\displaystyle s_{k-1,k}=0}$, then

${\displaystyle (R-s_{kk}I)y_{k}=f_{k}+\sum _{j=k+1}^{n}s_{kj}y_{j},}$

where ${\displaystyle y_{k}}$ is the ${\displaystyle k}$th column of ${\displaystyle Y}$. When ${\displaystyle s_{k-1,k}\neq 0}$, columns ${\displaystyle [y_{k-1}\mid y_{k}]}$ should be concatenated and solved for simultaneously.

4. Set ${\displaystyle X=UYV^{T}.}$

Computational cost

Using the QR algorithm, the real Schur decompositions in step 1 require approximately ${\displaystyle 10(m^{3}+n^{3})}$ flops, so that the overall computational cost is ${\displaystyle 10(m^{3}+n^{3})+2.5(mn^{2}+nm^{2})}$.^[2]

Simplifications and special cases

In the special case where ${\displaystyle B=-A^{T}}$ and ${\displaystyle C}$ is symmetric, the solution ${\displaystyle X}$ will also be symmetric. This symmetry can be exploited so that ${\displaystyle Y}$ is found more efficiently in step 3 of the algorithm.^[1]

The Hessenberg–Schur algorithm

The Hessenberg–Schur algorithm^[2] replaces the decomposition ${\displaystyle R=U^{T}AU}$ in step 1 with the decomposition ${\displaystyle H=Q^{T}AQ}$, where ${\displaystyle H}$ is an upper-Hessenberg matrix. This leads to a system of the form ${\displaystyle HY-YS^{T}=F}$ that can be solved using forward substitution. The advantage of this approach is that ${\displaystyle H=Q^{T}AQ}$ can be found using Householder reflections at a cost of ${\displaystyle (5/3)m^{3}}$ flops, compared to the ${\displaystyle 10m^{3}}$ flops required to compute the real Schur decomposition of ${\displaystyle A}$.

Software and implementation

The subroutines required for the Hessenberg-Schur variant of the Bartels–Stewart algorithm are implemented in the SLICOT library. These are used in the MATLAB control system toolbox.

Alternative approaches

For large systems, the ${\displaystyle {\mathcal {O}}(m^{3}+n^{3})}$ cost of the Bartels–Stewart algorithm can be prohibitive. When ${\displaystyle A}$ and ${\displaystyle B}$ are sparse or structured, so that linear solves and matrix vector multiplies involving them are efficient, iterative algorithms can potentially perform better. These include projection-based methods, which use Krylov subspace iterations, methods based on the alternating direction implicit (ADI) iteration, and hybridizations that involve both projection and ADI.^[3] Iterative methods can also be used to directly construct low rank approximations to ${\displaystyle X}$ when solving ${\displaystyle AX-XB=C}$.

References

1. Bartels, R. H.; Stewart, G. W. (1972). "Solution of the matrix equation AX + XB = C [F4]". Communications of the ACM. 15 (9): 820–826. doi:10.1145/361573.361582. ISSN 0001-0782.
2. Golub, G.; Nash, S.; Loan, C. Van (1979). "A Hessenberg–Schur method for the problem AX + XB= C". IEEE Transactions on Automatic Control. 24 (6): 909–913. doi:10.1109/TAC.1979.1102170. hdl:1813/7472. ISSN 0018-9286.
3. Simoncini, V. (2016). "Computational Methods for Linear Matrix Equations". SIAM Review. 58 (3): 377–441. doi:10.1137/130912839. hdl:11585/586011. ISSN 0036-1445. S2CID 17271167.