Bartels–Stewart algorithm: Difference between revisions

Revision as of 00:37, 21 September 2018

Revision as of 12:12, 21 September 2018

Line 3:

'''Bartels-Stewart Algorithm'''

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

== The algorithm ==

Line 12:

<math>R = U^TAU</math>,

<math>S = V^TBVTB^TV</math>.

The matrices <math>R^T</math> and <math>S</math> are block-upper triangular matrices, with square blocks of size no greater than <math>2</math>.

Line 18:

2. Set <math>F = U^TCV</math>.

3. Solve the simplified system <math>RY - YS^T = F</math>, where <math>Y = U^TXV^T</math>. This can be done using forward substitution on the blocks. Specifically, one has that if <math>s_{k-1, k} = 0</math>, then

<math>(R - s_{kk}I)y_k = f_{k} + \sum_{j = k+1}^n s_{kj}y_j</math>,

Line 27:

=== Computational cost ===

Using the [[QR algorithm]], the [[Schur decomposition| real Schur decompositions]] in step 1 requiresrequire approximately <math>10(m^3 + n^3)</math> flops, so that the overall computational cost is <math>10(m^3 + n^3) + 2.5(mn^2 + nm^2)</math>.

=== Simplifications and special cases ===

If <math>A</math> and <math>B</math> are symmetric matrices, then <math>R</math> and <math>S</math> from step 1 are diagonal. The solution matrix <math>Y</math> in step 3 can therefore be expressed explicitly as <math>Y_{jk} = F_{jk}/(R_{jj} - R_{kk})</math>.

In the special case where <math>AX + XAB=-A^T =C</math>, withand <math>C</math> is symmetric, the solution <math>X</math> will also be symmetric. This symmetry can be exploited so that <math>Y</math> is found more efficiently in step 3 of the algorithm.

If <math>A</math> and <math>B</math> are both symmetric matrices, then <math>R</math> and <math>S</math> from step 1 are diagonal. The solution matrix <math>Y</math> in step 3 can therefore be expressed explicitly as <math>Y_{jk} = F_{jk}/(R_{jj} - R_S_{kk})</math>.

== The Hessenberg-Schur algorithm ==

The Hessenberg-Schur algorithm replaces the decomposition <math>R = U^TAU</math> in step 1 with the decomposition <math>H = Q^TAQ</math>, where <math>H</math> is an [[Hessenberg matrix| upper-Hessenberg matrix]] and <math>S</math> is block-upper triangular. This leads to a system of the form <math> HY +- YS^T = F</math> that can be solved using forward substitution. The advantage of this approach is that <math>H = Q^TAQ</math> can be found using [[Householder transformation| Householder reflections]] at a cost of <math>(5/3)m^3</math> flops, compared to the <math>10m^3</math> flops required to compute the real Schur decomposition of <math>A</math>.

== 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 <math>\mathcal{O}(m^3 + n^3)</math> cost of the Bartels-Stewart algorithm can be prohibitive. When <math>A</math>and <math>B</math> 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 <math>X</math> when solving <math>AX-XB = C</math>. This is important when, for instance, <math>X</math> is too large to be stored in memory explicitly.

== References ==