# Alternating direction implicit method

In numerical linear algebra, the Alternating Direction Implicit (ADI) method is an iterative method used to solve Sylvester matrix equations. It is a popular method for solving the large matrix equations that arise in systems theory and control,[1] and can be formulated to construct solutions in a memory-efficient, factored form.[2][3] It is also used to numerically solve parabolic and elliptic partial differential equations, and is a classic method used for modeling heat conduction and solving the diffusion equation in two or more dimensions.[4] It is an example of an operator splitting method.[5]

### 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:[6]

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.[7][8] 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}}$.

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.[9] 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)=\emptyset }$, where ${\displaystyle \sigma (M)}$ is the spectrum of ${\displaystyle M}$.[1] 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+XA^{*}=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}$.[7][8] 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,[1] and convergence is also affected by the non-normality of ${\displaystyle A}$ or ${\displaystyle B}$ (sometimes advantageously).[10] Krylov subspace methods, such as the Rational Krylov Subspace Method,[11] are observed to typically converge more rapidly than ADI in this setting,[1][3] and this has led to the development of hybrid ADI-projection methods.[3]

### 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 {\left\|X-X^{(K)}\right\|_{2}}{\|X\|_{2}}}\leq \|r_{K}(A)\|_{2}\|r_{K}(B)^{-1}\|_{2},\quad r_{K}(M)=\prod _{j=1}^{K}{\frac {(M-\alpha _{j}I)}{(M-\beta _{j}I)}}.}$

where ${\displaystyle \|\cdot \|_{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 \|r_{K}(A)\|_{2}\|r_{K}(B)^{-1}\|_{2}}$. If ${\displaystyle A}$ and ${\displaystyle B}$ are normal matrices and have eigendecompositions ${\displaystyle A=V_{A}\Lambda _{A}V_{A}^{*}}$ and ${\displaystyle B=V_{B}\Lambda _{B}V_{B}^{*}}$, then

${\displaystyle \|r_{K}(A)\|_{2}\|r_{K}(B)^{-1}\|_{2}=\|r_{K}(\Lambda _{A})\|_{2}\|r_{K}(\Lambda _{B})^{-1}\|_{2}}$.

#### Near-optimal shift parameters

Near-optimal shift parameters are known in certain cases, such as when ${\displaystyle \Lambda _{A}\subset [a,b]}$ and ${\displaystyle \Lambda _{B}\subset [c,d]}$, where ${\displaystyle [a,b]}$ and ${\displaystyle [c,d]}$ are disjoint intervals on the real line.[7][8] The Lyapunov equation ${\displaystyle AX+XA^{*}=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 \Lambda _{A}\subset E}$ and ${\displaystyle \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):=\inf _{r}{\frac {\sup _{z\in E}|r(z)|}{\inf _{z\in F}|r(z)|}},}$

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

#### 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,[16] using the Ritz values of the matrices ${\displaystyle A}$, ${\displaystyle A^{-1}}$, ${\displaystyle B}$, and ${\displaystyle B^{-1}}$ to formulate a greedy approach,[17] and cyclic methods, where the same small collection of shift parameters are reused until a convergence tolerance is met.[17][10] When the same shift parameter is used at every iteration, ADI is equivalent to an algorithm called Smith's method.[18]

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}^{*}}$, where ${\displaystyle C_{1}\in \mathbb {C} ^{m\times r},C_{2}\in \mathbb {C} ^{n\times r}}$, with ${\displaystyle r=1,2}$.[1] 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,[3][2] can be used to compute ${\displaystyle ZY^{*}}$, where ${\displaystyle X\approx ZY^{*}}$. 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}$.[10][8]

Historically, the ADI method was developed to solve the 2D diffusion equation on a square domain using finite differences.[4] Unlike ADI for matrix equations, ADI for parabolic equations does not require the selection of shift parameters, since the shift appearing in each iteration is determined by parameters such as the timestep, diffusion coefficient, and grid spacing. The connection to ADI on matrix equations can be observed when one considers the action of the ADI iteration on the system at steady state.

### Example: 2D diffusion equation

Stencil figure for the alternating direction implicit method in finite difference equations

The traditional method for solving the heat conduction equation numerically is the Crank–Nicolson method. This method results in a very complicated set of equations in multiple dimensions, which are costly to solve. The advantage of the ADI method is that the equations that have to be solved in each step have a simpler structure and can be solved efficiently with the tridiagonal matrix algorithm.

Consider the linear diffusion equation in two dimensions,

${\displaystyle {\partial u \over \partial t}=\left({\partial ^{2}u \over \partial x^{2}}+{\partial ^{2}u \over \partial y^{2}}\right)=(u_{xx}+u_{yy})}$

The implicit Crank–Nicolson method produces the following finite difference equation:

${\displaystyle {u_{ij}^{n+1}-u_{ij}^{n} \over \Delta t}={1 \over 2(\Delta x)^{2}}\left(\delta _{x}^{2}+\delta _{y}^{2}\right)\left(u_{ij}^{n+1}+u_{ij}^{n}\right)}$

where:

${\displaystyle \Delta x=\Delta y}$

and ${\displaystyle \delta _{p}^{2}}$ is the central second difference operator for the p-th coordinate

${\displaystyle \delta _{p}^{2}u_{ij}=u_{ij+e_{p}}-2u_{ij}+u_{ij-e_{p}}}$

with ${\displaystyle e_{p}=10}$ or ${\displaystyle 01}$ for ${\displaystyle p=x}$ or ${\displaystyle y}$ respectively (and ${\displaystyle ij}$ a shorthand for lattice points ${\displaystyle (i,j)}$).

After performing a stability analysis, it can be shown that this method will be stable for any ${\displaystyle \Delta t}$.

A disadvantage of the Crank–Nicolson method is that the matrix in the above equation is banded with a band width that is generally quite large. This makes direct solution of the system of linear equations quite costly (although efficient approximate solutions exist, for example use of the conjugate gradient method preconditioned with incomplete Cholesky factorization).

The idea behind the ADI method is to split the finite difference equations into two, one with the x-derivative taken implicitly and the next with the y-derivative taken implicitly,

${\displaystyle {u_{ij}^{n+1/2}-u_{ij}^{n} \over \Delta t/2}={\left(\delta _{x}^{2}u_{ij}^{n+1/2}+\delta _{y}^{2}u_{ij}^{n}\right) \over \Delta x^{2}}}$
${\displaystyle {u_{ij}^{n+1}-u_{ij}^{n+1/2} \over \Delta t/2}={\left(\delta _{x}^{2}u_{ij}^{n+1/2}+\delta _{y}^{2}u_{ij}^{n+1}\right) \over \Delta y^{2}}}$

The system of equations involved is symmetric and tridiagonal (banded with bandwidth 3), and is typically solved using tridiagonal matrix algorithm.

It can be shown that this method is unconditionally stable and second order in time and space.[19] There are more refined ADI methods such as the methods of Douglas,[20] or the f-factor method[21] which can be used for three or more dimensions.

### Generalizations

The usage of the ADI method as an operator splitting scheme can be generalized. That is, we may consider general evolution equations

${\displaystyle {\dot {u}}=F_{1}u+F_{2}u,}$

where ${\displaystyle F_{1}}$ and ${\displaystyle F_{2}}$ are (possibly nonlinear) operators defined on a Banach space.[22][23] In the diffusion example above we have ${\displaystyle F_{1}={\partial ^{2} \over \partial x^{2}}}$ and ${\displaystyle F_{2}={\partial ^{2} \over \partial y^{2}}}$.

## References

1. Simoncini, V. (2016). "Computational Methods for Linear Matrix Equations". SIAM Review. 58 (3): 377–441. doi:10.1137/130912839. ISSN 0036-1445.
2. ^ a b Li, Jing-Rebecca; White, Jacob (2002). "Low Rank Solution of Lyapunov Equations". SIAM Journal on Matrix Analysis and Applications. 24 (1): 260–280. doi:10.1137/s0895479801384937. ISSN 0895-4798.
3. ^ a b c d Benner, Peter; Li, Ren-Cang; Truhar, Ninoslav (2009). "On the ADI method for Sylvester equations". Journal of Computational and Applied Mathematics. 233 (4): 1035–1045. Bibcode:2009JCoAM.233.1035B. doi:10.1016/j.cam.2009.08.108. ISSN 0377-0427.
4. ^ a b Peaceman, D. W.; Rachford Jr., H. H. (1955), "The numerical solution of parabolic and elliptic differential equations", Journal of the Society for Industrial and Applied Mathematics, 3 (1): 28–41, doi:10.1137/0103003, MR 0071874.
5. ^ *Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). "Section 20.3.3. Operator Splitting Methods Generally". Numerical Recipes: The Art of Scientific Computing (3rd ed.). New York: Cambridge University Press. ISBN 978-0-521-88068-8.
6. ^ Wachspress, Eugene L. (2008). "Trail to a Lyapunov equation solver". Computers & Mathematics with Applications. 55 (8): 1653–1659. doi:10.1016/j.camwa.2007.04.048. ISSN 0898-1221.
7. ^ a b c Lu, An; Wachspress, E.L. (1991). "Solution of Lyapunov equations by alternating direction implicit iteration". Computers & Mathematics with Applications. 21 (9): 43–58. doi:10.1016/0898-1221(91)90124-m. ISSN 0898-1221.
8. Beckermann, Bernhard; Townsend, Alex (2017). "On the Singular Values of Matrices with Displacement Structure". SIAM Journal on Matrix Analysis and Applications. 38 (4): 1227–1248. arXiv:1609.09494. doi:10.1137/16m1096426. ISSN 0895-4798.
9. ^ Golub, G. (1989). Matrix computations. Van Loan, C (Fourth ed.). Baltimore: Johns Hopkins University. ISBN 1421407949. OCLC 824733531.
10. ^ a b c Sabino, J (2007). "Solution of large-scale Lyapunov equations via the block modified Smith method". PhD diss., Rice Univ.
11. ^ Druskin, V.; Simoncini, V. (2011). "Adaptive rational Krylov subspaces for large-scale dynamical systems". Systems & Control Letters. 60 (8): 546–560. doi:10.1016/j.sysconle.2011.04.013. ISSN 0167-6911.
12. ^ 1944-, Saff, E. B.,. Logarithmic potentials with external fields. Totik, V. Berlin. ISBN 9783662033296. OCLC 883382758.
13. ^ Gonchar, A.A. (1969). "Zolotarev problems connected with rational functions". Mat. Sb. (N.S.). 78 (120):4: 640–654.
14. ^ Zolotarev, D.I. (1877). "Application of elliptic functions to questions of functions deviating least and most from zero". Zap. Imp. Akad. Nauk. St. Petersburg. 30: 1–59.
15. ^ Starke, Gerhard (July 1992). "Near-circularity for the rational Zolotarev problem in the complex plane". Journal of Approximation Theory. 70 (1): 115–130. doi:10.1016/0021-9045(92)90059-w. ISSN 0021-9045.
16. ^ Starke, Gerhard (June 1993). "Fejér-Walsh points for rational functions and their use in the ADI iterative method". Journal of Computational and Applied Mathematics. 46 (1–2): 129–141. doi:10.1016/0377-0427(93)90291-i. ISSN 0377-0427.
17. ^ a b Penzl, Thilo (January 1999). "A Cyclic Low-Rank Smith Method for Large Sparse Lyapunov Equations". SIAM Journal on Scientific Computing. 21 (4): 1401–1418. doi:10.1137/s1064827598347666. ISSN 1064-8275.
18. ^ Smith, R. A. (January 1968). "Matrix Equation XA + BX = C". SIAM Journal on Applied Mathematics. 16 (1): 198–201. doi:10.1137/0116017. ISSN 0036-1399.
19. ^ Douglas, Jr., J. (1955), "On the numerical integration of uxx+ uyy= ut by implicit methods", Journal of the Society of Industrial and Applied Mathematics, 3: 42–65, MR 0071875.
20. ^ Douglas Jr., Jim (1962), "Alternating direction methods for three space variables", Numerische Mathematik, 4 (1): 41–63, doi:10.1007/BF01386295, ISSN 0029-599X.
21. ^ Chang, M. J.; Chow, L. C.; Chang, W. S. (1991), "Improved alternating-direction implicit method for solving transient three-dimensional heat diffusion problems", Numerical Heat Transfer, Part B: Fundamentals, 19 (1): 69–84, doi:10.1080/10407799108944957, ISSN 1040-7790.
22. ^ Hundsdorfer, Willem; Verwer, Jan (2003). Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations. Berlin, Heidelberg: Springer Berlin Heidelberg. ISBN 978-3-662-09017-6.
23. ^ Lions, P. L.; Mercier, B. (December 1979). "Splitting Algorithms for the Sum of Two Nonlinear Operators". SIAM Journal on Numerical Analysis. 16 (6): 964–979. doi:10.1137/0716071.