Let $x^i$ be the current approximate solution, and let $r^i=b-Ax^i$ be the residual. The algorithm finds a search direction $q^i$ and a step size $\alpha^i$ to update the approximate solution, that is, $x^{i+1}=x^i+\alpha^iq^i$, such that the search directions $q^i,\ldots,q^0$ are mutually $A$-conjugate \[i.e., $(q^i)'Aq^j = 0$ for $i\ne j$\] and the new residual $r^{i+1}$ is orthogonal to all the previous ones, $r^i,\ldots,r^0$. One can show that the search direction is a linear combination of the current residual and the past search direction, yielding the following recurrence formulas: $$\begin{aligned}
x^{i+1}&=&x^i+\alpha^iq^i,
\\
r^{i+1}&=&r^i-\alpha^iAq^i,
\\
q^{i+1}&=&r^{i+1}+\beta^iq^i,\end{aligned}$$ where $\alpha^i= \langle r^i,r^i \rangle/ \langle
Aq^i,q^i \rangle$ and $\beta^i= \langle r^{i+1},r^{i+1} \rangle/ \langle
r^i,r^i \rangle$, and $ \langle\cdot,\cdot\rangle$ denotes the vector inner product. Letting $x^*$ be the exact solution, that is, $Ax^*=b$, then $x^i$ enjoys a linear convergence to $x^*$: $$\label{eqncgconverge} \bigl\|x^i-x^*\bigr\|_A\le2 \biggl(
\frac{\sqrt{\kappa(A)}-1}{\sqrt{\kappa
(A)}+1} \biggr)^{i}\bigl\|x^0-x^*\bigr\|_A,$$ where $\|\cdot\|_A= \langle A\cdot,\cdot\rangle^{{1/2}}$ is the $A$-norm of a vector.

Asymptotically, the time cost of one iteration is upper bounded by that of multiplying $A$ by $q^i$, which typically dominates other vector operations when $A$ is not sparse. Properties of the covariance matrix can be exploited to efficiently compute the matrix–vector products. For example, when the observations are on a 