Bartels–Stewart algorithm

In numerical linear algebra, the Bartels–Stewart algorithm is used to numerically solve the Sylvester matrix equation ${\displaystyle AX-XB=C}${\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}${\displaystyle A} and ${\displaystyle B}${\displaystyle B} to transform ${\displaystyle AX-XB=C}${\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}${\displaystyle X} is of small to moderate size.

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

1.Compute the real Schur decompositions

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

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

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

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

3. Solve the simplified system ${\displaystyle RY-YS^{T}=F}${\displaystyle RY-YS^{T}=F}, where ${\displaystyle Y=U^{T}XV}${\displaystyle Y=U^{T}XV}. This can be done using forward substitution on the blocks. Specifically, if ${\displaystyle s_{k-1,k}=0}${\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},}${\displaystyle (R-s_{kk}I)y_{k}=f_{k}+\sum _{j=k+1}^{n}s_{kj}y_{j},}

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

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

Using the QR algorithm, the real Schur decompositions in step 1 require approximately ${\displaystyle 10(m^{3}+n^{3})}${\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})}${\displaystyle 10(m^{3}+n^{3})+2.5(mn^{2}+nm^{2})}.^[2]

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

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

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.

For large systems, the ${\displaystyle {\mathcal {O}}(m^{3}+n^{3})}${\displaystyle {\mathcal {O}}(m^{3}+n^{3})} cost of the Bartels–Stewart algorithm can be prohibitive. When ${\displaystyle A}${\displaystyle A} and ${\displaystyle B}${\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}${\displaystyle X} when solving ${\displaystyle AX-XB=C}${\displaystyle AX-XB=C}.

• Algorithms
• Control theory
• Matrices (mathematics)
• Numerical linear algebra

• Articles with short description
• Short description matches Wikidata