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, and can be formulated to construct solutions in a memory-efficient, factored form. 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. It is an example of an operator splitting method. The method was developed at Humble Oil in the mid-1950s by Jim Douglas Jr, Henry Rachford, and Don Peaceman.

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, and can be formulated to construct solutions in a memory-efficient, factored form. 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. It is an example of an operator splitting method. The method was developed at Humble Oil in the mid-1950s by Jim Douglas Jr, Henry Rachford, and Don Peaceman. The ADI method is a two step iteration process that alternately updates the column and row spaces of an approximate solution to AX - XB = C . One ADI iteration consists of the following steps: 1. Solve for X^{(j + 1/2)} , where \left( A - \beta_{j +1} I\right) X^{(j+1/2)} = X^{(j)}\left( B - \beta_{j + 1} I \right) + C. 2. Solve for X^{(j + 1)} , where X^{(j+1)}\left( B - \alpha_{j + 1} I \right) = \left( A - \alpha_{j+1} I\right) X^{(j+1/2)} - C . The numbers (\alpha_{j+1}, \beta_{j+1}) are called shift parameters, and convergence depends strongly on the choice of these parameters. To perform K iterations of ADI, an initial guess X^{(0)} is required, as well as K shift parameters, \{ (\alpha_{j}, \beta_{j})\}_{j = 1}^{K} . If A \in \mathbb{C}^{m \times m} and B \in \mathbb{C}^{n \times n} , then AX - XB = C can be solved directly in \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 A and B can be applied cheaply. The equation AX-XB=C has a unique solution if and only if \sigma(A) \cap \sigma(B) = \emptyset , where \sigma(M) is the spectrum of M . However, the ADI method performs especially well when \sigma(A) and \sigma(B) are well-separated, and A and B are normal matrices. These assumptions are met, for example, by the Lyapunov equation AX + XA^* = C when A is positive definite. Under these assumptions, near-optimal shift parameters are known for several choices of A and 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 A or 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. 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 X^{(0)} = 0 results in the following bound on the relative error: \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 \| \cdot \|_2 is the operator norm. The ideal set of shift parameters \{ (\alpha_j, \beta_j)\}_{j = 1}^K defines a rational function r_K that minimizes the quantity \| r_K(A) \|_2 \| r_K(B)^{-1}\|_2 . If A and B are normal matrices and have eigendecompositions A = V_A\Lambda_AV_A^* and B = V_B\Lambda_BV_B^* , then \| r_K(A) \|_2 \| r_K(B)^{-1}\|_2 = \| r_K(\Lambda_A) \|_2 \| r_K(\Lambda_B)^{-1}\|_2 . Near-optimal shift parameters are known in certain cases, such as when \Lambda_A \subset [a, b] and \Lambda_B \subset [c, d] , where [a, b] and [c, d] are disjoint intervals on the real line. The Lyapunov equation AX + XA^* = C , for example, satisfies these assumptions when 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 E and F , where \Lambda_A \subset E and \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 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 (K, K) . This approximation problem is related to several results in potential theory, and was solved by Zolotarev in 1877 for E = [a, b] and F=-E. The solution is also known when E and F are disjoint disks in the complex plane. In many applications, A and B are very large, sparse matrices, and C can be factored as C = C_1C_2^* , where C_1 \in \mathbb{C}^{m \times r}, C_2 \in \mathbb{C}^{n \times r} , with r = 1, 2 . In such a setting, it may not be feasible to store the potentially dense matrix X explicitly. A variant of ADI, called factored ADI, can be used to compute ZY^* , where X \approx ZY^* . The effectiveness of factored ADI depends on whether X is well-approximated by a low rank matrix. This is known to be true under various assumptions about A and B . {\partial u\over \partial t} = \left({\partial^2 u\over \partial x^2 } + {\partial^2 u\over \partial y^2 } \right) \Delta x = \Delta y and \delta_p^2 is the central second difference operator for the p-th coordinate \delta_p^2 u_{ij}=u_{ij+e_p}-2u_{ij}+u_{ij-e_p} with e_p=(1,0) or (0,1) for p=x or y respectively (and ij a shorthand for lattice points (i,j) ). After performing a stability analysis, it can be shown that this method will be stable for any \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, {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} {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. There are more refined ADI methods such as the methods of Douglas, or the f-factor method which can be used for three or more dimensions. \dot u = F_1 u + F_2 u, where F_1 and F_2 are (possibly nonlinear) operators defined on a Banach space. In the diffusion example above we have F_1 = {\partial^2 \over \partial x^2} and F_2 = {\partial^2 \over \partial y^2} .

Further reading

(Provides a review of the ADI methods and variants over the years.)

Text derived from https://en.wikipedia.org/w/index.php?oldid=1322703346 under CC BY-SA 4.0; modifications have been made.