共轭梯度算法

Conjugate gradient method

In mathematics, the conjugate gradient method is an algorithm for the numerical solution of particular systems of linear equations, namely those whose matrix is symmetric and positive-definite. The conjugate gradient method is often implemented as an iterative algorithm, applicable to sparse systems that are too large to be handled by a direct implementation or other direct methods such as the Cholesky decomposition. Large sparse systems often arise when numerically solving partial differential equations or optimization problems.

The conjugate gradient method can also be used to solve unconstrained optimization problems such as energy minimization. It was mainly developed by Magnus Hestenes and Eduard Stiefel.[1]

The biconjugate gradient method provides a generalization to non-symmetric matrices. Various nonlinear conjugate gradient methods seek minima of nonlinear equations.

Description of the problem addressed by conjugate gradients[edit]

Suppose we want to solve the system of linear equations

{\displaystyle \mathbf {A} \mathbf {x} =\mathbf {b} }\mathbf {A} \mathbf {x} =\mathbf {b}

for the vector x, where the known n × n matrix A is symmetric (i.e., AT = A), positive-definite (i.e. xTAx > 0 for all non-zero vectors x in Rn), and real, and b is known as well. We denote the unique solution of this system by {\displaystyle \mathbf {x} _{*}}\mathbf{x}_*.

As a direct method[edit]

We say that two non-zero vectors u and v are conjugate (with respect to A) if

{\displaystyle \mathbf {u} ^{\mathsf {T}}\mathbf {A} \mathbf {v} =0.}\mathbf {u} ^{\mathsf {T}}\mathbf {A} \mathbf {v} =0.

Since A is symmetric and positive-definite, the left-hand side defines an inner product

{\displaystyle \langle \mathbf {u} ,\mathbf {v} \rangle _{\mathbf {A} }:=\langle \mathbf {A} \mathbf {u} ,\mathbf {v} \rangle =\langle \mathbf {u} ,\mathbf {A} ^{\mathsf {T}}\mathbf {v} \rangle =\langle \mathbf {u} ,\mathbf {A} \mathbf {v} \rangle =\mathbf {u} ^{\mathsf {T}}\mathbf {A} \mathbf {v} .}{\displaystyle \langle \mathbf {u} ,\mathbf {v} \rangle _{\mathbf {A} }:=\langle \mathbf {A} \mathbf {u} ,\mathbf {v} \rangle =\langle \mathbf {u} ,\mathbf {A} ^{\mathsf {T}}\mathbf {v} \rangle =\langle \mathbf {u} ,\mathbf {A} \mathbf {v} \rangle =\mathbf {u} ^{\mathsf {T}}\mathbf {A} \mathbf {v} .}

Two vectors are conjugate if and only if they are orthogonal with respect to this inner product. Being conjugate is a symmetric relation: if u is conjugate to v, then v is conjugate to u. Suppose that

{\displaystyle P=\{\mathbf {p} _{1},\dots ,\mathbf {p} _{n}\}}{\displaystyle P=\{\mathbf {p} _{1},\dots ,\mathbf {p} _{n}\}}

is a set of n mutually conjugate vectors (with respect to A). Then P forms a basis for {\displaystyle \mathbb {R} ^{n}}\mathbb {R} ^{n}, and we may express the solution x of {\displaystyle \mathbf {Ax} =\mathbf {b} }\mathbf {Ax} =\mathbf {b} in this basis:

{\displaystyle \mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {p} _{i}.}{\displaystyle \mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {p} _{i}.}

Based on this expansion we calculate:

{\displaystyle \mathbf {A} \mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {A} \mathbf {p} _{i}.}{\displaystyle \mathbf {A} \mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {A} \mathbf {p} _{i}.}

Left-multiplying by {\displaystyle \mathbf {p} _{k}^{\mathsf {T}}}{\displaystyle \mathbf {p} _{k}^{\mathsf {T}}}:

{\displaystyle \mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{i},}{\displaystyle \mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{i},}

substituting {\displaystyle \mathbf {Ax_{*}} =\mathbf {b} }{\displaystyle \mathbf {Ax_{*}} =\mathbf {b} } and {\displaystyle \mathbf {u} ^{\mathsf {T}}\mathbf {A} \mathbf {v} =\langle \mathbf {u} ,\mathbf {v} \rangle _{\mathbf {A} }}{\displaystyle \mathbf {u} ^{\mathsf {T}}\mathbf {A} \mathbf {v} =\langle \mathbf {u} ,\mathbf {v} \rangle _{\mathbf {A} }}:

{\displaystyle \mathbf {p} _{k}^{\mathsf {T}}\mathbf {b} =\sum _{i=1}^{n}\alpha _{i}\left\langle \mathbf {p} _{k},\mathbf {p} _{i}\right\rangle _{\mathbf {A} },}{\displaystyle \mathbf {p} _{k}^{\mathsf {T}}\mathbf {b} =\sum _{i=1}^{n}\alpha _{i}\left\langle \mathbf {p} _{k},\mathbf {p} _{i}\right\rangle _{\mathbf {A} },}

then {\displaystyle \mathbf {u} ^{\mathsf {T}}\mathbf {v} =\langle \mathbf {u} ,\mathbf {v} \rangle }{\displaystyle \mathbf {u} ^{\mathsf {T}}\mathbf {v} =\langle \mathbf {u} ,\mathbf {v} \rangle } and using {\displaystyle \forall i\neq k:\langle \mathbf {p} _{k},\mathbf {p} _{i}\rangle _{\mathbf {A} }=0}{\displaystyle \forall i\neq k:\langle \mathbf {p} _{k},\mathbf {p} _{i}\rangle _{\mathbf {A} }=0} yields

{\displaystyle \langle \mathbf {p} _{k},\mathbf {b} \rangle =\alpha _{k}\langle \mathbf {p} _{k},\mathbf {p} _{k}\rangle _{\mathbf {A} },}{\displaystyle \langle \mathbf {p} _{k},\mathbf {b} \rangle =\alpha _{k}\langle \mathbf {p} _{k},\mathbf {p} _{k}\rangle _{\mathbf {A} },}

which implies

{\displaystyle \alpha _{k}={\frac {\langle \mathbf {p} _{k},\mathbf {b} \rangle }{\langle \mathbf {p} _{k},\mathbf {p} _{k}\rangle _{\mathbf {A} }}}.}{\displaystyle \alpha _{k}={\frac {\langle \mathbf {p} _{k},\mathbf {b} \rangle }{\langle \mathbf {p} _{k},\mathbf {p} _{k}\rangle _{\mathbf {A} }}}.}

This gives the following method for solving the equation Ax = b: find a sequence of n conjugate directions, and then compute the coefficients αk.

As an iterative method[edit]

If we choose the conjugate vectors pk carefully, then we may not need all of them to obtain a good approximation to the solution x. So, we want to regard the conjugate gradient method as an iterative method. This also allows us to approximately solve systems where n is so large that the direct method would take too much time.

We denote the initial guess for x by x0. We can assume without loss of generality that x0 = 0 (otherwise, consider the system Az = b − Ax0 instead). Starting with x0we search for the solution and in each iteration we need a metric to tell us whether we are closer to the solution x (that is unknown to us). This metric comes from the fact that the solution x is also the unique minimizer of the following quadratic function

{\displaystyle f(\mathbf {x} )={\tfrac {1}{2}}\mathbf {x} ^{\mathsf {T}}\mathbf {A} \mathbf {x} -\mathbf {x} ^{\mathsf {T}}\mathbf {b} ,\qquad \mathbf {x} \in \mathbf {R} ^{n}\,.}{\displaystyle f(\mathbf {x} )={\tfrac {1}{2}}\mathbf {x} ^{\mathsf {T}}\mathbf {A} \mathbf {x} -\mathbf {x} ^{\mathsf {T}}\mathbf {b} ,\qquad \mathbf {x} \in \mathbf {R} ^{n}\,.}

The existence of a unique minimizer is apparent as its second derivative is given by a symmetric positive-definite matrix

{\displaystyle \mathrm {D} ^{2}f(\mathbf {x} )=\mathbf {A} \,,}{\displaystyle \mathrm {D} ^{2}f(\mathbf {x} )=\mathbf {A} \,,}

and that the minimizer (use Df(x)=0) solves the initial problem is obvious from its first derivative

{\displaystyle \mathrm {D} f(\mathbf {x} )=\mathbf {A} \mathbf {x} -\mathbf {b} \,.}{\displaystyle \mathrm {D} f(\mathbf {x} )=\mathbf {A} \mathbf {x} -\mathbf {b} \,.}

This suggests taking the first basis vector p0 to be the negative of the gradient of f at x = x0. The gradient of f equals Ax − b. Starting with a "guessed solution" x0 (we can always guess x0 = 0 if we have no reason to guess for anything else), this means we take p0 = b − Ax0. The other vectors in the basis will be conjugate to the gradient, hence the name conjugate gradient method.

Let rk be the residual at the kth step:

{\displaystyle \mathbf {r} _{k}=\mathbf {b} -\mathbf {Ax} _{k}.}\mathbf {r} _{k}=\mathbf {b} -\mathbf {Ax} _{k}.

Note that rk is the negative gradient of f at x = xk, so the gradient descent method would be to move in the direction rk. Here, we insist that the directions pk be conjugate to each other. We also require that the next search direction be built out of the current residue and all previous search directions, which is reasonable enough in practice.

The conjugation constraint is an orthonormal-type constraint and hence the algorithm bears resemblance to Gram-Schmidt orthonormalization.

This gives the following expression:

{\displaystyle \mathbf {p} _{k}=\mathbf {r} _{k}-\sum _{i<k}{\frac {\mathbf {p} _{i}^{\mathsf {T}}\mathbf {A} \mathbf {r} _{k}}{\mathbf {p} _{i}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{i}}}\mathbf {p} _{i}}\mathbf {p} _{k}=\mathbf {r} _{k}-\sum _{i<k}{\frac {\mathbf {p} _{i}^{\mathsf {T}}\mathbf {A} \mathbf {r} _{k}}{\mathbf {p} _{i}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{i}}}\mathbf {p} _{i}

(see the picture at the top of the article for the effect of the conjugacy constraint on convergence). Following this direction, the next optimal location is given by

{\displaystyle \mathbf {x} _{k+1}=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}}\mathbf {x} _{k+1}=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}

with

{\displaystyle \alpha _{k}={\frac {\mathbf {p} _{k}^{\mathsf {T}}\mathbf {b} }{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}={\frac {\mathbf {p} _{k}^{\mathsf {T}}(\mathbf {r} _{k}+\mathbf {Ax} _{k})}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}={\frac {\mathbf {p} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}},}{\displaystyle \alpha _{k}={\frac {\mathbf {p} _{k}^{\mathsf {T}}\mathbf {b} }{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}={\frac {\mathbf {p} _{k}^{\mathsf {T}}(\mathbf {r} _{k}+\mathbf {Ax} _{k})}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}={\frac {\mathbf {p} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}},}

where the last equality holds because pk and xk are conjugate. The expression for {\displaystyle \alpha _{k}}\alpha _{k} can be derived if one substitutes the expression for xk+1 into f and minimizing it w.r.t. {\displaystyle \alpha _{k}}\alpha _{k}

{\displaystyle {\begin{aligned}f(\mathbf {x} _{k+1})&=f(\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k})=:g(\alpha _{k})\\g'(\alpha _{k})&{\overset {!}{=}}0\quad \Rightarrow \quad \alpha _{k}={\frac {\mathbf {p} _{k}^{\mathsf {T}}\mathbf {b} }{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}\,.\end{aligned}}}{\displaystyle {\begin{aligned}f(\mathbf {x} _{k+1})&=f(\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k})=:g(\alpha _{k})\\g'(\alpha _{k})&{\overset {!}{=}}0\quad \Rightarrow \quad \alpha _{k}={\frac {\mathbf {p} _{k}^{\mathsf {T}}\mathbf {b} }{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}\,.\end{aligned}}}

The resulting algorithm[edit]

The above algorithm gives the most straightforward explanation of the conjugate gradient method. Seemingly, the algorithm as stated requires storage of all previous searching directions and residue vectors, as well as many matrix-vector multiplications, and thus can be computationally expensive. However, a closer analysis of the algorithm shows that ri is orthonormal to rj , i.e. {\displaystyle \mathbf {r} _{i}^{\mathsf {T}}\mathbf {r} _{j}=0}{\displaystyle \mathbf {r} _{i}^{\mathsf {T}}\mathbf {r} _{j}=0} , for i ≠ j. And pi is A-orthogonal to pj , i.e. {\displaystyle \mathbf {p} _{i}^{\mathsf {T}}A\mathbf {p} _{j}=0}{\displaystyle \mathbf {p} _{i}^{\mathsf {T}}A\mathbf {p} _{j}=0} , for i ≠ j. This can be regarded that as the algorithm progresses, pi and ri span the same Krylov subspace. Where ri form the orthogonal basis with respect to standard inner product, and pi form the orthogonal basis with respect to inner product induced by A. Therefore, xk can be regarded as the projection of x on the Krylov subspace.

The algorithm is detailed below for solving Ax = b where A is a real, symmetric, positive-definite matrix. The input vector x0 can be an approximate initial solution or 0. It is a different formulation of the exact procedure described above.

{\displaystyle {\begin{aligned}&\mathbf {r} _{0}:=\mathbf {b} -\mathbf {Ax} _{0}\\&\mathbf {p} _{0}:=\mathbf {r} _{0}\\&k:=0\\&{\text{repeat}}\\&\qquad \alpha _{k}:={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {Ap} _{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 {Ap} _{k}\\&\qquad {\hbox{if }}r_{k+1}{\text{ is sufficiently small, then exit loop}}\\&\qquad \beta _{k}:={\frac {\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {r} _{k+1}}{\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}}\\&\qquad \mathbf {p} _{k+1}:=\mathbf {r} _{k+1}+\beta _{k}\mathbf {p} _{k}\\&\qquad k:=k+1\\&{\text{end repeat}}\\&{\text{The result is }}\mathbf {x} _{k+1}\end{aligned}}}{\displaystyle {\begin{aligned}&\mathbf {r} _{0}:=\mathbf {b} -\mathbf {Ax} _{0}\\&\mathbf {p} _{0}:=\mathbf {r} _{0}\\&k:=0\\&{\text{repeat}}\\&\qquad \alpha _{k}:={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {Ap} _{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 {Ap} _{k}\\&\qquad {\hbox{if }}r_{k+1}{\text{ is sufficiently small, then exit loop}}\\&\qquad \beta _{k}:={\frac {\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {r} _{k+1}}{\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}}\\&\qquad \mathbf {p} _{k+1}:=\mathbf {r} _{k+1}+\beta _{k}\mathbf {p} _{k}\\&\qquad k:=k+1\\&{\text{end repeat}}\\&{\text{The result is }}\mathbf {x} _{k+1}\end{aligned}}}

This is the most commonly used algorithm. The same formula for βk is also used in the Fletcher–Reeves nonlinear conjugate gradient method.

Computation of alpha and beta[edit]

In the algorithm, αk is chosen such that {\displaystyle \mathbf {r} _{k+1}}\mathbf {r} _{k+1} is orthogonal to rk. The denominator is simplified from

{\displaystyle \alpha _{k}={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {r} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {Ap} _{k}}}}\alpha _{k}={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {r} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {Ap} _{k}}}

since {\displaystyle \mathbf {r} _{k+1}=\mathbf {p} _{k+1}-\mathbf {\beta } _{k}\mathbf {p} _{k}}{\displaystyle \mathbf {r} _{k+1}=\mathbf {p} _{k+1}-\mathbf {\beta } _{k}\mathbf {p} _{k}}. The βk is chosen such that {\displaystyle \mathbf {p} _{k+1}}\mathbf {p} _{k+1} is conjugated to pk. Initially, βk is

{\displaystyle \beta _{k}=-{\frac {\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}}\beta _{k}=-{\frac {\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}

using

{\displaystyle \mathbf {r} _{k+1}=\mathbf {r} _{k}-\alpha _{k}\mathbf {A} \mathbf {p} _{k}}{\displaystyle \mathbf {r} _{k+1}=\mathbf {r} _{k}-\alpha _{k}\mathbf {A} \mathbf {p} _{k}}

and equivalently

{\displaystyle \mathbf {A} \mathbf {p} _{k}={\frac {1}{\alpha _{k}}}(\mathbf {r} _{k}-\mathbf {r} _{k+1}),}{\displaystyle \mathbf {A} \mathbf {p} _{k}={\frac {1}{\alpha _{k}}}(\mathbf {r} _{k}-\mathbf {r} _{k+1}),}

the numerator of βk is rewritten as

{\displaystyle \mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}={\frac {1}{\alpha _{k}}}\mathbf {r} _{k+1}^{\mathsf {T}}(\mathbf {r} _{k}-\mathbf {r} _{k+1})=-{\frac {1}{\alpha _{k}}}\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {r} _{k+1}}\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}={\frac {1}{\alpha _{k}}}\mathbf {r} _{k+1}^{\mathsf {T}}(\mathbf {r} _{k}-\mathbf {r} _{k+1})=-{\frac {1}{\alpha _{k}}}\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {r} _{k+1}

because {\displaystyle \mathbf {r} _{k+1}}\mathbf {r} _{k+1} and rk are orthogonal by design. The denominator is rewritten as

{\displaystyle \mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}=(\mathbf {r} _{k}+\beta _{k-1}\mathbf {p} _{k-1})^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}={\frac {1}{\alpha _{k}}}\mathbf {r} _{k}^{\mathsf {T}}(\mathbf {r} _{k}-\mathbf {r} _{k+1})={\frac {1}{\alpha _{k}}}\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}=(\mathbf {r} _{k}+\beta _{k-1}\mathbf {p} _{k-1})^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}={\frac {1}{\alpha _{k}}}\mathbf {r} _{k}^{\mathsf {T}}(\mathbf {r} _{k}-\mathbf {r} _{k+1})={\frac {1}{\alpha _{k}}}\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}

using that the search directions pk are conjugated and again that the residuals are orthogonal. This gives the β in the algorithm after cancelling αk.

Example code in MATLAB / GNU Octave[edit]
function [x] = conjgrad(A, b, x)
    r = b - A * x;
    p = r;
    rsold = r' * r;

    for i = 1:length(b)
        Ap = A * p;
        alpha = rsold / (p' * Ap);
        x = x + alpha * p;
        r = r - alpha * Ap;
        rsnew = r' * r;
        if sqrt(rsnew) < 1e-10
              break;
        end
        p = r + (rsnew / rsold) * p;
        rsold = rsnew;
    end
end

Numerical example[edit]

Consider the linear system Ax = b given by

{\displaystyle \mathbf {A} \mathbf {x} ={\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}x_{1}\\x_{2}\end{bmatrix}}={\begin{bmatrix}1\\2\end{bmatrix}},}\mathbf {A} \mathbf {x} ={\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}x_{1}\\x_{2}\end{bmatrix}}={\begin{bmatrix}1\\2\end{bmatrix}},

we will perform two steps of the conjugate gradient method beginning with the initial guess

{\displaystyle \mathbf {x} _{0}={\begin{bmatrix}2\\1\end{bmatrix}}}\mathbf {x} _{0}={\begin{bmatrix}2\\1\end{bmatrix}}

in order to find an approximate solution to the system.

Solution[edit]

For reference, the exact solution is

{\displaystyle \mathbf {x} ={\begin{bmatrix}{\frac {1}{11}}\\\\{\frac {7}{11}}\end{bmatrix}}\approx {\begin{bmatrix}0.0909\\\\0.6364\end{bmatrix}}}{\displaystyle \mathbf {x} ={\begin{bmatrix}{\frac {1}{11}}\\\\{\frac {7}{11}}\end{bmatrix}}\approx {\begin{bmatrix}0.0909\\\\0.6364\end{bmatrix}}}

Our first step is to calculate the residual vector r0 associated with x0. This residual is computed from the formula r0 = b - Ax0, and in our case is equal to

{\displaystyle \mathbf {r} _{0}={\begin{bmatrix}1\\2\end{bmatrix}}-{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}2\\1\end{bmatrix}}={\begin{bmatrix}-8\\-3\end{bmatrix}}.}\mathbf {r} _{0}={\begin{bmatrix}1\\2\end{bmatrix}}-{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}2\\1\end{bmatrix}}={\begin{bmatrix}-8\\-3\end{bmatrix}}.

Since this is the first iteration, we will use the residual vector r0 as our initial search direction p0; the method of selecting pk will change in further iterations.

We now compute the scalar α0 using the relationship

{\displaystyle \alpha _{0}={\frac {\mathbf {r} _{0}^{\mathsf {T}}\mathbf {r} _{0}}{\mathbf {p} _{0}^{\mathsf {T}}\mathbf {Ap} _{0}}}={\frac {{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}{{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}}={\frac {73}{331}}.}\alpha _{0}={\frac {\mathbf {r} _{0}^{\mathsf {T}}\mathbf {r} _{0}}{\mathbf {p} _{0}^{\mathsf {T}}\mathbf {Ap} _{0}}}={\frac {​{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}{​{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}}={\frac {73}{331}}.

We can now compute x1 using the formula

{\displaystyle \mathbf {x} _{1}=\mathbf {x} _{0}+\alpha _{0}\mathbf {p} _{0}={\begin{bmatrix}2\\1\end{bmatrix}}+{\frac {73}{331}}{\begin{bmatrix}-8\\-3\end{bmatrix}}={\begin{bmatrix}0.2356\\0.3384\end{bmatrix}}.}\mathbf {x} _{1}=\mathbf {x} _{0}+\alpha _{0}\mathbf {p} _{0}={\begin{bmatrix}2\\1\end{bmatrix}}+{\frac {73}{331}}{\begin{bmatrix}-8\\-3\end{bmatrix}}={\begin{bmatrix}0.2356\\0.3384\end{bmatrix}}.

This result completes the first iteration, the result being an "improved" approximate solution to the system, x1. We may now move on and compute the next residual vector r1 using the formula

{\displaystyle \mathbf {r} _{1}=\mathbf {r} _{0}-\alpha _{0}\mathbf {A} \mathbf {p} _{0}={\begin{bmatrix}-8\\-3\end{bmatrix}}-{\frac {73}{331}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}={\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}.}\mathbf {r} _{1}=\mathbf {r} _{0}-\alpha _{0}\mathbf {A} \mathbf {p} _{0}={\begin{bmatrix}-8\\-3\end{bmatrix}}-{\frac {73}{331}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}={\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}.

Our next step in the process is to compute the scalar β0 that will eventually be used to determine the next search direction p1.

{\displaystyle \beta _{0}={\frac {\mathbf {r} _{1}^{\mathsf {T}}\mathbf {r} _{1}}{\mathbf {r} _{0}^{\mathsf {T}}\mathbf {r} _{0}}}={\frac {{\begin{bmatrix}-0.2810&0.7492\end{bmatrix}}{\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}}{{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}}=0.0088.}\beta _{0}={\frac {\mathbf {r} _{1}^{\mathsf {T}}\mathbf {r} _{1}}{\mathbf {r} _{0}^{\mathsf {T}}\mathbf {r} _{0}}}={\frac {​{\begin{bmatrix}-0.2810&0.7492\end{bmatrix}}{\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}}{​{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}}=0.0088.

Now, using this scalar β0, we can compute the next search direction p1 using the relationship

{\displaystyle \mathbf {p} _{1}=\mathbf {r} _{1}+\beta _{0}\mathbf {p} _{0}={\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}+0.0088{\begin{bmatrix}-8\\-3\end{bmatrix}}={\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}.}\mathbf {p} _{1}=\mathbf {r} _{1}+\beta _{0}\mathbf {p} _{0}={\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}+0.0088{\begin{bmatrix}-8\\-3\end{bmatrix}}={\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}.

We now compute the scalar α1 using our newly acquired p1 using the same method as that used for α0.

{\displaystyle \alpha _{1}={\frac {\mathbf {r} _{1}^{\mathsf {T}}\mathbf {r} _{1}}{\mathbf {p} _{1}^{\mathsf {T}}\mathbf {Ap} _{1}}}={\frac {{\begin{bmatrix}-0.2810&0.7492\end{bmatrix}}{\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}}{{\begin{bmatrix}-0.3511&0.7229\end{bmatrix}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}}}=0.4122.}\alpha _{1}={\frac {\mathbf {r} _{1}^{\mathsf {T}}\mathbf {r} _{1}}{\mathbf {p} _{1}^{\mathsf {T}}\mathbf {Ap} _{1}}}={\frac {​{\begin{bmatrix}-0.2810&0.7492\end{bmatrix}}{\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}}{​{\begin{bmatrix}-0.3511&0.7229\end{bmatrix}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}}}=0.4122.

Finally, we find x2 using the same method as that used to find x1.

{\displaystyle \mathbf {x} _{2}=\mathbf {x} _{1}+\alpha _{1}\mathbf {p} _{1}={\begin{bmatrix}0.2356\\0.3384\end{bmatrix}}+0.4122{\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}={\begin{bmatrix}0.0909\\0.6364\end{bmatrix}}.}\mathbf {x} _{2}=\mathbf {x} _{1}+\alpha _{1}\mathbf {p} _{1}={\begin{bmatrix}0.2356\\0.3384\end{bmatrix}}+0.4122{\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}={\begin{bmatrix}0.0909\\0.6364\end{bmatrix}}.

The result, x2, is a "better" approximation to the system's solution than x1 and x0. If exact arithmetic were to be used in this example instead of limited-precision, then the exact solution would theoretically have been reached after n = 2 iterations (n being the order of the system).

Convergence properties[edit]

The conjugate gradient method can theoretically be viewed as a direct method, as it produces the exact solution after a finite number of iterations, which is not larger than the size of the matrix, in the absence of round-off error. However, the conjugate gradient method is unstable with respect to even small perturbations, e.g., most directions are not in practice conjugate, and the exact solution is never obtained. Fortunately, the conjugate gradient method can be used as an iterative method as it provides monotonically improving approximations {\displaystyle \mathbf {x} _{k}}\mathbf {x} _{k} to the exact solution, which may reach the required tolerance after a relatively small (compared to the problem size) number of iterations. The improvement is typically linear and its speed is determined by the condition number {\displaystyle \kappa (A)}\kappa (A) of the system matrix {\displaystyle A}A: the larger {\displaystyle \kappa (A)}\kappa (A) is, the slower the improvement.[2]

If {\displaystyle \kappa (A)}\kappa (A) is large, preconditioning is used to replace the original system {\displaystyle \mathbf {Ax} -\mathbf {b} =0}\mathbf {Ax} -\mathbf {b} =0 with {\displaystyle \mathbf {M} ^{-1}(\mathbf {Ax} -\mathbf {b} )=0}\mathbf {M} ^{-1}(\mathbf {Ax} -\mathbf {b} )=0 such that {\displaystyle \kappa (\mathbf {M} ^{-1}\mathbf {A} )}\kappa (\mathbf {M} ^{-1}\mathbf {A} ) is smaller than {\displaystyle \kappa (\mathbf {A} )}\kappa (\mathbf {A} ), see below.

Convergence theorem[edit]

Define a subset of polynomials as

{\displaystyle \Pi _{k}^{*}:=\left\lbrace \ p\in \Pi _{k}\ :\ p(0)=1\ \right\rbrace \,,}{\displaystyle \Pi _{k}^{*}:=\left\lbrace \ p\in \Pi _{k}\ :\ p(0)=1\ \right\rbrace \,,}

where {\displaystyle \Pi _{k}}{\displaystyle \Pi _{k}} is the set of polynomials of maximal degree {\displaystyle k}k.

Let {\displaystyle \left(\mathbf {x} _{k}\right)_{k}}{\displaystyle \left(\mathbf {x} _{k}\right)_{k}} be the iterative approximations of the exact solution {\displaystyle \mathbf {x} _{*}}{\displaystyle \mathbf {x} _{*}}, and define the errors as {\displaystyle \mathbf {e} _{k}:=\mathbf {x} _{k}-\mathbf {x} _{*}}{\displaystyle \mathbf {e} _{k}:=\mathbf {x} _{k}-\mathbf {x} _{*}}. Now, the rate of convergence can be approximated as [3]

{\displaystyle {\begin{aligned}\left\|\mathbf {e} _{k}\right\|_{\mathbf {A} }&=\min _{p\in \Pi _{k}^{*}}\left\|p(\mathbf {A} )\mathbf {e} _{0}\right\|_{\mathbf {A} }\\&\leq \min _{p\in \Pi _{k}^{*}}\,\max _{\lambda \in \sigma (\mathbf {A} )}|p(\lambda )|\ \left\|\mathbf {e} _{0}\right\|_{\mathbf {A} }\\&\leq 2\left({\frac {{\sqrt {\kappa (\mathbf {A} )}}-1}{{\sqrt {\kappa (\mathbf {A} )}}+1}}\right)^{k}\ \left\|\mathbf {e} _{0}\right\|_{\mathbf {A} }\,,\end{aligned}}}{\displaystyle {\begin{aligned}\left\|\mathbf {e} _{k}\right\|_{\mathbf {A} }&=\min _{p\in \Pi _{k}^{*}}\left\|p(\mathbf {A} )\mathbf {e} _{0}\right\|_{\mathbf {A} }\\&\leq \min _{p\in \Pi _{k}^{*}}\,\max _{\lambda \in \sigma (\mathbf {A} )}|p(\lambda )|\ \left\|\mathbf {e} _{0}\right\|_{\mathbf {A} }\\&\leq 2\left({\frac {​{\sqrt {\kappa (\mathbf {A} )}}-1}{​{\sqrt {\kappa (\mathbf {A} )}}+1}}\right)^{k}\ \left\|\mathbf {e} _{0}\right\|_{\mathbf {A} }\,,\end{aligned}}}

where {\displaystyle \sigma (\mathbf {A} )}{\displaystyle \sigma (\mathbf {A} )} denotes the spectrum, and {\displaystyle \kappa (\mathbf {A} )}{\displaystyle \kappa (\mathbf {A} )} denotes the condition number.

Note, the important limit when {\displaystyle \kappa (\mathbf {A} )}{\displaystyle \kappa (\mathbf {A} )} tends to {\displaystyle \infty }\infty

{\displaystyle {\frac {{\sqrt {\kappa (\mathbf {A} )}}-1}{{\sqrt {\kappa (\mathbf {A} )}}+1}}\approx 1-{\frac {2}{\sqrt {\kappa (\mathbf {A} )}}}\quad {\text{for}}\quad \kappa (\mathbf {A} )\gg 1\,.}{\displaystyle {\frac {​{\sqrt {\kappa (\mathbf {A} )}}-1}{​{\sqrt {\kappa (\mathbf {A} )}}+1}}\approx 1-{\frac {2}{\sqrt {\kappa (\mathbf {A} )}}}\quad {\text{for}}\quad \kappa (\mathbf {A} )\gg 1\,.}

This limit shows a faster convergence rate compared to the iterative methods of Jacobi or Gauss-Seidel which scale as {\displaystyle \approx 1-{\frac {2}{\kappa (\mathbf {A} )}}}{\displaystyle \approx 1-{\frac {2}{\kappa (\mathbf {A} )}}}.

The preconditioned conjugate gradient method[edit]

In most cases, preconditioning is necessary to ensure fast convergence of the conjugate gradient method. The preconditioned conjugate gradient method takes the following form:

{\displaystyle \mathbf {r} _{0}:=\mathbf {b} -\mathbf {Ax} _{0}}\mathbf {r} _{0}:=\mathbf {b} -\mathbf {Ax} _{0}
{\displaystyle \mathbf {z} _{0}:=\mathbf {M} ^{-1}\mathbf {r} _{0}}\mathbf {z} _{0}:=\mathbf {M} ^{-1}\mathbf {r} _{0}
{\displaystyle \mathbf {p} _{0}:=\mathbf {z} _{0}}\mathbf {p} _{0}:=\mathbf {z} _{0}
{\displaystyle k:=0\,}k:=0\,
repeat
{\displaystyle \alpha _{k}:={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {z} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {Ap} _{k}}}}\alpha _{k}:={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {z} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {Ap} _{k}}}
{\displaystyle \mathbf {x} _{k+1}:=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}}\mathbf {x} _{k+1}:=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}
{\displaystyle \mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{k}\mathbf {Ap} _{k}}\mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{k}\mathbf {Ap} _{k}
if r k+1 is sufficiently small  then exit loop  end if
{\displaystyle \mathbf {z} _{k+1}:=\mathbf {M} ^{-1}\mathbf {r} _{k+1}}\mathbf {z} _{k+1}:=\mathbf {M} ^{-1}\mathbf {r} _{k+1}
{\displaystyle \beta _{k}:={\frac {\mathbf {z} _{k+1}^{\mathsf {T}}\mathbf {r} _{k+1}}{\mathbf {z} _{k}^{\mathsf {T}}\mathbf {r} _{k}}}}\beta _{k}:={\frac {\mathbf {z} _{k+1}^{\mathsf {T}}\mathbf {r} _{k+1}}{\mathbf {z} _{k}^{\mathsf {T}}\mathbf {r} _{k}}}
{\displaystyle \mathbf {p} _{k+1}:=\mathbf {z} _{k+1}+\beta _{k}\mathbf {p} _{k}}\mathbf {p} _{k+1}:=\mathbf {z} _{k+1}+\beta _{k}\mathbf {p} _{k}
{\displaystyle k:=k+1\,}k:=k+1\,
end repeat
The result is  x k+1

The above formulation is equivalent to applying the conjugate gradient method without preconditioning to the system[1]

{\displaystyle \mathbf {E} ^{-1}\mathbf {A} (\mathbf {E} ^{-1})^{\mathsf {T}}\mathbf {\hat {x}} =\mathbf {E} ^{-1}\mathbf {b} }\mathbf {E} ^{-1}\mathbf {A} (\mathbf {E} ^{-1})^{\mathsf {T}}\mathbf {\hat {x}} =\mathbf {E} ^{-1}\mathbf {b}

where

{\displaystyle \mathbf {EE} ^{\mathsf {T}}=\mathbf {M} ,\qquad \mathbf {\hat {x}} =\mathbf {E} ^{\mathsf {T}}\mathbf {x} .}\mathbf {EE} ^{\mathsf {T}}=\mathbf {M} ,\qquad \mathbf {\hat {x}} =\mathbf {E} ^{\mathsf {T}}\mathbf {x} .

The preconditioner matrix M has to be symmetric positive-definite and fixed, i.e., cannot change from iteration to iteration. If any of these assumptions on the preconditioner is violated, the behavior of the preconditioned conjugate gradient method may become unpredictable.

An example of a commonly used preconditioner is the incomplete Cholesky factorization.

The flexible preconditioned conjugate gradient method[edit]

In numerically challenging applications, sophisticated preconditioners are used, which may lead to variable preconditioning, changing between iterations. Even if the preconditioner is symmetric positive-definite on every iteration, the fact that it may change makes the arguments above invalid, and in practical tests leads to a significant slow down of the convergence of the algorithm presented above. Using the Polak–Ribière formula

{\displaystyle \beta _{k}:={\frac {\mathbf {z} _{k+1}^{\mathsf {T}}\left(\mathbf {r} _{k+1}-\mathbf {r} _{k}\right)}{\mathbf {z} _{k}^{\mathsf {T}}\mathbf {r} _{k}}}}\beta _{k}:={\frac {\mathbf {z} _{k+1}^{\mathsf {T}}\left(\mathbf {r} _{k+1}-\mathbf {r} _{k}\right)}{\mathbf {z} _{k}^{\mathsf {T}}\mathbf {r} _{k}}}

instead of the Fletcher–Reeves formula

{\displaystyle \beta _{k}:={\frac {\mathbf {z} _{k+1}^{\mathsf {T}}\mathbf {r} _{k+1}}{\mathbf {z} _{k}^{\mathsf {T}}\mathbf {r} _{k}}}}\beta _{k}:={\frac {\mathbf {z} _{k+1}^{\mathsf {T}}\mathbf {r} _{k+1}}{\mathbf {z} _{k}^{\mathsf {T}}\mathbf {r} _{k}}}

may dramatically improve the convergence in this case.[4] This version of the preconditioned conjugate gradient method can be called[5] flexible, as it allows for variable preconditioning. The implementation of the flexible version requires storing an extra vector. For a fixed preconditioner, {\displaystyle \mathbf {z} _{k+1}^{\mathsf {T}}\mathbf {r} _{k}=0,}\mathbf {z} _{k+1}^{\mathsf {T}}\mathbf {r} _{k}=0, so both formulas for βkare equivalent in exact arithmetic, i.e., without the round-off error.

The mathematical explanation of the better convergence behavior of the method with the Polak–Ribière formula is that the method is locally optimal in this case, in particular, it does not converge slower than the locally optimal steepest descent method.[6]

Example code in MATLAB / GNU Octave[edit]

function [x, k] = cgp(x0, A, C, b, mit, stol, bbA, bbC)
% Synopsis:
% x0: initial point
% A: Matrix A of the system Ax=b
% C: Preconditioning Matrix can be left or right
% mit: Maximum number of iterations
% stol: residue norm tolerance
% bbA: Black Box that computes the matrix-vector product for A * u
% bbC: Black Box that computes:
%      for left-side preconditioner : ha = C \ ra
%      for right-side preconditioner: ha = C * ra
% x: Estimated solution point
% k: Number of iterations done 
%
% Example:
% tic;[x, t] = cgp(x0, S, speye(1), b, 3000, 10^-8, @(Z, o) Z*o, @(Z, o) o);toc
% Elapsed time is 0.550190 seconds.
%
% Reference:
%  Métodos iterativos tipo Krylov para sistema lineales
%  B. Molina y M. Raydan - {{ISBN|908-261-078-X}}
        if (nargin < 8), error('Not enough input arguments. Try help.'); end;
        if (isempty(A)), error('Input matrix A must not be empty.'); end;
        if (isempty(C)), error('Input preconditioner matrix C must not be empty.'); end;
        x = x0;
        ha = 0;
        hp = 0;
        hpp = 0;
        ra = 0;
        rp = 0;
        rpp = 0;
        u = 0;
        k = 0;

        ra = b - bbA(A, x0); % <--- ra = b - A * x0;
        while (norm(ra, inf) > stol),
                ha = bbC(C, ra); % <--- ha = C \ ra;
                k = k + 1;
                if (k == mit), warning('GCP:MAXIT', 'mit reached, no conversion.'); return; end;
                hpp = hp;
                rpp = rp;
                hp = ha;
                rp = ra;
                t = rp' * hp;
                if (k == 1),
                        u = hp;
                else
                        u = hp + (t / (rpp' * hpp)) * u;
                end;
                Au = bbA(A, u); % <--- Au = A * u;
                a = t / (u' * Au);
                x = x + a * u;
                ra = rp - a * Au;
        end;

Vs. the locally optimal steepest descent method[edit]

In both the original and the preconditioned conjugate gradient methods one only needs to set {\displaystyle \beta _{k}:=0}\beta _{k}:=0 in order to make them locally optimal, using the line searchsteepest descent methods. With this substitution, vectors p are always the same as vectors z, so there is no need to store vectors p. Thus, every iteration of these steepest descent methods is a bit cheaper compared to that for the conjugate gradient methods. However, the latter converge faster, unless a (highly) variable preconditioner is used, see above.

Derivation of the method[edit]

The conjugate gradient method can be derived from several different perspectives, including specialization of the conjugate direction method for optimization, and variation of the Arnoldi/Lanczos iteration for eigenvalue problems. Despite differences in their approaches, these derivations share a common topic—proving the orthogonality of the residuals and conjugacy of the search directions. These two properties are crucial to developing the well-known succinct formulation of the method.

Conjugate gradient on the normal equations[edit]

The conjugate gradient method can be applied to an arbitrary n-by-m matrix by applying it to normal equations ATA and right-hand side vector ATb, since ATA is a symmetric positive-semidefinite matrix for any A. The result is conjugate gradient on the normal equations (CGNR).

A T Ax =  A T b

As an iterative method, it is not necessary to form ATA explicitly in memory but only to perform the matrix-vector and transpose matrix-vector multiplications. Therefore, CGNR is particularly useful when A is a sparse matrix since these operations are usually extremely efficient. However the downside of forming the normal equations is that the condition number κ(ATA) is equal to κ2(A) and so the rate of convergence of CGNR may be slow and the quality of the approximate solution may be sensitive to roundoff errors. Finding a good preconditioner is often an important part of using the CGNR method.

Several algorithms have been proposed (e.g., CGLS, LSQR). The LSQR algorithm purportedly has the best numerical stability when A is ill-conditioned, i.e., A has a large condition number.

See also[edit]

Notes[edit]

  1. Jump up^ Straeter, T. A. "On the Extension of the Davidon-Broyden Class of Rank One, Quasi-Newton Minimization Methods to an Infinite Dimensional Hilbert Space with Applications to Optimal Control Problems"NASA Technical Reports Server. NASA. Retrieved 10 October 2011.
  2. Jump up^ Saad, Yousef (2003). Iterative methods for sparse linear systems (2nd ed.). Philadelphia, Pa.: Society for Industrial and Applied Mathematics. p. 195. ISBN 978-0-89871-534-7.
  3. Jump up^ 1948-, Hackbusch, W.,. Iterative solution of large sparse systems of equations (Second edition ed.). Switzerland. ISBN 9783319284835OCLC 952572240.
  4. Jump up^ Golub, Gene H.; Ye, Qiang (1999). "Inexact Preconditioned Conjugate Gradient Method with Inner-Outer Iteration". SIAM Journal on Scientific Computing21 (4): 1305. doi:10.1137/S1064827597323415.
  5. Jump up^ Notay, Yvan (2000). "Flexible Conjugate Gradients". SIAM Journal on Scientific Computing22 (4): 1444. doi:10.1137/S1064827599362314.
  6. Jump up^ Knyazev, Andrew V.; Lashuk, Ilya (2008). "Steepest Descent and Conjugate Gradient Methods with Variable Preconditioning". SIAM Journal on Matrix Analysis and Applications29 (4): 1267. arXiv:math/0605767Freely accessibledoi:10.1137/060675290.

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值