Biconjugate gradient stabilized method: Difference between revisions

Revision as of 07:11, 13 May 2011

Revision as of 07:53, 13 May 2011

Line 14:

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

## {{math|<var>'''v'''_i</var> {{=}} '''<var>Ap</var>'''_<var>i</var>}}

## {{math|<var>α</var> {{=}} <var>ρ_i</var>/('''<var>rr̂</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>'''}}

Line 50:

==Derivation==

===BiCG in polynomial form===

In BiCG, the search directions {{math|<var>'''p'''_i</var>}} and {{math|'''<var>pp̂</var>̂'''_<var>i</var>}} and the residuals {{math|<var>'''r'''_i</var>}} and {{math|'''<var>rr̂</var>̂'''_<var>i</var>}} are updated using the following recurrence relations:

:{{math|<var>'''p'''_i</var> {{=}} '''<var>r</var>'''_{<var>i</var>−1} + <var>β_i</var>'''<var>p</var>'''_{<var>i</var>−1}}},

:{{math|'''<var>pp̂</var>̂'''_<var>i</var> {{=}} '''<var>rr̂</var>̂'''_{<var>i</var>−1} + <var>β_i</var>'''<var>pp̂</var>̂'''_{<var>i</var>−1}}},

:{{math|<var>'''r'''_i</var> {{=}} '''<var>r</var>'''_{<var>i</var>−1} − <var>α_i</var>'''<var>Ap</var>'''_<var>i</var>}},

:{{math|'''<var>rr̂</var>̂'''_<var>i</var> {{=}} '''<var>rr̂</var>̂'''_{<var>i</var>−1} − <var>α_i</var>'''<var>A</var>'''^T'''<var>pp̂</var>̂'''_<var>i</var>}}.

The constants {{math|<var>α_i</var>}} and {{math|<var>β_i</var>}} are chosen to be

:{{math|<var>α_i</var> {{=}} <var>ρ_i</var>/('''<var>pp̂</var>̂'''_<var>i</var>, <var>'''Ap'''_i</var>)}},

:{{math|<var>β_i</var> {{=}} <var>ρ_i</var>/<var>ρ</var>_{<var>i</var>−1}}}

where {{math|<var>ρ_i</var> {{=}} ('''<var>rr̂</var>̂'''_{<var>i</var>−1}, '''<var>r</var>'''_{<var>i</var>−1})}} so that the residuals and the search directions satisfy biorthogonality and biconjugacy, respectively, i.e., for {{math|<var>i</var> ≠ <var>j</var>}},

:{{math|('''<var>rr̂</var>̂'''_<var>i</var>, <var>'''r'''_j</var>) {{=}} 0}},

:{{math|('''<var>pp̂</var>̂'''_<var>i</var>, <var>'''Ap'''_j</var>) {{=}} 0}}.

It is straightforward to show that

:{{math|<var>'''r'''_i</var> {{=}} <var>P_i</var>('''<var>A</var>''')'''<var>r</var>'''₀}},

:{{math|'''<var>rr̂</var>̂'''_<var>i</var> {{=}} <var>P_i</var>('''<var>A</var>'''^T)'''<var>rr̂</var>̂'''₀}},

:{{math|'''<var>p</var>'''_{<var>i</var>+1} {{=}} <var>T_i</var>('''<var>A</var>''')'''<var>r</var>'''₀}},

:{{math|'''<var>pp̂</var>̂'''_{<var>i</var>+1} {{=}} <var>T_i</var>('''<var>A</var>'''^T)'''<var>rr̂</var>̂'''₀}}

where {{math|<var>P_i</var>('''<var>A</var>''')}} and {{math|<var>T_i</var>('''<var>A</var>''')}} are {{math|<var>i</var>}}th-degree polynomials in {{math|'''<var>A</var>'''}}. These polynomials satisfy the following recurrence relations:

Line 82:

It is unnecessary to explicitly keep track of the residuals and search directions of BiCG. In other words, the BiCG iterations can be performed implicitly. In BiCGSTAB, one wishes to have recurrence relations for

:{{math|'''<var>rr̃</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>''')}} 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>}}.

Line 96:

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

:{{math|'''<var>pp̃</var>̃'''_{<var>i</var>+1} {{=}} <var>Q_i</var>('''<var>A</var>''')<var>T_i</var>('''<var>A</var>''')'''<var>r</var>'''₀}}.

Written in vector form, the recurrence relations for {{math|'''<var>pp̃</var>̃'''_<var>i</var>}} and {{math|'''<var>rr̃</var>̃'''_<var>i</var>}} are

:{{math|'''<var>pp̃</var>̃'''_<var>i</var> {{=}} '''<var>rr̃</var>̃'''_{<var>i</var>−1} + <var>β_i</var>('''<var>I</var>''' − <var>ω</var>_{<var>i</var>−1}'''<var>A</var>''')'''<var>pp̃</var>̃'''_{<var>i</var>−1}}},

:{{math|'''<var>rr̃</var>̃'''_<var>i</var> {{=}} ('''<var>I</var>''' − <var>ω_i'''A'''</var>)('''<var>rr̃</var>̃'''_{<var>i</var>−1} − <var>α_i'''A'''</var>'''<var>pp̃</var>̃'''_<var>i</var>)}}.

To derive a recurrence relation for {{math|<var>'''x'''_i</var>}}, define

:{{math|<var>'''s'''_i</var> {{=}} '''<var>rr̃</var>̃'''_{<var>i</var>−1} − <var>α_i'''A'''</var>'''<var>pp̃</var>̃'''_<var>i</var>}}.

The recurrence relation for {{math|'''<var>rr̃</var>̃'''_<var>i</var>}} can then be written as

:{{math|'''<var>rr̃</var>̃'''_<var>i</var> {{=}} '''<var>rr̃</var>̃'''_{<var>i</var>−1} − <var>α_i'''A'''</var>'''<var>pp̃</var>̃'''_<var>i</var> − <var>ω_i'''As'''_i</var>}},

which corresponds to

:{{math|'''<var>x</var>'''_<var>i</var> {{=}} '''<var>x</var>'''_{<var>i</var>−1} + <var>α_i</var>'''<var>pp̃</var>̃'''_<var>i</var> + <var>ω_i'''s'''_i</var>}}.

===Determination of BiCGSTAB constants===

Line 120:

In BiCG, {{math|<var>β_i</var> {{=}} <var>ρ_i</var>/<var>ρ</var>_{<var>i</var>−1}}} with

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

Since BiCGSTAB does not explicitly keep track of {{math|'''<var>rr̂</var>̂'''_<var>i</var>}} or {{math|'''<var>r</var>'''_<var>i</var>}}, {{math|<var>ρ_i</var>}} is not immediately computable from this formula. However, it can be related to the scalar

:{{math|<var>ρρ̃</var>̃_<var>i</var> {{=}} (<var>Q</var>_{<var>i</var>−1}('''<var>A</var>'''^T)'''<var>rr̂</var>̂'''₀, <var>P</var>_{<var>i</var>−1}('''<var>A</var>''')'''<var>r</var>'''₀) {{=}} ('''<var>rr̂</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>rr̂</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>rr̂</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 products {{math|(<var>P</var>_{<var>i</var>−1}('''<var>A</var>'''^T)'''<var>rr̂</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>rr̂</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>}},

and thus

:{{math||<var>β_i</var> {{=}} <var>ρ_i</var>/<var>ρ</var>_{<var>i</var>−1} {{=}} (<var>ρρ̃</var>̃_<var>i</var>/<var>ρρ̃</var>̃_{<var>i</var>−1})(<var>α</var>_{<var>i</var>−1}/<var>ω</var>_{<var>i</var>−1})}}.

A simple formula for {{math|<var>α_i</var>}} can be similarly derived. In BiCG,

:{{math|<var>α_i</var> {{=}} <var>ρ_i</var>/('''<var>pp̂</var>̂'''_<var>i</var>, <var>'''Ap'''_i</var>) {{=}} (<var>P</var>_{<var>i</var>−1}('''<var>A</var>'''^T)'''<var>rr̂</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>rr̂</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 products 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>rr̂</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>rr̂</var>̂'''₀, <var>'''A'''T</var>_{<var>i</var>−1}('''<var>A</var>''')'''<var>r</var>'''₀) {{=}} <var>ρρ̃</var>̃_<var>i</var>/('''<var>rr̂</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>rr̂</var>̂'''₀, '''<var>ApAp̃</var>̃'''_<var>i</var>)}}.

Finally, BiCGSTAB selects {{math|<var>ω_i</var>}} to minimize {{math|'''<var>rr̃</var>̃'''_<var>i</var> {{=}} ('''I''' − <var>ω_i'''A'''</var>)<var>'''s'''_i</var>}} in {{math|2}}-norm as a function of {{math|<var>ω_i</var>}}. This is achieved when

:{{math|(('''I''' − <var>ω_i'''A'''</var>)<var>'''s'''_i</var>, <var>'''As'''_i</var>) {{=}} 0}},