Bartels–Stewart 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}.}$

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}$.

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})}$ 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}$.