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 .

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^ r_0\,$
# $p_0^* \leftarrow r_0^*M^\,$
# for $k=0, 1, \ldots$ do
## $\alpha_k \leftarrow \,$
## $x_ \leftarrow x_k + \alpha_k \cdot p_k\,$
## $x_^* \leftarrow x_k^* + \overline\cdot p_k^*\,$
## $r_ \leftarrow r_k - \alpha_k \cdot A p_k\,$
## $r_^* \leftarrow r_k^*- \overline \cdot p_k^*\, A$
## $\beta_k \leftarrow \,$
## $p_ \leftarrow M^ r_ + \beta_k \cdot p_k\,$
## $p_^* \leftarrow r_^*M^ + \overline\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.

Unpreconditioned version of the algorithm

# Choose initial guess $x_0\,$,
# $r_0 \leftarrow b-A\, x_0\,$
# $\hat_0 \leftarrow \hat - \hat_0A$
# $p_0 \leftarrow r_0\,$
# $\hat_0 \leftarrow \hat_0\,$
# for $k=0, 1, \ldots$ do
## $\alpha_k \leftarrow \,$
## $x_ \leftarrow x_k + \alpha_k \cdot p_k\,$
## $\hat_ \leftarrow \hat_k + \alpha_k \cdot \hat_k\,$
## $r_ \leftarrow r_k - \alpha_k \cdot A p_k\,$
## $\hat_ \leftarrow \hat_k- \alpha_k \cdot \hat_k A$
## $\beta_k \leftarrow \,$
## $p_ \leftarrow r_ + \beta_k \cdot p_k\,$
## $\hat_ \leftarrow \hat_ + \beta_k \cdot \hat_k\,$

Discussion

The biconjugate gradient method is numerically unstable (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^\left(b - A x_j \right),$
:$x_k^*:= x_j^*+\left(b^*- x_j^* A \right) P_k A^,$

where j using the related projection

:$P_k:= \mathbf_k \left(\mathbf_k^* A \mathbf_k \right)^ \mathbf_k^* A,$

with

:\mathbf_k=\left _0, u_1, \dots, u_ \right
:\mathbf_k=\left _0, v_1, \dots, v_ \right

These related projections may be iterated themselves as

:$P_= P_k+ \left( 1-P_k\right) u_k \otimes .$

A relation to Quasi-Newton methods is given by $P_k= A_k^ A$ and $x_= x_k- A_^\left(A x_k -b \right)$, where
:$A_^= A_k^+ \left( 1-A_k^A\right) u_k \otimes .$

The new directions

:$p_k = \left(1-P_k \right) u_k,$
:$p_k^* = v_k^* A \left(1- P_k \right) A^$

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 \left( 1- P_k \right) A^ r_j,$
:$r_k^*= r_j^* \left( 1- P_k \right)$

where i,j.

The biconjugate gradient method now makes a special choice and uses the setting
:$u_k = M^ r_k,\,$
:$v_k^* = r_k^* \, M^.\,$

With this particular choice, explicit evaluations of $P_k$ and are avoided, and the algorithm takes the form stated above.

Properties

* 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^*Ap_j=r_i^*M^r_j=0$ for $i \neq j$.
* if $P_\,$ is a polynomial with \deg\left(P_\right)+j, then $r_k^*P_\left(M^A\right)u_j=0$. The algorithm thus produces projections onto the Krylov subspace.
* if $P_\,$ is a polynomial with i+\deg\left(P_\right), then $v_i^*P_\left(AM^\right)r_k=0$.

See also

* Biconjugate gradient stabilized method
* Conjugate gradient method

References

*
*

{{Numerical linear algebra

Numerical linear algebra
Gradient methods