Innumerical linear algebra, thealternating-direction implicit (ADI) method is an iterative method used to solveSylvester matrix equations. It is a popular method for solving the large matrix equations that arise insystems theory andcontrol,^[1] and can be formulated to construct solutions in a memory-efficient, factored form.^[2][3] It is also used to numerically solveparabolic andelliptic partial differential equations, and is a classic method used for modelingheat conduction and solving thediffusion equation in two or more dimensions.^[4] It is an example of anoperator splitting method.^[5]

The method was developed atHumble Oil in the mid-1950s by Jim Douglas Jr, Henry Rachford, and Don Peaceman.^[6]

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:^[7]

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.^[8][9] 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 theBartels-Stewart method.^[10] 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 thespectrum 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}$ arenormal matrices. These assumptions are met, for example, by the Lyapunov equation${\displaystyle AX+X{A}^{\ast }=C}$ when${\displaystyle A}$ ispositive definite. Under these assumptions, near-optimal shift parameters are known for several choices of${\displaystyle A}$ and${\displaystyle B}$.^[8][9] 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).^[11]Krylov subspace methods, such as the Rational Krylov Subspace Method,^[12] 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]

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 theoperator norm. The ideal set of shift parameters${\displaystyle \{({\alpha }_j,{\beta }_j){\}}_{j=1}^K}$ defines arational 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}$ arenormal matrices and haveeigendecompositions${\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 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.^[8][9] TheLyapunov equation${\displaystyle AX+X{A}^{\ast }=C}$, for example, satisfies these assumptions when${\displaystyle A}$ ispositive definite. In this case, the shift parameters can be expressed in closed form usingelliptic 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)}$.^[9] This approximation problem is related to several results inpotential theory,^[13][14] and was solved byZolotarev in 1877 for${\displaystyle E}$ = [a, b] and${\displaystyle F=-E.}$^[15] The solution is also known when${\displaystyle E}$ and${\displaystyle F}$ are disjoint disks in the complex plane.^[16]

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

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}$.^[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 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}$.^[11][9]

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.

The traditional method for solving the heat conduction equation numerically is theCrank–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 thetridiagonal matrix algorithm.

Consider the linear diffusion equation in two dimensions,

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

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

${\displaystyle \frac{u_{ij}^{n+1}-u_{ij}^n}{\mathrm{\Delta }t}=\frac{1}{2(\mathrm{\Delta }x{)}^2}\left({\delta }_x^2+{\delta }_y^2\right)\left(u_{ij}^{n+1}+u_{ij}^n\right)}$

where:

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

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

${\displaystyle {\delta }_p^2{u}_{ij}=u_{ij+e_p}-2u_{ij}+u_{ij-e_p}}$

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

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

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

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

${\displaystyle \frac{u_{ij}^{n+1/2}-u_{ij}^n}{\mathrm{\Delta }t/2}=\frac{\left({\delta }_x^2{u}_{ij}^{n+1/2}+{\delta }_y^2{u}_{ij}^n\right)}{\mathrm{\Delta }{x}^2}}$

${\displaystyle \frac{u_{ij}^{n+1}-u_{ij}^{n+1/2}}{\mathrm{\Delta }t/2}=\frac{\left({\delta }_x^2{u}_{ij}^{n+1/2}+{\delta }_y^2{u}_{ij}^{n+1}\right)}{\mathrm{\Delta }{y}^2}}$

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

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

The usage of the ADI method as anoperator 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 aBanach space.^[23][24] In the diffusion example above we have${\displaystyle F_1=\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{x}^2}}$ and${\displaystyle F_2=\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{y}^2}}$.

It is possible to simplify the conventional ADI method into Fundamental ADI method, which only has the similar operators at the left-hand sides while being operator-free at the right-hand sides. This may be regarded as the fundamental (basic) scheme of ADI method,^[25][26] with no more operator (to be reduced) at the right-hand sides, unlike most traditional implicit methods that usually consist of operators at both sides of equations. The FADI method leads to simpler, more concise and efficient update equations without degrading the accuracy of conventional ADI method.

Many classical implicit methods by Peaceman-Rachford, Douglas-Gunn, D'Yakonov, Beam-Warming, Crank-Nicolson, etc., may be simplified to fundamental implicit schemes with operator-free right-hand sides.^[26] In their fundamental forms, the FADI method of second-order temporal accuracy can be related closely to the fundamental locally one-dimensional (FLOD) method, which can be upgraded to second-order temporal accuracy, such as for three-dimensional Maxwell's equations^[27][28] incomputational electromagnetics. For two- and three-dimensional heat conduction and diffusion equations, both FADI and FLOD methods may be implemented in simpler, more efficient and stable manner compared to their conventional methods.^[29][30]

• Usadi, Adam; Dawson, Clint (March 2006)."50 Years of ADI Methods: Celebrating the Contributions of Jim Douglas, Don Peaceman, and Henry Rachford".SIAM News.39 (2). RetrievedMarch 28, 2025 – via ResearchGate.net. (Provides a review of the ADI methods and variants over the years.)

• Partial differential equations
• Numerical differential equations
• Iterative methods

• CS1: long volume value