Bartels–Stewart algorithm

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, it was the first numerically stable method that could by systematically applied to solve such equations. The algorithm works by using the 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 and back-substitutions. In 1979, G. Golub, C. Van Loan and S.Nash introduced an improved version of the algorithm, known as the Hessenberg-Schur algorithm. It remains the 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 ${\displaystyle AX-XB=C}$  has a unique solution. The Bartels-Stewart algorithm computes ${\displaystyle X}$  by applying the following steps:

1.Compute the Schur decompositions

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

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

The matrices ${\displaystyle R^{T}}$  and ${\displaystyle S}$  are block-upper triangular matrices, with square blocks of size no greater than ${\displaystyle 2}$ .

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

3. Solve the simplified system ${\displaystyle RY-YS^{T}=F}$ , where ${\displaystyle Y=U^{T}XV^{T}}$ . This can be done using forward substitution on the blocks. Specifically, one has that 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}$ .

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

Computational cost

Using the QR algorithm, the Schur decompositions in step 1 requires 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})}$ .

The Hessenberg-Schur algorithm

The Hessenberg-Schur algorithm 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 and ${\displaystyle S}$  is block-upper triangular. 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 Schur decomposition of ${\displaystyle A}$ .

Implementations and availability

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 cheap, iterative algorithms can potentially perform better. These include projection-based methods, which use Krylov subspace iterations, methods based on the alternating-implicit direction (ADI) iteration, and hybridizations that involve both projection and ADI. Iterative methods can also be used to directly construct low rank approximations to ${\displaystyle X}$  when solving ${\displaystyle AX-XB=C}$ . This is important when, for instance, ${\displaystyle X}$  is too large to be stored in memory explicitly.

References

External links

• www.example.com