Bartels–Stewart algorithm explained

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 A and B to transform 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 X is of small to moderate size.

The algorithm

Let

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

1.Compute the real Schur decompositions

R = U ^{TAU,} S = V ^{TB} ^{TV.}

The matrices

R and S are block-upper triangular matrices, with diagonal blocks of size 1 x 1 or 2 x 2 .

2. Set

F = U ^{TCV.}

3. Solve the simplified system

RY - YS ^{T} = F , where Y = U ^{TXV} . This can be done using forward substitution on the blocks. Specifically, if s _{k-1,} = 0 , then(R - s _{kk} I)y _{k} = f _{k} + \sum ^{n} _{j = k +1} s _{kj} y _{j,}

where

y _{k} is the k th column of Y . When s _{k-1,} ≠ 0 , columns [ y _{k-1} \mid y _{k} ] should be concatenated and solved for simultaneously.

4. Set

X = UYV ^{T.}

Computational cost

Using the QR algorithm, the real Schur decompositions in step 1 require approximately

10 (m ^{3} + n ^{3)} flops, so that the overall computational cost is 10 (m ^{3} + n ^{3)} + 2 . 5 (mn ^{2} + nm ^{2)} .

Simplifications and special cases

In the special case where

B =- A ^{T} and C is symmetric, the solution X will also be symmetric. This symmetry can be exploited so that Y is found more efficiently in step 3 of the algorithm.

The Hessenberg–Schur algorithm

The Hessenberg–Schur algorithm replaces the decomposition

R = U ^{TAU} in step 1 with the decomposition H = Q ^{TAQ} , where H is an upper-Hessenberg matrix. This leads to a system of the form HY - YS ^{T} = F that can be solved using forward substitution. The advantage of this approach is that H = Q ^{TAQ} can be found using Householder reflections at a cost of ( 5 / 3 )m ^{3} flops, compared to the 10 m ^{3} flops required to compute the real Schur decomposition of 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

l{O}(m ^{3} + n ^{3)} cost of the Bartels–Stewart algorithm can be prohibitive. When A and 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 X when solving AX - XB = C .

Notes and 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. 10.1145/361573.361582. 0001-0782. free.
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. 10.1109/TAC.1979.1102170. 0018-9286. 1813/7472. free.
3. Simoncini. V.. 17271167. 2016. Computational Methods for Linear Matrix Equations. SIAM Review. en-US. 58. 3. 377–441. 10.1137/130912839. 0036-1445. 11585/586011. free.