Biconjugate gradient method

In mathematics, more specifically in numerical linear algebra, the biconjugate gradient method is an algorithm to solve systems of linear equations

A x = b .

Unlike the conjugate gradient method, this algorithm does not require the matrix $A$ to be self-adjoint, but instead one needs to perform multiplications by the conjugate transpose Template:Math.

The Algorithm

Choose initial guess $x_{0}$ , two other vectors $x_{0}^{*}$ and $b^{*}$ and a preconditioner $M$
$r_{0} \leftarrow b - A x_{0}$
$r_{0}^{*} \leftarrow b^{*} - x_{0}^{*} A^{*}$
$p_{0} \leftarrow M^{- 1} r_{0}$
$p_{0}^{*} \leftarrow r_{0}^{*} M^{- 1}$
for k=0,1,… do
1. $α_{k} \leftarrow \frac{r_{k}^{*} M^{- 1} r_{k}}{p_{k}^{*} A p_{k}}$
2. $x_{k + 1} \leftarrow x_{k} + α_{k} \cdot p_{k}$
3. $x_{k + 1}^{*} \leftarrow x_{k}^{*} + \overline{α_{k}} \cdot p_{k}^{*}$
4. $r_{k + 1} \leftarrow r_{k} - α_{k} \cdot A p_{k}$
5. $r_{k + 1}^{*} \leftarrow r_{k}^{*} - \overline{α_{k}} \cdot p_{k}^{*} A^{*}$
6. $β_{k} \leftarrow \frac{r_{k + 1}^{*} M^{- 1} r_{k + 1}}{r_{k}^{*} M^{- 1} r_{k}}$
7. $p_{k + 1} \leftarrow M^{- 1} r_{k + 1} + β_{k} \cdot p_{k}$
8. $p_{k + 1}^{*} \leftarrow r_{k + 1}^{*} M^{- 1} + \overline{β_{k}} \cdot p_{k}^{*}$

In the above formulation, the computed $r_{k}$ and $r_{k}^{*}$ satisfy

r_{k} = b - A x_{k},

r_{k}^{*} = b^{*} - x_{k}^{*} A^{*}

and thus are the respective residuals corresponding to $x_{k}$ and $x_{k}^{*}$ , as approximate solutions to the systems

A x = b,

x^{*} A^{*} = b^{*};

$x^{*}$ is the adjoint, and $\overline{α}$ is the complex conjugate.

The biconjugate gradient method is numerically unstable Template:Citation needed (compare to the biconjugate gradient stabilized method), but very important from a theoretical point of view. Define the iteration steps by

x_{k} := x_{j} + P_{k} A^{- 1} (b - A x_{j}),

x_{k}^{*} := x_{j}^{*} + (b^{*} - x_{j}^{*} A) P_{k} A^{- 1},

where $j < k$ using the related projection

P_{k} := 𝐮_{k} {(𝐯_{k}^{*} A 𝐮_{k})}^{- 1} 𝐯_{k}^{*} A,

with

𝐮_{k} = [u_{0}, u_{1}, \dots, u_{k - 1}],

𝐯_{k} = [v_{0}, v_{1}, \dots, v_{k - 1}] .

These related projections may be iterated themselves as

P_{k + 1} = P_{k} + (1 - P_{k}) u_{k} \otimes \frac{v_{k}^{*} A (1 - P_{k})}{v_{k}^{*} A (1 - P_{k}) u_{k}} .

A relation to Quasi-Newton methods is given by $P_{k} = A_{k}^{- 1} A$ and $x_{k + 1} = x_{k} - A_{k + 1}^{- 1} (A x_{k} - b)$ , where

A_{k + 1}^{- 1} = A_{k}^{- 1} + (1 - A_{k}^{- 1} A) u_{k} \otimes \frac{v_{k}^{*} (1 - A A_{k}^{- 1})}{v_{k}^{*} A (1 - A_{k}^{- 1} A) u_{k}} .

The new directions

p_{k} = (1 - P_{k}) u_{k},

p_{k}^{*} = v_{k}^{*} A (1 - P_{k}) A^{- 1}

are then orthogonal to the residuals:

v_{i}^{*} r_{k} = p_{i}^{*} r_{k} = 0,

r_{k}^{*} u_{j} = r_{k}^{*} p_{j} = 0,

which themselves satisfy

r_{k} = A (1 - P_{k}) A^{- 1} r_{j},

r_{k}^{*} = r_{j}^{*} (1 - P_{k})

where $i, j < k$ .

The biconjugate gradient method now makes a special choice and uses the setting

u_{k} = M^{- 1} r_{k},

v_{k}^{*} = r_{k}^{*} M^{- 1} .

With this particular choice, explicit evaluations of $P_{k}$ and Template:Math are avoided, and the algorithm takes the form stated above.

If $A = A^{*}$ is self-adjoint, $x_{0}^{*} = x_{0}$ and $b^{*} = b$ , then $r_{k} = r_{k}^{*}$ , $p_{k} = p_{k}^{*}$ , and the conjugate gradient method produces the same sequence $x_{k} = x_{k}^{*}$ at half the computational cost.
The sequences produced by the algorithm are biorthogonal, i.e., $p_{i}^{*} A p_{j} = r_{i}^{*} M^{- 1} r_{j} = 0$ for $i \neq j$ .
if $P_{j^{'}}$ is a polynomial with $\deg (P_{j^{'}}) + j < k$ , then $r_{k}^{*} P_{j^{'}} (M^{- 1} A) u_{j} = 0$ . The algorithm thus produces projections onto the Krylov subspace.
if $P_{i^{'}}$ is a polynomial with $i + \deg (P_{i^{'}}) < k$ , then $v_{i}^{*} P_{i^{'}} (A M^{- 1}) r_{k} = 0$ .