Alternating-direction implicit method

The method was developed at Humble Oil in the mid-1950s by Jim Douglas Jr, Henry Rachford, and Don Peaceman.

ADI for matrix equations

The method

The ADI method is a two step iteration process that alternately updates the column and row spaces of an approximate solution to ${\displaystyle AX-XB=C}$. One ADI iteration consists of the following steps:

1. Solve for ${\displaystyle X^{(j+1/2)}}$, where ${\displaystyle \left(A-{\beta }_{j+1}I\right)X^{(j+1/2)}=X^{(j)}\left(B-{\beta }_{j+1}I\right)+C.}$

2. Solve for ${\displaystyle X^{(j+1)}}$, where ${\displaystyle X^{(j+1)}\left(B-{\alpha }_{j+1}I\right)=\left(A-{\alpha }_{j+1}I\right)X^{(j+1/2)}-C}$.

The numbers ${\displaystyle ({\alpha }_{j+1},{\beta }_{j+1})}$ are called shift parameters, and convergence depends strongly on the choice of these parameters. To perform ${\displaystyle K}$ iterations of ADI, an initial guess ${\displaystyle X^{(0)}}$ is required, as well as ${\displaystyle K}$ shift parameters, ${\displaystyle \{({\alpha }_j,{\beta }_j){\}}_{j=1}^K}$.

When to use ADI

If ${\displaystyle A\in {\mathbb{C}}^{m\times m}}$ and ${\displaystyle B\in {\mathbb{C}}^{n\times n}}$, then ${\displaystyle AX-XB=C}$ can be solved directly in ${\displaystyle \mathcal{O}(m^3+n^3)}$ using the Bartels-Stewart method. It is therefore only beneficial to use ADI when matrix-vector multiplication and linear solves involving ${\displaystyle A}$ and ${\displaystyle B}$ can be applied cheaply.

The equation ${\displaystyle AX-XB=C}$ has a unique solution if and only if ${\displaystyle \sigma (A)\cap \sigma (B)=\mathrm{\varnothing }}$, where ${\displaystyle \sigma (M)}$ is the spectrum of ${\displaystyle M}$. However, the ADI method performs especially well when ${\displaystyle \sigma (A)}$ and ${\displaystyle \sigma (B)}$ are well-separated, and ${\displaystyle A}$ and ${\displaystyle B}$ are normal matrices. These assumptions are met, for example, by the Lyapunov equation ${\displaystyle AX+X{A}^{\ast }=C}$ when ${\displaystyle A}$ is positive definite. Under these assumptions, near-optimal shift parameters are known for several choices of ${\displaystyle A}$ and ${\displaystyle B}$. Additionally, a priori error bounds can be computed, thereby eliminating the need to monitor the residual error in implementation.

The ADI method can still be applied when the above assumptions are not met. The use of suboptimal shift parameters may adversely affect convergence, and convergence is also affected by the non-normality of ${\displaystyle A}$ or ${\displaystyle B}$ (sometimes advantageously). Krylov subspace methods, such as the Rational Krylov Subspace Method, are observed to typically converge more rapidly than ADI in this setting, and this has led to the development of hybrid ADI-projection methods.

Shift-parameter selection and the ADI error equation

The problem of finding good shift parameters is nontrivial. This problem can be understood by examining the ADI error equation. After ${\displaystyle K}$ iterations, the error is given by

${\displaystyle X-X^{(K)}=\prod _{j=1}^K\frac{(A-{\alpha }_{j}I)}{(A-{\beta }_{j}I)}\left(X-X^{(0)}\right)\prod _{j=1}^K\frac{(B-{\beta }_{j}I)}{(B-{\alpha }_{j}I)}.}$

Choosing ${\displaystyle X^{(0)}=0}$ results in the following bound on the relative error:

${\displaystyle \frac{{\Vert X-X^{(K)}\Vert }_2}{\Vert X{\Vert }_2}\le \Vert r_K(A){\Vert }_2\Vert r_K(B{)}^{-1}{\Vert }_2,r_K(M)=\prod _{j=1}^K\frac{(M-{\alpha }_{j}I)}{(M-{\beta }_{j}I)}.}$

where ${\displaystyle \Vert \cdot {\Vert }_2}$ is the operator norm. The ideal set of shift parameters ${\displaystyle \{({\alpha }_j,{\beta }_j){\}}_{j=1}^K}$ defines a rational function ${\displaystyle r_K}$ that minimizes the quantity ${\displaystyle \Vert r_K(A){\Vert }_2\Vert r_K(B{)}^{-1}{\Vert }_2}$. If ${\displaystyle A}$ and ${\displaystyle B}$ are normal matrices and have eigendecompositions ${\displaystyle A=V_A{\mathrm{\Lambda }}_A{V}_A^{\ast }}$ and ${\displaystyle B=V_B{\mathrm{\Lambda }}_B{V}_B^{\ast }}$, then

${\displaystyle \Vert r_K(A){\Vert }_2\Vert r_K(B{)}^{-1}{\Vert }_2=\Vert r_K({\mathrm{\Lambda }}_A){\Vert }_2\Vert r_K({\mathrm{\Lambda }}_B{)}^{-1}{\Vert }_2}$.

Near-optimal shift parameters

Near-optimal shift parameters are known in certain cases, such as when ${\displaystyle {\mathrm{\Lambda }}_A\subset [a,b]}$ and ${\displaystyle {\mathrm{\Lambda }}_B\subset [c,d]}$, where ${\displaystyle [a,b]}$ and ${\displaystyle [c,d]}$ are disjoint intervals on the real line. The Lyapunov equation ${\displaystyle AX+X{A}^{\ast }=C}$, for example, satisfies these assumptions when ${\displaystyle A}$ is positive definite. In this case, the shift parameters can be expressed in closed form using elliptic integrals, and can easily be computed numerically.

More generally, if closed, disjoint sets ${\displaystyle E}$ and ${\displaystyle F}$, where ${\displaystyle {\mathrm{\Lambda }}_A\subset E}$ and ${\displaystyle {\mathrm{\Lambda }}_B\subset F}$, are known, the optimal shift parameter selection problem is approximately solved by finding an extremal rational function that attains the value

${\displaystyle Z_K(E,F):=\underset{r}{inf}\frac{\underset{z\in E}{sup}|r(z)|}{\underset{z\in F}{inf}|r(z)|},}$

where the infimum is taken over all rational functions of degree ${\displaystyle (K,K)}$. This approximation problem is related to several results in potential theory, and was solved by Zolotarev in 1877 for ${\displaystyle E}$ = [a, b] and ${\displaystyle F=-E.}$ The solution is also known when ${\displaystyle E}$ and ${\displaystyle F}$ are disjoint disks in the complex plane.

Heuristic shift-parameter strategies

When less is known about ${\displaystyle \sigma (A)}$ and ${\displaystyle \sigma (B)}$, or when ${\displaystyle A}$ or ${\displaystyle B}$ are non-normal matrices, it may not be possible to find near-optimal shift parameters. In this setting, a variety of strategies for generating good shift parameters can be used. These include strategies based on asymptotic results in potential theory, using the Ritz values of the matrices ${\displaystyle A}$, ${\displaystyle A^{-1}}$, ${\displaystyle B}$, and ${\displaystyle B^{-1}}$ to formulate a greedy approach, and cyclic methods, where the same small collection of shift parameters are reused until a convergence tolerance is met. When the same shift parameter is used at every iteration, ADI is equivalent to an algorithm called Smith's method.

Factored ADI

In many applications, ${\displaystyle A}$ and ${\displaystyle B}$ are very large, sparse matrices, and ${\displaystyle C}$ can be factored as ${\displaystyle C=C_1{C}_2^{\ast }}$, where ${\displaystyle C_1\in {\mathbb{C}}^{m\times r},C_2\in {\mathbb{C}}^{n\times r}}$, with ${\displaystyle r=1,2}$. In such a setting, it may not be feasible to store the potentially dense matrix ${\displaystyle X}$ explicitly. A variant of ADI, called factored ADI, can be used to compute ${\displaystyle Z{Y}^{\ast }}$, where ${\displaystyle X\approx Z{Y}^{\ast }}$. The effectiveness of factored ADI depends on whether ${\displaystyle X}$ is well-approximated by a low rank matrix. This is known to be true under various assumptions about ${\displaystyle A}$ and ${\displaystyle B}$.