conjugate residual method

The conjugate residual method is an iterative numeric method used for solving systems of linear equations. It's a Krylov subspace method very similar to the much more popular conjugate gradient method, with similar construction and convergence properties.

This method is used to solve linear equations of the form

:$\backslash mathbf\; A\; \backslash mathbf\; x\; =\; \backslash mathbf\; b$

where A is an invertible and Hermitian matrix, and b is nonzero.

The conjugate residual method differs from the closely related conjugate gradient method. It involves more numerical operations and requires more storage.

Given an (arbitrary) initial estimate of the solution $\backslash mathbf\; x\_0$, the method is outlined below:

:$$

\begin{align}

& \mathbf{x}_0 := \text{Some initial guess} \\

& \mathbf{r}_0 := \mathbf{b} - \mathbf{A x}_0 \\

& \mathbf{p}_0 := \mathbf{r}_0 \\

& \text{Iterate, with } k \text{ starting at } 0:\\

& \qquad \alpha_k := \frac{\mathbf{r}_k^\mathrm{T} \mathbf{A r}_k}{(\mathbf{A p}_k)^\mathrm{T} \mathbf{A p}_k} \\

& \qquad \mathbf{x}_{k+1} := \mathbf{x}_k + \alpha_k \mathbf{p}_k \\

& \qquad \mathbf{r}_{k+1} := \mathbf{r}_k - \alpha_k \mathbf{A p}_k \\

& \qquad \beta_k := \frac{\mathbf{r}_{k+1}^\mathrm{T} \mathbf{A r}_{k+1}}{\mathbf{r}_k^\mathrm{T} \mathbf{A r}_k} \\

& \qquad \mathbf{p}_{k+1} := \mathbf{r}_{k+1} + \beta_k \mathbf{p}_k \\

& \qquad \mathbf{A p}_{k + 1} := \mathbf{A r}_{k+1} + \beta_k \mathbf{A p}_k \\

& \qquad k := k + 1

\end{align}

the iteration may be stopped once $\backslash mathbf\; x\_k$ has been deemed converged. The only difference between this and the conjugate gradient method is the calculation of $\backslash alpha\_k$ and $\backslash beta\_k$ (plus the optional incremental calculation of $\backslash mathbf\{A\; p\}\_k$ at the end).

Note: the above algorithm can be transformed so to make only one symmetric matrix-vector multiplication in each iteration.