Alternating-direction implicit method: Difference between revisions

Revision as of 03:32, 8 July 2019

Latest revision as of 09:33, 15 April 2025

(43 intermediate revisions by 18 users not shown)

Line 1:

In [[numerical linear algebra]], the '''Alternating Directionalternating-direction Implicitimplicit (ADI) method''' is an iterative method used to solve [[Sylvester equation|Sylvester]] matrix equations]]. It is a popular method for solving the large matrix equations that arise in [[systems theory]] and [[Control theory|control]],<ref name=":1">{{Cite journal |last=Simoncini |first=V. |s2cid=17271167 |date=2016 |title=Computational Methods for Linear Matrix Equations |journal=SIAM Review |language=en |volume=58 |issue=3 |pages=377–441 |doi=10.1137/130912839 |issn=0036-1445|hdl=11585/586011 |hdl-access=free }}</ref> and can be formulated to construct solutions in a memory-efficient, factored form.<ref name=":2">{{Cite journal |lastlast1=Li |firstfirst1=Jing-Rebecca |author1-link=Jing-Rebecca Li |last2=White |first2=Jacob |date=2002 |title=Low Rank Solution of Lyapunov Equations |journal=SIAM Journal on Matrix Analysis and Applications |language=en |volume=24 |issue=1 |pages=260–280 |doi=10.1137/s0895479801384937 |issn=0895-4798}}</ref><ref name=":3">{{Cite journal |lastlast1=Benner |firstfirst1=Peter |last2=Li |first2=Ren-Cang |last3=Truhar |first3=Ninoslav |date=2009 |title=On the ADI method for Sylvester equations |journal=Journal of Computational and Applied Mathematics |volume=233 |issue=4 |pages=1035–1045 |doi=10.1016/j.cam.2009.08.108 |issn=0377-0427 |bibcode=2009JCoAM.233.1035B |doi-access=free}}</ref> It is also used to numerically solve [[Parabolic partial differential equation|parabolic]] and [[Elliptic partial differential equation|elliptic]] partial differential equations, and is a classic method used for modeling [[heat conduction]] and solving the [[diffusion equation]] in two or more dimensions.<ref name=":0">{{Citation |title=The numerical solution of parabolic and elliptic differential equations |year=1955 |last1=Peaceman |last2=Rachford Jr. |first1=D. W. |first2=H. H. |journal=Journal of the Society for Industrial and Applied Mathematics |volume=3 |issue=1 |pages=28–41 |doi=10.1137/0103003 |mr=0071874 |hdl=10338.dmlcz/135399 |hdl-access=free}}.</ref> It is an example of an [[operator splitting]] method.<ref>*{{Cite book | last1=Press | first1=WHW. H. | last2=Teukolsky | first2=SAS. A. | last3=Vetterling | first3=WTW. T. | last4=Flannery | first4=BPB. P. | year=2007 | title=Numerical Recipes: The Art of Scientific Computing | edition=3rd | publisher=Cambridge University Press | ___location=New York | isbn=978-0-521-88068-8 | chapter=Section 20.3.3. Operator Splitting Methods Generally | chapter-url=http://apps.nrbook.com/empanel/index.html#pg=1052 | access-date=2011-08-18 | archive-date=2011-08-11 | archive-url=https://web.archive.org/web/20110811154417/http://apps.nrbook.com/empanel/index.html#pg=1052 | url-status=dead }}</ref>

</ref>

The method was developed at [[Humble Oil]] in the mid-1950s by Jim Douglas Jr, Henry Rachford, and Don Peaceman.<ref name=PurdueObit-Douglas>{{cite web|accessdate=March 25, 2025

|url=https://www.math.purdue.edu/news/2016/douglas_obit.html

|publisher=Purdue University|title=Prof. Jim Douglas

|date=May 2, 2016}}</ref>

== 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 <math>AX - XB = C</math>. One ADI iteration consists of the following steps:<ref>{{Cite journal|last=Wachspress|first=Eugene L.|date=2008|title=Trail to a Lyapunov equation solver|journal=Computers & Mathematics with Applications|volume=55|issue=8|pages=1653–1659|doi=10.1016/j.camwa.2007.04.048|issn=0898-1221|doi-access=free}}</ref><blockquote>1. Solve for <math>X^{(j + 1/2)}</math>, where <math>\left( A - \beta_{j +1} I\right) X^{(j+1/2)} = X^{(j)}\left( B - \beta_{j + 1} I \right) + C.</math> </blockquote><blockquote>2. Solve for <math> X^{(j + 1)}</math>, where <math> X^{(j+1)}\left( B - \alpha_{j + 1} I \right) = \left( A - \alpha_{j+1} I\right) X^{(j+1/2)} - C</math>.</blockquote>

The numbers <math>(\alpha_{j+1}, \beta_{j+1})</math> are called shift parameters, and convergence depends strongly on the choice of these parameters.<ref name=":4">{{Cite journal|lastlast1=Lu|firstfirst1=An|last2=Wachspress|first2=E.L.|date=1991|title=Solution of Lyapunov equations by alternating direction implicit iteration|journal=Computers & Mathematics with Applications|volume=21|issue=9|pages=43–58|doi=10.1016/0898-1221(91)90124-m|issn=0898-1221|doi-access=}}</ref><ref name=":5">{{Cite journal|lastlast1=Beckermann|firstfirst1=Bernhard|last2=Townsend|first2=Alex|date=2017|title=On the Singular Values of Matrices with Displacement Structure|journal=SIAM Journal on Matrix Analysis and Applications|language=en|volume=38|issue=4|pages=1227–1248|doi=10.1137/16m1096426|issn=0895-4798|arxiv=1609.09494|s2cid=3828461}}</ref> To perform <math>K</math> iterations of ADI, an initial guess <math>X^{(0)}</math> is required, as well as <math>K</math> shift parameters, <math>\{ (\alpha_{j}, \beta_{j})\}_{j = 1}^{K}</math>.

=== When to use ADI ===

If <math>A \in \mathbb{C}^{m \times m}</math> and <math>B \in \mathbb{C}^{n \times n}</math>, then <math> AX - XB = C</math> can be solved directly in <math> \mathcal{O}(m^3 + n^3)</math> using the [[Bartels–Stewart algorithm|Bartels-Stewart method]].<ref>{{Cite book|title=Matrix computations|lastlast1=Golub,|first1= G.|publisher=Johns Hopkins University|otherslast2=Van Loan,|first2= C|year=1989|isbn=1421407949|edition=Fourth |___location=Baltimore|pages=|oclc=824733531}}</ref> It is therefore only beneficial to use ADI when matrix-vector multiplication and linear solves involving <math>A</math> and <math>B</math> can be applied cheaply.

The equation <math> AX-XB=C</math> has a unique solution if and only if <math> \sigma(A) \cap \sigma(B) = \emptyset</math>, where <math> \sigma(M) </math> is the [[Spectrum of a matrix|spectrum]] of <math>M</math>.<ref name=":1" /> However, the ADI method performs especially well when <math>\sigma(A)</math> and <math>\sigma(B)</math> are well-separated, and <math>A</math> and <math>B</math> are [[Normal matrix|normal matrices]]. These assumptions are met, for example, by the Lyapunov equation <math>AX + XA^* = C</math> when <math>A</math> is [[Positive-definite matrix|positive definite]]. Under these assumptions, near-optimal shift parameters are known for several choices of <math>A</math> and <math>B</math>.<ref name=":4" /><ref name=":5" /> 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,<ref name=":1" /> and convergence is also affected by the non-normality of <math>A</math> or <math>B</math> (sometimes advantageously).<ref name=":6">{{Cite thesis|last=Sabino|first=J|date=2007|title=Solution of large-scale Lyapunov equations via the block modified Smith method|journal=PHDPhD Diss., Rice Univ.|volume=|pages=|hdl=1911/20641|type=Thesis}}</ref> [[Krylov subspace]] methods, such as the Rational Krylov Subspace Method,<ref>{{Cite journal|lastlast1=Druskin|firstfirst1=V.|last2=Simoncini|first2=V.|date=2011|title=Adaptive rational Krylov subspaces for large-scale dynamical systems|journal=Systems & Control Letters|volume=60|issue=8|pages=546–560|doi=10.1016/j.sysconle.2011.04.013|issn=0167-6911}}</ref> are observed to typically converge more rapidly than ADI in this setting,<ref name=":1" /><ref name=":3" /> and this has led to the development of hybrid ADI-projection methods.<ref name=":3" />

=== Shift-parameter Parameter Selectionselection 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 <math>K</math> iterations, the error is given by

Line 37 ⟶ 41:

</math>

where the infimum is taken over all rational functions of degree <math>(K, K)</math>.<ref name=":5" /> This approximation problem is related to several results in [[potential theory]],<ref>{{Cite book|title=Logarithmic potentials with external fields|lastlast1=1944-Saff|firstfirst1=Saff, E. B.|otherslast2=Totik, |first2=V.|publisher= Springer Science & Business Media |isbn=9783662033296|___location=Berlin|oclc=883382758|date = 2013-11-11}}</ref><ref>{{Cite journal|last=Gonchar|first=A.A.|date=1969|title=Zolotarev problems connected with rational functions|url=|journal=Mat. Sb.Mathematics (N.S.)of the USSR-Sbornik|volume=78 (120):47|issue=4|pages=640–654623–635|bibcode=1969SbMat...7..623G|doi=10.1070/SM1969v007n04ABEH001107}}</ref> and was solved by [[Yegor Ivanovich Zolotarev|Zolotarev]] in 1877 for <math>E</math> = [a, b] and <math> F=-E.</math><ref>{{Cite journal|last=Zolotarev|first=D.I.|date=1877|title=Application of elliptic functions to questions of functions deviating least and most from zero|url=|journal=Zap. Imp. Akad. Nauk. St. Petersburg|volume=30|pages=1–59|via=}}</ref> The solution is also known when <math> E</math> and <math> F</math> are disjoint disks in the complex plane.<ref>{{Cite journal|last=Starke|first=Gerhard|date=July 1992|title=Near-circularity for the rational Zolotarev problem in the complex plane|journal=Journal of Approximation Theory|volume=70|issue=1|pages=115–130|doi=10.1016/0021-9045(92)90059-w|issn=0021-9045|doi-access=free}}</ref>

==== Heuristic shift -parameter strategies ====

When less is known about <math>\sigma(A)</math> and <math>\sigma(B)</math>, or when <math>A</math> or <math>B</math> 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,<ref>{{Cite journal|last=Starke|first=Gerhard|date=June 1993|title=Fejér-Walsh points for rational functions and their use in the ADI iterative method|journal=Journal of Computational and Applied Mathematics|volume=46|issue=1–2|pages=129–141|doi=10.1016/0377-0427(93)90291-i|issn=0377-0427|doi-access=free}}</ref> using the Ritz values of the matrices <math>A</math>, <math>A^{-1}</math>, <math>B</math>, and <math>B^{-1}</math> to formulate a greedy approach,<ref name=":7">{{Cite journal|last=Penzl|first=Thilo|date=January 1999|title=A Cyclic Low-Rank Smith Method for Large Sparse Lyapunov Equations|journal=SIAM Journal on Scientific Computing|language=en|volume=21|issue=4|pages=1401–1418|doi=10.1137/s1064827598347666|bibcode=1999SJSC...21.1401P |issn=1064-8275}}</ref> and cyclic methods, where the same small collection of shift parameters are reused until a convergence tolerance is met.<ref name=":7" /><ref name=":6" /> When the same shift parameter is used at every iteration, ADI is equivalent to an algorithm called Smith's method.<ref>{{Cite journal|last=Smith|first=R. A.|date=January 1968|title=Matrix Equation XA + BX = C|journal=SIAM Journal on Applied Mathematics|language=en|volume=16|issue=1|pages=198–201|doi=10.1137/0116017|issn=0036-1399}}</ref>

=== Factored ADI ===

Line 74 ⟶ 78:

and <math>\delta_p^2</math> is the central second difference operator for the ''p''-th coordinate

: <math>\delta_p^2 u_{ij}=u_{ij+e_p}-2u_{ij}+u_{ij-e_p}</math>

with <math>e_p=10(1,0)</math> or <math>01(0,1)</math> for <math>p=x</math> or <math>y</math> respectively (and <math>ij</math> a shorthand for lattice points <math>(i,j)</math>).

After performing a [[Von Neumann stability analysis|stability analysis]], it can be shown that this method will be stable for any <math>\Delta t</math>.

Line 90 ⟶ 94:

The system of equations involved is [[symmetric matrix|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.<ref>{{Citation | last1=Douglas, Jr. | first1=J. Jr.| title=On the numerical integration of u_xx+ u_yy= u_t by implicit methods | mr=0071875 | year=1955 | journal=Journal of the Society for Industrial and Applied Mathematics | volume=3 | pages=42–65}}.

</ref> There are more refined ADI methods such as the methods of Douglas,<ref>{{Citation | last1=Douglas Jr. | first1=Jim Jr.| title=Alternating direction methods for three space variables | doi=10.1007/BF01386295 | year=1962 | journal=Numerische Mathematik | issn=0029-599X | volume=4 | issue=1 | pages=41–63| s2cid=121455963 }}.</ref> or the f-factor method<ref>{{Citation | last1=Chang | first1=M. J. | last2=Chow | first2=L. C. | last3=Chang | first3=W. S. | title=Improved alternating-direction implicit method for solving transient three-dimensional heat diffusion problems | doi=10.1080/10407799108944957 | year=1991 | journal=Numerical Heat Transfer, Part B: Fundamentals | issn=1040-7790 | volume=19 | issue=1 | pages=69–84| bibcode=1991NHTB...19...69C | url=https://zenodo.org/record/1234469 }}.</ref> 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

: <math> \dot u = F_1 u + F_2 u, </math>

where <math> F_1 </math> and <math> F_2 </math> are (possibly nonlinear) operators defined on a [[Banach space|Banach]] space.<ref>{{cite book|last2=Verwer|first2=Jan|first1=Willem|last1=Hundsdorfer|title=Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations|date=2003|publisher=Springer Berlin Heidelberg|___location=Berlin, Heidelberg|isbn=978-3-662-09017-6}}</ref><ref>{{cite journal|last1=Lions|first1=P. L.|last2=Mercier|first2=B.|title=Splitting Algorithms for the Sum of Two Nonlinear Operators|journal=SIAM Journal on Numerical Analysis|date=December 1979|volume=16|issue=6|pages=964–979|doi=10.1137/0716071|bibcode=1979SJNA...16..964L}}</ref> In the diffusion example above we have <math> F_1 = {\partial^2 \over \partial x^2} </math> and <math> F_2 = {\partial^2 \over \partial y^2} </math>.

== Fundamental ADI (FADI) ==

=== Simplification of ADI to FADI ===

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,<ref>{{Cite journal|last=Tan|first=E. L.|date=2007|title=Efficient Algorithm for the Unconditionally Stable 3-D ADI-FDTD Method|url=https://www.ntu.edu.sg/home/eeltan/papers/2007%20Efficient%20Algorithm%20for%20the%20Unconditionally%20Stable%203-D%20ADI–FDTD%20Method.pdf|journal=IEEE Microwave and Wireless Components Letters|volume=17|issue=1|pages=7–9|doi=10.1109/LMWC.2006.887239|hdl=10356/138245 |s2cid=29025478}}</ref><ref name=":8">{{Cite journal|last=Tan|first=E. L.|date=2008|title=Fundamental Schemes for Efficient Unconditionally Stable Implicit Finite-Difference Time-Domain Methods|url=https://www.ntu.edu.sg/home/eeltan/papers/2008%20Fundamental%20Schemes%20for%20Efficient%20Unconditionally%20Stable%20Implicit%20Finite-Difference%20Time-Domain%20Methods.pdf|journal=IEEE Transactions on Antennas and Propagation|volume=56|issue=1|pages=170–177|doi=10.1109/TAP.2007.913089|arxiv=2011.14043|bibcode=2008ITAP...56..170T|hdl=10356/138249 |s2cid=37135325}}</ref> 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.

=== Relations to other implicit methods ===

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.<ref name=":8" /> 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 <ref>{{Cite journal|last=Tan|first=E. L.|date=2007|title=Unconditionally Stable LOD-FDTD Method for 3-D Maxwell's Equations|url=https://www.ntu.edu.sg/home/eeltan/papers/2007%20Unconditionally%20Stable%20LOD-FDTD%20Method%20for%203-D%20Maxwell’s%20Equations.pdf|journal=IEEE Microwave and Wireless Components Letters|volume=17|issue=2|pages=85–87|doi=10.1109/LMWC.2006.890166|hdl=10356/138296 |s2cid=22940993}}</ref><ref>{{Cite journal|last1=Gan|first1=T. H.|last2=Tan|first2=E. L.|date=2013|title=Unconditionally Stable Fundamental LOD-FDTD Method with Second-Order Temporal Accuracy and Complying Divergence|url=https://www.ntu.edu.sg/home/eeltan/papers/2013%20Unconditionally%20Stable%20Fundamental%20LOD-FDTD%20Method%20With%20Second-Order%20Temporal%20Accuracy%20and%20Complying%20Divergence.pdf|journal=IEEE Transactions on Antennas and Propagation|volume=61|issue=5|pages=2630–2638|doi=10.1109/TAP.2013.2242036|bibcode=2013ITAP...61.2630G|s2cid=7578037}}</ref> in [[computational 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.<ref>{{Cite journal|last1=Tay|first1=W. C.|last2=Tan|first2=E. L.|last3=Heh|first3=D. Y.|date=2014|title=Fundamental Locally One-Dimensional Method for 3-D Thermal Simulation|journal=IEICE Transactions on Electronics|volume=E-97-C|issue=7|pages=636–644|doi=10.1587/transele.E97.C.636|bibcode=2014IEITE..97..636T|hdl=10220/20410 |hdl-access=free}}</ref><ref>{{Cite journal|last1=Heh|first1=D. Y.|last2=Tan|first2=E. L.|last3=Tay|first3=W. C.|date=2016|title=Fast Alternating Direction Implicit Method for Efficient Transient Thermal Simulation of Integrated Circuits|journal=International Journal of Numerical Modelling: Electronic Networks, Devices and Fields|volume=29|issue=1|pages=93–108|doi=10.1002/jnm.2049|hdl=10356/137201 |s2cid=61039449 |hdl-access=free}}</ref>

==Further reading==

*{{cite journal|accessdate=March 28, 2025

|url=https://www.researchgate.net/publication/265623760 |via=ResearchGate.net

|title=50 Years of ADI Methods: Celebrating the Contributions of Jim Douglas, Don Peaceman, and Henry Rachford

|first1=Adam |last1=Usadi |first2=Clint |last2=Dawson

|journal=SIAM News |volume=39|number=2|date=March 2006}} (Provides a review of the ADI methods and variants over the years.)

== References ==

Line 106 ⟶ 125:

[[Category:Partial differential equations]]

[[Category:Numerical differential equations]]

[[Category:Iterative methods]]