Biconjugate gradient stabilized method: Difference between revisions

Revision as of 08:40, 12 March 2014

Latest revision as of 17:57, 29 July 2025

(41 intermediate revisions by 29 users not shown)

Line 1:

{{Short description|Concept in mathematics}}

In [[numerical linear algebra]], the '''biconjugate gradient stabilized method''', often abbreviated as '''BiCGSTAB''', is an [[iterative method]] developed by [[Henk van der Vorst|H. A. van der Vorst]] for the numerical solution of nonsymmetric [[System of linear equations|linear system]]s. It is a variant of the [[biconjugate gradient method]] (BiCG) and has faster and smoother convergence than the original BiCG as well as other variants such as the [[conjugate gradient squared method]] (CGS). It is a [[Krylov subspace]] method.

{{Technical|date=May 2015}}

In [[numerical linear algebra]], the '''biconjugate gradient stabilized method''', often abbreviated as '''BiCGSTAB''', is an [[iterative method]] developed by [[Henk van der Vorst|H. A. van der Vorst]] for the numerical solution of nonsymmetric [[System of linear equations|linear system]]s. It is a variant of the [[biconjugate gradient method]] (BiCG) and has faster and smoother convergence than the original BiCG as well as other variants such as the [[conjugate gradient squared method]] (CGS). It is a [[Krylov subspace]] method. Unlike the original BiCG method, it doesn't require multiplication by the transpose of the system matrix.

==Algorithmic steps==

===Unpreconditioned BiCGSTAB===

In the following sections, {{math|1=('''<var>x</var>''','''<var>y</var>''') = '''<var>x</var>'''^T '''<var>y</var>'''}} denotes the [[dot product]] of vectors. To solve a linear system {{math|'''<var>Ax</var>''' {{=}} '''<var>b</var>'''}}, BiCGSTAB starts with an initial guess {{math|'''<var>x</var>'''₀}} and proceeds as follows:

# {{math|'''<var>r</var>'''₀ {{=}} '''<var>b</var>''' − '''<var>Ax</var>'''₀}}

# Choose an arbitrary vector {{math|'''<var>r̂</var>'''₀}} such that {{math|('''<var>r̂</var>'''₀, '''<var>r</var>'''₀) ≠ 0}}, e.g., {{math|'''<var>r̂</var>'''₀ {{=}} '''<var>r</var>'''₀}}

# {{math|<var>ρ</var>₀ {{=}} ('''<var>αr̂</var>'''₀, {{=}} '''<var>ωr</var>'''₀) {{=}} 1}}

# {{math|'''<var>vp</var>'''₀ {{=}} '''<var>pr</var>'''₀ {{=}} '''0'''}}

# For {{math|<var>i</var> {{=}} 1, 2, 3, …}}

## {{math|'''<var>ρ_iv</var>''' {{=}} ('''<var>r̂Ap</var>'''₀, '''<var>ri</var>'''_{<var>i<var>−1})}}

## {{math|<var>βα</var> {{=}} (<var>ρ_i</var>/<var>ρ</var>_{<var>i</var>−1})/('''<var>αr̂</var>/'''<varsub>ω0</var>_{, <var>i'''v'''</var>−1})}}

## {{math|<var>'''ph'''_i</var> {{=}} '''<var>rx</var>'''_{<var>i</var>−1} + <var>β</var>(α'''<var>p'''</var>'''_{<var>i</var>−1} −}}
## {{math|<var>ω'''s'''</var>_{{{=}} <var>i<var>−1}'''<var>vr'''</var>'''_{<var>i</var>−1}) − <var>α'''v'''</var>}}

## If {{math|<var>'''vh'''_i</var> {{=}} is accurate enough, i.e., if <var>'''s'''</var>Ap is small enough, then set {{math|</var>'''x'''<sub>i<var/sub>i</var> {{=}} <var>'''h'''</subvar>}} and quit

## {{math|<var>α</var> {{=}} <var>ρ_i</var>/('''<var>r̂</var>'''₀, <var>'''v'''_i</var>)}}

## {{math|'''<var>s</var>''' {{=}} '''<var>r</var>'''_{<var>i<var>−1} − <var>α'''v'''_i</var>}}

## {{math|'''<var>t</var>''' {{=}} '''<var>As</var>'''}}

## {{math|<var>ω_i</var> {{=}} (<var>'''t'''</var>t, </var>''', s'''</var>s<)/(<var>''')/(t'''</var>t, </var>''', t'''<var>t</var>''')}}

## {{math|<var>'''x'''_i</var> {{=}} '''<var>x</var>'''_{<var>i<var>−1} + <var>α'''ph'''_i</var> + <var>ω_i'''s'''</var>}}

## If {{math|<var>'''xr'''_i</var> {{=}} is<var>'''s'''</var> accurate− enough then quit<var>ω'''t'''</var>}}

## If {{math|<var>'''rx'''_i</var>}} is accurate enough, i.e., if {{=}} '''math|<var>s</var>''' − <var>ωr'''_i'''t'''</var>}} is small enough, then quit

## {{math|<var>ρ_i</var> {{=}} ('''<var>r̂</var>'''₀, '''<var>r</var>'''_<var>i</var>)}}

## {{math|<var>αβ</var> {{=}} (<var>ρ_i</var>/('''<var>r̂ρ</var>'''_{0</subvar>, i</var>'''v'''−1}i)(<var>α</subvar>/<var>ω</var>)}}

## {{math|<var>'''p'''_i</var> {{=}} '''<var>r</var>'''_<var>i</var> + <var>β</var>('''<var>p</var>'''_{<var>i</var>−1} − <var>ω</var>'''<var>v</var>''')}}

In some cases, choosing the vector {{math|'''<var>r̂</var>'''₀}} randomly improves numerical stability.<ref>{{Cite journal |last=Schoutrop |first=Chris |last2=Boonkkamp |first2=Jan ten Thije |last3=Dijk |first3=Jan van |date=July 2022 |title=Reliability Investigation of BiCGStab and IDR Solvers for the Advection-Diffusion-Reaction Equation |url=https://doi.org/10.4208/cicp.OA-2021-0182 |journal=Communications in Computational Physics |language=en |volume=32 |issue=1 |pages=156–188 |doi=10.4208/cicp.oa-2021-0182 |issn=1815-2406}}</ref>

===Preconditioned BiCGSTAB===

Line 27 ⟶ 33:

# {{math|'''<var>r</var>'''₀ {{=}} '''<var>b</var>''' − '''<var>Ax</var>'''₀}}

# Choose an arbitrary vector {{math|'''<var>r̂</var>'''₀}} such that {{math|('''<var>r̂</var>'''₀, '''<var>r</var>'''₀) ≠ 0}}, e.g., {{math|'''<var>r̂</var>'''₀ {{=}} '''<var>r</var>'''₀}}

# {{math|<var>ρ</var>₀ {{=}} ('''<var>αr̂</var>'''₀, {{=}} '''<var>ωr</var>'''₀) {{=}} 1}}

# {{math|'''<var>vp</var>'''₀ {{=}} '''<var>pr</var>'''₀ {{=}} '''0'''}}

# For {{math|<var>i</var> {{=}} 1, 2, 3, …}}

## {{math|'''<var>ρ_iy</var>''' {{=}} ({{SubSup|'''<var>r̂K</var>'''|2|−1}}{{SubSup|'''<subvar>0K</subvar>, '''|1|−1}}'''<var>rp</var>'''_{<var>i</var>−1})}}

## {{math|<var>β'''v'''</var> {{=}} ('''<var>ρ_iAy</var>/<var>ρ</var>_{<var>i<var>−1})(<var>α</var>/<var>ω</var>_{<var>i<var>−1})'''}}

## {{math|<var>'''p'''_iα</var> {{=}} '''<var>rρ</var>'''_{<var>i</var>−1} + <var>β</var>('''<var>pr̂</var>'''_{<var>i<var>−10} −, <var>ω</var>_{<var>i<var>−1}'''<var>v</var>'''_{</var>i<var>−1})}}

## {{math|<var>'''h'''<var>y</var>''' {{=}} '''<var>Kx</var>'''<supsub>−1</sup>'''<var>pi</var>'''−1</sub> + <var>iα'''y'''</var></sub> }}

## {{math|<var>'''v'''<sub>i</subvar>s</var>''' {{=}} '''<var>Ayr</var>'''_{<var>i</var>−1} − <var>α'''v'''</var>}}

## If {{math|<var>α'''h'''</var>}} is accurate enough then {{=}} math|<var>ρ_i</var>/('''<var>r̂</var>x'''_0i,</var> {{=}} <var>'''vh'''_i</var>)}} and quit

## {{math|'''<var>sz</var>''' {{=}} {{SubSup|'''<var>rK</var>'''_{|2|−1}}{{SubSup|'''<var>i<var>−1K} − <var>α'''v|1|−1}}'''<sub>i</subvar>s</var>'''}}

## {{math|'''<var>z</var>''' {{=}} '''<var>K</var>'''⁻¹'''<var>s</var>'''}}

## {{math|'''<var>t</var>''' {{=}} '''<var>Az</var>'''}}

## {{math|<var>ω_i</var> {{=}} ({{SubSup|'''<var>K</var>'''|1|−1}}'''<var>t</var>''', {{SubSup|'''<var>K</var>'''|1|−1}}'''<var>s</var>''')/({{SubSup|'''<var>K</var>'''|1|−1}}'''<var>t</var>''', {{SubSup|'''<var>K</var>'''|1|−1}}'''<var>t</var>''')}}

## {{math|<var>'''x'''_i</var> {{=}} '''<var>x</var>'''_{<var>i<var>−1} + <var>α'''yh'''</var> + <var>ω_i'''z'''</var>}}

## {{math|<var>'''r'''_i</var> {{=}} '''<var>s</var>''' − <var>ω'''t'''</var>}}

## If {{math|<var>'''x'''_i</var>}} is accurate enough then quit

## {{math|<var>'''r'''ρ_i</var> {{=}} ('''<var>sr̂</var>''' − <var>ω_i0, '''t<var>r</var>'''_<var>i</var>)}}

## {{math|'''<var>sβ</var>''' {{=}} '''(<var>rρ_i</var>/<var>ρ</var>'''_{<var>i</var>−1} − )(<var>α'''v'''<sub/var>i</sub<var>ω</var>)}}

## {{math|<var>'''p'''_i</var> {{=}} '''<var>r</var>'''_<var>i</var> + <var>β</var>('''<var>p</var>'''_{<var>i</var>−1} − <var>ω</var>'''<var>v</var>''')}}

This formulation is equivalent to applying unpreconditioned BiCGSTAB to the explicitly preconditioned system

Line 77 ⟶ 85:

:{{math|<var>P_i</var>('''<var>A</var>''') {{=}} <var>P</var>_{<var>i</var>−1}('''<var>A</var>''') − <var>α_i'''A'''T</var>_{<var>i</var>−1}('''<var>A</var>''')}},

:{{math|<var>T_i</var>('''<var>A</var>''') {{=}} <var>P_i</var>('''<var>A</var>''') −+ <var>β</var>_{<var>i</var>+1}<var>T</var>_{<var>i</var>−1}('''<var>A</var>''')}}.

===Derivation of BiCGSTAB from BiCG===

Line 84 ⟶ 92:

:{{math|'''<var>r̃</var>'''_<var>i</var> {{=}} <var>Q_i</var>('''<var>A</var>''')<var>P_i</var>('''<var>A</var>''')'''<var>r</var>'''₀}}

where {{math|<var>Q_i</var>('''<var>A</var>''') {{=}} ('''<var>I</var>''' − <var>ω</var>₁'''<var>A</var>''')('''<var>I</var>''' − <var>ω</var>₂'''<var>A</var>''')⋯('''<var>I</var>''' − <var>ω_i'''A'''</var>)}} with suitable constants {{math|<var>ω_j</var>}} instead of {{math|<var>'''r'''_i</var> {{=}} <var>P_i</var>('''<var>A</var>''')<var>'''r'''₀</var>}} in the hope that {{math|<var>Q_i</var>('''<var>A</var>''')}} will enable faster and smoother convergence in {{math|<var>'''r̃'''_i</var>}} than {{math|<var>'''r'''_i</var>}}.

It follows from the recurrence relations for {{math|<var>P_i</var>('''<var>A</var>''')}} and {{math|<var>T_i</var>('''<var>A</var>''')}} and the definition of {{math|<var>Q_i</var>('''<var>A</var>''')}} that

Line 92 ⟶ 100:

which entails the necessity of a recurrence relation for {{math|<var>Q_i</var>('''<var>A</var>''')<var>T_i</var>('''<var>A</var>''')'''<var>r</var>'''₀}}. This can also be derived from the BiCG relations:

:{{math|<var>Q_i</var>('''<var>A</var>''')<var>T_i</var>('''<var>A</var>''')'''<var>r</var>'''₀ {{=}} <var>Q_i</var>('''<var>A</var>''')<var>P_i</var>('''<var>A</var>''')'''<var>r</var>'''₀ + <var>β</var>_{<var>i</var>+1}('''<var>I</var>''' − <var>ω_i'''A'''</var>)<var>Q</var>_{<var>i</var>−1}('''<var>A</var>''')<var>PT</var>_{<var>i</var>−1}('''<var>A</var>''')'''<var>r</var>'''₀}}.

Similarly to defining {{math|<var>'''r̃'''_i</var>}}, BiCGSTAB defines

Line 126 ⟶ 134:

:{{math|<var>ρ̃</var>_<var>i</var> {{=}} (<var>Q</var>_{<var>i</var>−1}('''<var>A</var>'''^T)'''<var>r̂</var>'''₀, <var>P</var>_{<var>i</var>−1}('''<var>A</var>''')'''<var>r</var>'''₀) {{=}} ('''<var>r̂</var>'''₀, <var>Q</var>_{<var>i</var>−1}('''<var>A</var>''')<var>P</var>_{<var>i</var>−1}('''<var>A</var>''')'''<var>r</var>'''₀) {{=}} ('''<var>r̂</var>'''₀, '''<var>r</var>'''_{<var>i</var>−1})}}.

Due to biorthogonality, {{math|'''<var>r</var>'''_{<var>i</var>−1} {{=}} <var>P</var>_{<var>i</var>−1}('''<var>A</var>''')'''<var>r</var>'''₀}} is orthogonal to {{math|<var>U</var>_{<var>i</var>−2}('''<var>A</var>'''^T)'''<var>r̂</var>'''₀}} where {{math|<var>U</var>_{<var>i</var>−2}('''<var>A</var>'''^T)}} is any polynomial of degree {{math|<var>i</var> − 2}} in {{math|'''<var>A</var>'''^T}}. Hence, only the highest-order terms of {{math|<var>P</var>_{<var>i</var>−1}('''<var>A</var>'''^T)}} and {{math|<var>Q</var>_{<var>i</var>−1}('''<var>A</var>'''^T)}} matter in the [[dot productsproduct]]s {{math|(<var>P</var>_{<var>i</var>−1}('''<var>A</var>'''^T)'''<var>r̂</var>'''₀, <var>P</var>_{<var>i</var>−1}('''<var>A</var>''')'''<var>r</var>'''₀)}} and {{math| (<var>Q</var>_{<var>i</var>−1}('''<var>A</var>'''^T)'''<var>r̂</var>'''₀, <var>P</var>_{<var>i</var>−1}('''<var>A</var>''')'''<var>r</var>'''₀)}}. The leading coefficients of {{math|<var>P</var>_{<var>i</var>−1}('''<var>A</var>'''^T)}} and {{math|<var>Q</var>_{<var>i</var>−1}('''<var>A</var>'''^T)}} are {{math|(−1)^{<var>i</var>−1}<var>α</var>₁<var>α</var>₂⋯<var>α</var>_{<var>i</var>−1}}} and {{math|(−1)^{<var>i</var>−1}<var>ω</var>₁<var>ω</var>₂⋯<var>ω</var>_{<var>i</var>−1}}}, respectively. It follows that

:{{math|<var>ρ_i</var> {{=}} (<var>α</var>₁/<var>ω</var>₁)(<var>α</var>₂/<var>ω</var>₂)⋯(<var>α</var>_{<var>i</var>−1}/<var>ω</var>_{<var>i</var>−1})<var>ρ̃</var>_<var>i</var>}},

Line 138 ⟶ 146:

:{{math|<var>α_i</var> {{=}} <var>ρ_i</var>/('''<var>p̂</var>'''_<var>i</var>, <var>'''Ap'''_i</var>) {{=}} (<var>P</var>_{<var>i</var>−1}('''<var>A</var>'''^T)'''<var>r̂</var>'''₀, <var>P</var>_{<var>i</var>−1}('''<var>A</var>''')'''<var>r</var>'''₀)/(<var>T</var>_{<var>i</var>−1}('''<var>A</var>'''^T)'''<var>r̂</var>'''₀, <var>'''A'''T</var>_{<var>i</var>−1}('''<var>A</var>''')'''<var>r</var>'''₀)}}.

Similarly to the case above, only the highest-order terms of {{math|<var>P</var>_{<var>i</var>−1}('''<var>A</var>'''^T)}} and {{math|<var>T</var>_{<var>i</var>−1}('''<var>A</var>'''^T)}} matter in the [[dot productsproduct]]s thanks to biorthogonality and biconjugacy. It happens that {{math|<var>P</var>_{<var>i</var>−1}('''<var>A</var>'''^T)}} and {{math|<var>T</var>_{<var>i</var>−1}('''<var>A</var>'''^T)}} have the same leading coefficient. Thus, they can be replaced simultaneously with {{math|<var>Q</var>_{<var>i</var>−1}('''<var>A</var>'''^T)}} in the formula, which leads to

:{{math|<var>α_i</var> {{=}} (<var>Q</var>_{<var>i</var>−1}('''<var>A</var>'''^T)'''<var>r̂</var>'''₀, <var>P</var>_{<var>i</var>−1}('''<var>A</var>''')'''<var>r</var>'''₀)/(<var>Q</var>_{<var>i</var>−1}('''<var>A</var>'''^T)'''<var>r̂</var>'''₀, <var>'''A'''T</var>_{<var>i</var>−1}('''<var>A</var>''')'''<var>r</var>'''₀) {{=}} <var>ρ̃</var>_<var>i</var>/('''<var>r̂</var>'''₀, '''<var>A</var>'''<var>Q</var>_{<var>i</var>−1}('''<var>A</var>''')<var>T</var>_{<var>i</var>−1}('''<var>A</var>''')'''<var>r</var>'''₀) {{=}} <var>ρ̃</var>_<var>i</var>/('''<var>r̂</var>'''₀, '''<var>Ap̃</var>'''_<var>i</var>)}}.

Line 151 ⟶ 159:

==Generalization==

BiCGSTAB can be viewed as a combination of BiCG and [[GMRES]] where each BiCG step is followed by a GMRES({{math|1}}) (i.e., GMRES restarted at each step) step to repair the irregular convergence behavior of CGS, as an improvement of which BiCGSTAB was developed. However, due to the use of degree-one minimum residual polynomials, such repair may not be effective if the matrix {{math|'''<var>A</var>'''}} has large complex eigenpairs. In such cases, BiCGSTAB is likely to stagnate, as confirmed by numerical experiments.

One may expect that higher-degree minimum residual polynomials may better handle this situation. This gives rise to algorithms including BiCGSTAB2{{ref|bicgstab2}} and the more general BiCGSTAB({{math|<var>l</var>}}){{ref|bicgstab(l)}}. In BiCGSTAB({{math|<var>l</var>}}), a GMRES({{math|<var>l</var>}}) step follows every {{math|<var>l</var>}} BiCG steps. BiCGSTAB2 is equivalent to BiCGSTAB({{math|<var>l</var>}}) with {{math|<var>l</var> {{=}} 2}}.

Line 161 ⟶ 169:

==References==

{{reflist}}

* {{cite doi|10.1137/0913035}}

* {{Cite journal | doi = 10.1137/0913035 | title = Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems | year = 1992 | last1 = Van der Vorst | first1 = H. A. | journal = [[SIAM Journal on Scientific Computing|SIAM J. Sci. Stat. Comput.]] | volume = 13 | issue = 2 | pages = 631–644 | hdl = 10338.dmlcz/104566 | hdl-access = free }}

* {{cite book

|last = Saad

|first = Y.

|title = Iterative Methods for Sparse Linear Systems

|chapter-url = https://archive.org/details/iterativemethods0000saad

|chapter-url-access = registration

|edition = 2nd

|chapter = §7.4.2 BICGSTAB

|pages = [https://archive.org/details/iterativemethods0000saad/page/231 231–234]

|year = 2003

|publisher = SIAM

|isbn = 978-0-89871-534-7

|doi = 10.2277/0898715342

}}

* {{note|bicgstab2}}{{Cite journal| doi = 10.1137/0914062| title = Variants of BICGSTAB for Matrices with Complex Spectrum| year = 1993| last1 = Gutknecht | first1 = M. H.| journal = [[SIAM Journal on Scientific Computing|SIAM J. Sci. Comput.]]| volume = 14| issue = 5| pages = 1020–1033 }}

* {{note|bicgstab2}}{{cite doi|10.1137/0914062}}

* {{note|bicgstab(l)}}{{cite journal

| last1 = Sleijpen
| first1 = G. L. G.

| last2 = Fokkema
| first2 = D. R.

| date = November 1993

| title = BiCGstab(''l'') for linear equations involving unsymmetric matrices with complex spectrum

| journal = [[Electronic Transactions on Numerical Analysis]]

| volume = 1

| pages = 11–32

Line 186 ⟶ 198:

| ___location = Kent, OH

| issn = 1068-9613

| url = http://www.emis.ams.orgde/journals/ETNA/vol.1.1993/pp11-32.dir/pp11-32.pdf

| format = PDF

}}