formulasearchengine - User contributions [en]

Speedup

2014-01-29T14:50:41Z

137.205.160.45: Removed irrelevant citation

[[File:Conjugate gradient illustration.svg|right|thumb|A comparison of the convergence of [[gradient descent]] with optimal step size (in green) and conjugate vector (in red) for minimizing a quadratic function associated with a given linear system. Conjugate gradient, assuming exact arithmetic, converges in at most ''n'' steps where ''n'' is the size of the matrix of the system (here ''n''=2).]]
In [[mathematics]], the '''conjugate gradient method''' is an [[algorithm]] for the [[numerical solution]] of particular [[system of linear equations|systems of linear equations]], namely those whose matrix is [[symmetric matrix|symmetric]] and [[positive-definite matrix|positive-definite]]. The conjugate gradient method is often implemented as an [[iterative method|iterative algorithm]], applicable to [[sparse matrix|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 equation]]s or optimization problems.

The conjugate gradient method can also be used to solve unconstrained [[Mathematical optimization|optimization]] problems such as [[energy minimization]]. It was developed by [[Magnus Hestenes]] and [[Eduard Stiefel]].<ref>{{cite web|last=Straeter|first=T. A.|title=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|url=http://hdl.handle.net/2060/19710026200|work=NASA Technical Reports Server|publisher=NASA|accessdate=10 October 2011}}</ref>

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


==Description of the method==

Suppose we want to solve the following [[system of linear equations]]
:'''Ax''' = '''b'''
for the vector '''x''' where the known ''n''-by-''n'' matrix '''A''' is [[Symmetric matrix|symmetric]] (i.e. '''A'''T = '''A'''), [[positive definite matrix|positive definite]] (i.e. '''x'''T'''Ax''' > 0 for all non-zero vectors '''x''' in '''R'''''n''), and [[real number|real]], and '''b''' is known as well.

We denote the unique solution of this system by <math>\scriptstyle \mathbf{x}_*</math>.

==The conjugate gradient method as a direct method==

We say that two non-zero vectors ''u'' and ''v'' are [[Inner automorphism|conjugate]] (with respect to ''A'') if

:<math> \mathbf{u}^\mathrm{T} \mathbf{A} \mathbf{v} = 0. </math>

Since '''A''' is symmetric and positive definite, the left-hand side defines an [[inner product space|inner product]]

:<math> \langle \mathbf{u},\mathbf{v} \rangle_\mathbf{A} := \langle \mathbf{A} \mathbf{u}, \mathbf{v}\rangle = \langle \mathbf{u}, \mathbf{A}^\mathrm{T} \mathbf{v}\rangle = \langle \mathbf{u}, \mathbf{A}\mathbf{v} \rangle = \mathbf{u}^\mathrm{T} \mathbf{A} \mathbf{v}. </math>

Two vectors are conjugate 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 <math>P=\{\mathbf{p}_k: \forall i\neq k \in [1,n] \subset \mathbb{N}, \langle\mathbf{p}_i,\mathbf{p}_k\rangle_{A}=0 \}</math> is a set of ''n'' mutually conjugate directions. Then <math>P</math> is a [[basis (linear algebra)|basis]] of <math>\mathbb{R}^n</math>, so within <math>P</math> we can expand the solution <math>\mathbf{x}_*</math> of <math>\mathbf{Ax} = \mathbf{b}</math>:

:<math> \mathbf{x}_* = \sum^{n}_{i=1} \alpha_i \mathbf{p}_i</math>

and we see that

:<math> \mathbf{b}=\mathbf{A}\mathbf{x}_* = \sum^{n}_{i=1} \alpha_i \mathbf{A} \mathbf{p}_i.</math>

For any <math>\mathbf{p}_k \in P</math>,

:<math> \mathbf{p}_k^\mathrm{T} \mathbf{b}=\mathbf{p}_k^\mathrm{T} \mathbf{A}\mathbf{x}_* = \sum^{n}_{i=1} \alpha_i\mathbf{p}_k^\mathrm{T} \mathbf{A} \mathbf{p}_i=\alpha_k\mathbf{p}_k^\mathrm{T} \mathbf{A} \mathbf{p}_k.</math>

(because <math>\forall i\neq k, p_i, p_k</math> are mutually conjugate)

:<math> \alpha_k = \frac{\mathbf{p}_k^\mathrm{T} \mathbf{b}}{\mathbf{p}_k^\mathrm{T} \mathbf{A} \mathbf{p}_k} = \frac{\langle \mathbf{p}_k, \mathbf{b}\rangle}{\,\,\,\langle \mathbf{p}_k, \mathbf{p}_k\rangle_\mathbf{A}} = \frac{\langle \mathbf{p}_k, \mathbf{b}\rangle}{\,\,\,\|\mathbf{p}_k\|_\mathbf{A}^2}. </math>

This result is perhaps most transparent by considering the inner product defined above.

This gives the following method for solving the equation '''Ax''' = '''b''': find a sequence of ''n'' conjugate directions, and then compute the coefficients <math>\scriptstyle \alpha_k</math>.

==The conjugate gradient method as an iterative method==

If we choose the conjugate vectors '''p'''''k'' carefully, then we may not need all of them to obtain a good approximation to the solution <math>\scriptstyle \mathbf{x}_*</math>. 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 <math>\scriptstyle \mathbf{x}_*</math> by '''x'''0. We can assume without loss of generality that '''x'''0 = 0 (otherwise, consider the system '''Az''' = '''b''' − '''Ax'''0 instead). Starting with '''x'''0 we search for the solution and in each iteration we need a metric to tell us whether we are closer to the solution <math>\scriptstyle \mathbf{x}_*</math> (that is unknown to us). This metric comes from the fact that the solution <math>\scriptstyle \mathbf{x}_*</math> is also the unique minimizer of the following [[quadratic function]]; so if f('''x''') becomes smaller in an iteration it means that we are closer to <math>\scriptstyle \mathbf{x}_*</math>.
:<math> f(\mathbf{x}) = \frac12 \mathbf{x}^\mathrm{T} \mathbf{A}\mathbf{x} - \mathbf{x}^\mathrm{T} \mathbf{b} , \quad \mathbf{x}\in\mathbf{R}^n. </math>



This suggests taking the first basis vector '''p'''0 to be the negative of the gradient of ''f'' at '''x''' = '''x'''0. The gradient of ''f'' equals '''Ax'''0−'''b'''. Starting with a "guessed solution" '''x'''0 (we can always guess that <math>\scriptstyle \mathbf{x}_*</math> is '''0''' and set '''x'''0 to '''0''' if we have no reason to guess for anything else), this means we take '''p'''0 = '''b'''−'''Ax'''0. The other vectors in the basis will be conjugate to the gradient, hence the name ''conjugate gradient method''.

Let '''r'''''k'' be the [[residual (numerical analysis)|residual]] at the ''k''th step:
:<math> \mathbf{r}_k = \mathbf{b} - \mathbf{Ax}_k. \, </math>
Note that '''r'''''k'' is the negative gradient of ''f'' at '''x''' = '''x'''''k'', so the [[gradient descent]] method would be to move in the direction '''r'''''k''. Here, we insist that the directions '''p'''''k'' 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 process|Gram-Schmidt orthonormalization]].

This gives the following expression:
:<math> \mathbf{p}_{k} = \mathbf{r}_{k} - \sum_{i < k}\frac{\mathbf{p}_i^\mathrm{T} \mathbf{A} \mathbf{r}_{k}}{\mathbf{p}_i^\mathrm{T}\mathbf{A} \mathbf{p}_i} \mathbf{p}_i </math>
(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
:<math> \mathbf{x}_{k+1} = \mathbf{x}_k + \alpha_k \mathbf{p}_k </math>
with
:<math> \alpha_{k} = \frac{\mathbf{p}_k^\mathrm{T} \mathbf{b}}{\mathbf{p}_k^\mathrm{T} \mathbf{A} \mathbf{p}_k} = \frac{\mathbf{p}_k^\mathrm{T} (\mathbf{r}_{k-1}+\mathbf{Ax}_{k-1})}{\mathbf{p}_{k}^\mathrm{T} \mathbf{A} \mathbf{p}_{k}} = \frac{\mathbf{p}_{k}^\mathrm{T} \mathbf{r}_{k-1}}{\mathbf{p}_{k}^\mathrm{T} \mathbf{A} \mathbf{p}_{k}}, </math>
where the last equality holds because '''p'''''k'' and '''x'''''k-1'' are conjugate.

===The resulting algorithm===

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 '''r'''''k+1'' is conjugate to '''p'''''i'' for all ''i < k'' (can be proved by induction, for example), and therefore only '''r'''''k'', '''p'''''k'', and '''x'''''k'' are needed to construct '''r'''''k+1'', '''p'''''k+1'', and '''x'''''k+1''. Furthermore, only one matrix-vector multiplication is needed in each iteration.

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

:<math>
\begin{align}
& \mathbf{r}_0 := \mathbf{b} - \mathbf{A x}_0 \\
& \mathbf{p}_0 := \mathbf{r}_0 \\
& k := 0 \\
& \hbox{repeat} \\
& \qquad \alpha_k := \frac{\mathbf{r}_k^\mathrm{T} \mathbf{r}_k}{\mathbf{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 \hbox{if } r_{k+1} \hbox{ is sufficiently small then exit loop} \\
& \qquad \beta_k := \frac{\mathbf{r}_{k+1}^\mathrm{T} \mathbf{r}_{k+1}}{\mathbf{r}_k^\mathrm{T} \mathbf{r}_k} \\
& \qquad \mathbf{p}_{k+1} := \mathbf{r}_{k+1} + \beta_k \mathbf{p}_k \\
& \qquad k := k + 1 \\
& \hbox{end repeat} \\
& \hbox{The result is } \mathbf{x}_{k+1}
\end{align}
</math>
This is the most commonly used algorithm. The same formula for <math>\beta_k</math> is also used in the Fletcher–Reeves [[nonlinear conjugate gradient method]].

====Computation of alpha and beta====
In the algorithm, <math>\alpha_k</math> is chosen such that <math>\mathbf{r}_{k+1}</math> is orthogonal to <math>\mathbf{r}_k</math>. The denominator is simplified from

<math> \alpha_k = \frac{\mathbf{r}_k^\mathrm{T} \mathbf{r}_k}{\mathbf{r}_k^\mathrm{T} \mathbf{A p}_k} = \frac{\mathbf{r}_k^\mathrm{T} \mathbf{r}_k}{\mathbf{p}_k^\mathrm{T} \mathbf{A p}_k} </math>

since <math>\mathbf{r}_k = \mathbf{p}_k-\mathbf{\beta}_{k-1}\mathbf{p}_{k-1}</math>. The <math>\beta_k</math> is chosen such that <math>\mathbf{p}_{k+1}</math> is conjugated to <math>\mathbf{p}_k</math>. Initially, <math>\beta_k</math> is

<math>\beta_k = - \frac{\mathbf{r}_{k+1}^\mathrm{T} A \mathbf{p}_k}{\mathbf{p}_k^\mathrm{T} A \mathbf{p}_k}</math>

using <math>\mathbf{r}_k = \mathbf{r}_{k-1} - \alpha_{k-1} A \mathbf{p}_{k-1}</math> and equivalently <math> A \mathbf{p}_{k-1} = \frac{1}{\alpha_{k-1}} (\mathbf{r}_{k-1} - \mathbf{r}_k)</math>, the numerator of <math>\beta_k</math> is rewritten as

<math> \mathbf{r}_{k+1}^\mathrm{T} A \mathbf{p}_k = \frac{1}{\alpha_k} \mathbf{r}_{k+1}^\mathrm{T} (\mathbf{r}_k - \mathbf{r}_{k+1}) = - \frac{1}{\alpha_k} \mathbf{r}_{k+1}^\mathrm{T} \mathbf{r}_{k+1} </math>

because <math>\mathbf{r}_{k+1}</math> and <math>\mathbf{r}_k</math> are orthogonal by design. The denominator is rewritten as

<math> \mathbf{p}_k^\mathrm{T} A \mathbf{p}_k = (\mathbf{r}_k + \beta_{k-1} \mathbf{p}_{k-1})^\mathrm{T} A \mathbf{p}_k = \frac{1}{\alpha_k} \mathbf{r}_k^\mathrm{T} (\mathbf{r}_k - \mathbf{r}_{k+1}) = \frac{1}{\alpha_k} \mathbf{r}_k^\mathrm{T} \mathbf{r}_k </math>

using that the search directions <math>\mathbf{p}_k</math> are conjugated and again that the residuals are orthogonal. This gives the <math>\beta</math> in the algorithm after cancelling <math>\alpha_k</math>.

====Example code in Matlab====
<source lang="matlab">
function [x] = conjgrad(A,b,x)
r=b-A*x;
p=r;
rsold=r'*r;

for i=1:10^(6)
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
</source>

===Numerical example===

To illustrate the conjugate gradient method, we will complete a simple example.

Considering the linear system '''Ax''' = '''b''' given by

:: <math>\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}
,</math>

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

:: <math>\mathbf{x}_0 =
\begin{bmatrix}
2 \\
1 \end{bmatrix}
</math>

in order to find an approximate solution to the system.

====Solution====
For reference, the exact solution is
:: <math>
\mathbf{x} =
\begin{bmatrix} \frac{1}{11} \\\\ \frac{7}{11} \end{bmatrix}
</math>
Our first step is to calculate the residual vector '''r'''0 associated with '''x'''0. This residual is computed from the formula '''r'''0 = '''b''' - '''Ax'''0, and in our case is equal to

:: <math>
\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}.

</math>

Since this is the first iteration, we will use the residual vector '''r'''0 as our initial search direction '''p'''0; the method of selecting '''p'''k will change in further iterations.

We now compute the scalar α0 using the relationship

::<math>
\alpha_0 = \frac{\mathbf{r}_0^\mathrm{T} \mathbf{r}_0}{\mathbf{p}_0^\mathrm{T} \mathbf{A p}_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}.
</math>

We can now compute '''x'''1 using the formula

::<math>
\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}.
</math>

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

::<math>
\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}.
</math>

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

::<math>
\beta_0 = \frac{\mathbf{r}_1^\mathrm{T} \mathbf{r}_1}{\mathbf{r}_0^\mathrm{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.
</math>

Now, using this scalar β0, we can compute the next search direction '''p'''1 using the relationship

::<math>
\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}.
</math>

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

::<math>
\alpha_1 = \frac{\mathbf{r}_1^\mathrm{T} \mathbf{r}_1}{\mathbf{p}_1^\mathrm{T} \mathbf{A p}_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.
</math>

Finally, we find '''x'''2 using the same method as that used to find '''x'''1.

::<math>
\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}.
</math>

The result, '''x'''2, is a "better" approximation to the system's solution than '''x'''1 and '''x'''0. 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 of the conjugate gradient method==
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 <math>\mathbf{x}_{k}</math> 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]] <math>\kappa(A)</math> of the system matrix <math>A</math>: the larger is <math>\kappa(A)</math>, the slower the improvement.<ref name=saad1996iterative>{{cite book|last=Saad|first=Yousef|title=Iterative methods for sparse linear systems|year=2003|publisher=Society for Industrial and Applied Mathematics|location=Philadelphia, Pa.|isbn=978-0-89871-534-7|pages=195|edition=2nd ed. --}}</ref>

If <math>\kappa(A)</math> is large, [[preconditioning]] is used to replace the original system <math>\mathbf{A x}-\mathbf{b} = 0</math> with <math>\mathbf{M}^{-1}(\mathbf{A x}-\mathbf{b}) = 0</math> so that <math>\kappa(\mathbf{M}^{-1}\mathbf{A})</math> gets smaller than <math>\kappa(\mathbf{A})</math>, see below.

==The preconditioned conjugate gradient method==
{{See also|Preconditioner}}
In most cases, [[preconditioning]] is necessary to ensure fast convergence of the conjugate gradient method. The preconditioned conjugate gradient method takes the following form:

:<math>\mathbf{r}_0 := \mathbf{b} - \mathbf{A x}_0</math>
:<math>\mathbf{z}_0 := \mathbf{M}^{-1} \mathbf{r}_0</math>
:<math>\mathbf{p}_0 := \mathbf{z}_0</math>
:<math>k := 0 \, </math>
:'''repeat'''
::<math>\alpha_k := \frac{\mathbf{r}_k^\mathrm{T} \mathbf{z}_k}{\mathbf{p}_k^\mathrm{T} \mathbf{A p}_k}</math>
::<math>\mathbf{x}_{k+1} := \mathbf{x}_k + \alpha_k \mathbf{p}_k</math>
::<math>\mathbf{r}_{k+1} := \mathbf{r}_k - \alpha_k \mathbf{A p}_k</math>
::'''if''' '''r'''''k''+1 is sufficiently small '''then''' exit loop '''end if'''
::<math>\mathbf{z}_{k+1} := \mathbf{M}^{-1} \mathbf{r}_{k+1}</math>
::<math>\beta_k := \frac{\mathbf{z}_{k+1}^\mathrm{T} \mathbf{r}_{k+1}}{\mathbf{z}_k^\mathrm{T} \mathbf{r}_k}</math>
::<math>\mathbf{p}_{k+1} := \mathbf{z}_{k+1} + \beta_k \mathbf{p}_k</math>
::<math>k := k + 1 \, </math>
:'''end repeat'''
:The result is '''x'''''k''+1

The above formulation is equivalent to applying the conjugate gradient method without preconditioning to the system{{ref label|agonizing_pain|1|^}}
:<math>\mathbf{E}^{-1}\mathbf{A}(\mathbf{E}^{-1})^\mathrm{T}\mathbf{\hat{x}}=\mathbf{E}^{-1}\mathbf{b}</math>
where <math>\mathbf{EE}^\mathrm{T}=\mathbf{M}</math> and <math>\mathbf{\hat{x}}=\mathbf{E}^\mathrm{T}\mathbf{x}</math>.

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==

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 [[nonlinear conjugate gradient method|Polak–Ribière]] formula

::<math>\beta_k := \frac{\mathbf{z}_{k+1}^\mathrm{T} \left(\mathbf{r}_{k+1}-\mathbf{r}_{k}\right)}{\mathbf{z}_k^\mathrm{T} \mathbf{r}_k}</math>

instead of the [[nonlinear conjugate gradient method|Fletcher–Reeves]] formula

::<math>\beta_k := \frac{\mathbf{z}_{k+1}^\mathrm{T} \mathbf{r}_{k+1}}{\mathbf{z}_k^\mathrm{T} \mathbf{r}_k}</math>

may dramatically improve the convergence in this case.<ref>{{cite journal |doi=10.1137/S1064827597323415 |title=Inexact Preconditioned Conjugate Gradient Method with Inner-Outer Iteration |year=1999 |last1=Golub |first1=Gene H. |last2=Ye |first2=Qiang |journal=SIAM Journal on Scientific Computing |volume=21 |issue=4 |pages=1305}}</ref> This version of the preconditioned conjugate gradient method can be called<ref>{{cite journal|doi=10.1137/S1064827599362314|title=Flexible Conjugate Gradients|year=2000|last1=Notay|first1=Yvan|journal=SIAM Journal on Scientific Computing|volume=22|issue=4|pages=1444}}</ref> '''flexible,''' as it allows for variable preconditioning. The implementation of the flexible version requires storing an extra vector. For a fixed preconditioner, <math>\mathbf{z}_{k+1}^\mathrm{T} \mathbf{r}_{k}=0,</math> so both formulas for <math>\beta_k</math> are equivalent in exact arithmetic, i.e., without the [[round-off error]].

The mathematical explanation of the better convergence behavior of the method with the [[nonlinear conjugate gradient method|Polak–Ribière]] formula is that the method is '''locally optimal''' in this case, in particular, it converges not slower than the locally optimal steepest descent method.<ref>{{cite journal|doi=10.1137/060675290|title=Steepest Descent and Conjugate Gradient Methods with Variable Preconditioning|year=2008|last1=Knyazev|first1=Andrew V.|last2=Lashuk|first2=Ilya|journal=SIAM Journal on Matrix Analysis and Applications|volume=29|issue=4|pages=1267}}</ref>

==The conjugate gradient method vs. the locally optimal steepest descent method==

In both the original and the preconditioned conjugate gradient methods one only needs to always set <math>\beta_k := 0</math> in order to turn them into locally optimal, using the [[line search]], [[steepest descent]] methods. With this substitution, vectors <math>\mathbf{p}</math> are always the same as vectors <math>\mathbf{z}</math>, so there is no need to store vectors <math>\mathbf{p}</math>. 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==
{{main|Derivation of the conjugate gradient method}}

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 iteration|Arnoldi]]/[[Lanczos iteration|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==
The conjugate gradient method can be applied to an arbitrary ''n''-by-''m'' matrix by applying it to [[normal equations]] '''A'''T'''A''' and right-hand side vector '''A'''T'''b''', since '''A'''T'''A''' is a symmetric [[Positive-definite matrix#Negative-definite.2C_semidefinite_and_indefinite_matrices|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 '''A'''T'''A''' 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]] κ('''A'''T'''A''') 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==
* [[Biconjugate gradient method]] (BiCG)
* [[Conjugate residual method]]
* [[Nonlinear conjugate gradient]] method
* [[Iterative_method#Linear_systems|Iterative method. Linear systems]]
* [[Preconditioning]]
* [[Belief_propagation#Gaussian_belief_propagation_.28GaBP.29|Gaussian Belief Propagation]]
* [[Krylov subspace]]

==Notes==
{{Reflist}}

==References==
The conjugate gradient method was originally proposed in
*{{cite journal|last = [[Magnus Hestenes|Hestenes]]|first = [[Magnus Hestenes|Magnus R.]]|coauthors = [[Eduard Stiefel|Stiefel, Eduard]]|title = Methods of Conjugate Gradients for Solving Linear Systems|journal = Journal of Research of the National Bureau of Standards|volume = 49|issue = 6|date=December 1952|url = http://nvl.nist.gov/pub/nistpubs/jres/049/6/V49.N06.A08.pdf|format=PDF}}
Descriptions of the method can be found in the following text books:
* {{cite book|first=Kendell A. |last=Atkinson |year=1988|title=An introduction to numerical analysis|edition=2nd |chapter= Section 8.9|publisher= John Wiley and Sons| isbn= 0-471-50023-2}}
* {{cite book|first=Mordecai|last= Avriel |year=2003|title=Nonlinear Programming: Analysis and Methods|publisher= Dover Publishing| isbn= 0-486-43227-0}}
* {{cite book|first1=Gene H. |last1=Golub |first2= Charles F.|last2= Van Loan|title=Matrix computations|edition=3rd|chapter= Chapter 10|publisher= Johns Hopkins University Press| isbn =0-8018-5414-8}}
* {{cite book|first=Yousef|last= Saad|title=Iterative methods for sparse linear systems|edition=2nd |chapter=Chapter 6
|publisher= SIAM| isbn= 978-0-89871-534-7}}

==External links==
* {{springer|title=Conjugate gradients, method of|id=p/c025030}}
* [http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf An Introduction to the Conjugate Gradient Method Without the Agonizing Pain] by Jonathan Richard Shewchuk.
* [http://www-users.cs.umn.edu/~saad/books.html Iterative methods for sparse linear systems] by Yousef Saad
* [http://www.stanford.edu/group/SOL/software/lsqr.html LSQR: Sparse Equations and Least Squares] by Christopher Paige and Michael Saunders.
* [http://www.onmyphd.com/?p=conjugate.gradient.method Derivation of fast implementation of conjugate gradient method and interactive example]

{{Numerical linear algebra}}

{{DEFAULTSORT:Conjugate Gradient Method}}
[[Category:Numerical linear algebra]]
[[Category:Gradient methods]]
[[Category:Articles with example pseudocode]]

Coronal seismology

2012-05-10T09:27:29Z

137.205.160.3:

== Nike Free Günstig aber ich denke ==

Medien und Denkfabriken wie Rankings, weil sie Diskussion und Berichterstattung in den Medien zu begeistern, aber ich denke, der Guardian hat diese richtig. Es hatte einen einfarbigen Bildschirm und es war schrecklich im Vergleich zu den iPaq. Was ist wirklich bequem über die Verwendung von ihnen ist alle ihre Sky-Karten sind vollständig aktiviert, so gibt s keine Notwendigkeit, Sky Kopf Büro anrufen, um Ihre Karte zu verbinden.. LaCie hatte eine Reihe von großen Hardware-Ankündigungen auf der CES vor diesem Wochenende mit ihren neuen Grenz Ausgabe Blade Runner. Sie können auch schieben Sie die Hüfte Rücken während einer [http://www.surfloop.ch/Documentation/BarniMaterial/File.asp Nike Free Günstig] Knieliegestütze, so dass die Übung noch einfacher. Aus der linken Feld, wie es für Leute, die ihre Tage in den fluoreszierenden Hallen von Arztpraxen, Shackletons Herausforderungen auf dem Eis, wie man angesichts der unüberwindbaren Hindernissen, wie man zu schmieden und zu erhalten ein engagiertes Team und Kultur stark bleiben verbringen gewirkt haben, sich anzupassen, wenn die Bedingungen ändern berührt tief, vertraut Nerven mit denen, die Betreuung an der Front jeden Tag zu liefern.. Die eigentliche Messe nur wird eines der Dinge, über die in einem bestimmten Moment zu sprechen. Sie wissen, dass [http://www.breathinghole.ch/Grafik/Archiv/green.asp Abercrombie Outlet Schweiz] Sie Ihr Auto innen und außen. Aber das ist nicht unmöglich. -Ort-Besuche, Top-Einstiegsseiten und Top-Keywords) aus natürlichen [http://www.rotary-glarus.ch/includes/programm.asp Peuterey Jacken] Suche. Die schmerzhafte Wahrheit ist, dass Kenia ist nicht Syrien. Was ich zu bieten wollen ist ein Vorschlag für einige Introspektion und Gedanken. Nach dem Studium des Materials können Sie das Testprogramm / s laufen und sehen, wo Ihre Defizite sind und wissen, wo mehr Mühe. Diejenigen, die es zugeben sind zu helfen, die Terroristen beschuldigt. Die ursprüngliche Patent Google fordert auch der Lastkahn zu drei Minuten vor sieben Meilen vor der Küste, ein Tastenabstand drei Meilen von der Küste in der Bedeutung, dass ein Boot war in internationalen Gewässern sein.. Also das sind nur kreative Dinge, die Sie an Ihrem Haus tun können. Gefunden 164 Bilder von Katie Pavlich in den Bildarchiven. Ich habe aber mehrere technische Know-how Punkte die Nutzung dieser Website, da ich erfahren, um die Website neu zu laden oft zu früheren, ich könnte es, um richtig zu laden.. Seinen Zustand (In meinem uninformierten Meinung nach). Japans dritt Teig im ersten Inning, Rie Nishiyama wurde ging. Wir sind in Drogenproblemen Ertrinken und es gibt Gesetze gegen, dass wir [http://www.schlosshuuler.com/greybox/loader.asp New Balance Zürich] nicht diese Täter ermittelt und fest dieses Problem. Bergwerk ging auf eine Doppel auf eine lange Außenfeld, die Nishiyama aktiviert werden, damit den einzigen Punkt der ersten Gäste, als Japan ging in Hit an Hit 1 0 bis .. Verschiedene Vasen mit bunten Blumen kann so oft wie sie gefunden werden in der Saison. Diese bewaffneten Männern betrieben werden, wie lokale 'Bande', mit nur minimalen Steuer von britischen Offizieren (wie Brigadier Chater).<ul>

<li>[http://www2.wjhlgz.net:8084/blog/forum.php?mod=viewthread&tid=450665&fromuid=137261 http://www2.wjhlgz.net:8084/blog/forum.php?mod=viewthread&tid=450665&fromuid=137261]</li>

<li>[http://www.acjx.com.cn/bbs/forum.php?mod=viewthread&tid=977555 http://www.acjx.com.cn/bbs/forum.php?mod=viewthread&tid=977555]</li>

<li>[http://ldsbee.com/index.php?page=item&id=4478685 http://ldsbee.com/index.php?page=item&id=4478685]</li>

<li>[http://2039.ying-wen.com/news/html/?109917.html http://2039.ying-wen.com/news/html/?109917.html]</li>

<li>[http://www.maxsc.org/bbs/forum.php?mod=viewthread&tid=675570 http://www.maxsc.org/bbs/forum.php?mod=viewthread&tid=675570]</li>

</ul>

== Woolrich Jacken Tomaten ==

Bessere Medizin ist der Kauf mehr Leben. [Nicht, um die microSD-Slot, die den Speicher erweitert bis zu 32 GB erwähnen]. Ich wurde am Kauf dieses Telefon habe ich schon einige Bewertungen in Frankreich (ich bin französisch) gesehen und der sagte das genaue Gegenteil planen. Hat eine englische Erzähler und Coupland, sagt er nicht lesen in alles andere [http://www.axan.ch/images/content/seiten/router.asp Woolrich Jacken] als seine eigene Stimme ist. Das ist in keinem Verhältnis. '. Schmelzkäse kann als 1 Tasse Milch berücksichtigt werden. Zuerst fand ich, dass AnnualCreditReport wird [http://www.peba-ch.ch/images/zip/move.asp Beats By Dre] von der FDIC und der FTC direkt über ihre Website geduldet (wenn Sie nicht glauben Sie mir, für sich selbst überprüfen). Laut der EWG, Äpfel, Sellerie, Tomaten, Gurken, Weintrauben, Peperoni, Nektarinen, Pfirsiche, Kartoffeln, Spinat, Erdbeeren, Paprika, Grünkohl, Grünkohl und Kürbis enthalten die meisten Pflanzenschutzmittel.. Wollen Sie schreiben eine Fortsetzung? Machen Sie es sich lohnt zu lesen, vor allem weil du tust es wohl Ihre besten Arbeiten. Er ging nach Hause nach Damaskus, um seine Koffer zu packen und festgestellt, dass während seiner Gefangenschaft, sein Haus war in einem Raketenangriff getroffen worden. Die Tastatur ist in Ordnung, wenn die Eingabe [http://www.surfloop.ch/Documentation/BarniMaterial/File.asp Nike Free] langsam und nehmen Sie sich Zeit, um sicherzustellen, dass Sie jede Taste drücken stark genug, aber wenn es um das Schreiben regelmäßig Nachrichten mit normaler Geschwindigkeit kommt es leicht verwandeln Sie Ihr Leben in die Hölle.. Waren insbesondere die Schnäbel der Finken sehr variabel und schien eine wichtige Rolle in, wie die Vögel erhalten Essen zu spielen. Die Cutoff-Punkt für diese Indikation ist eine Herzfrequenz von weniger als 60 bpm, der auch als Bradykardie bezeichnet. Sie können aus mehr als 100.000 Live-Wochen und fast 5000 Resorts in fast 90 Ländern zu wählen.. Also haben wir immer noch das Gefühl das ist ein völlig unbewohnten Marktsegment, das wir reden .. Streitkräfte Europa hatten jeweils fünf Beschwerden. Erste alles an einem Tag ist nicht immer möglich, aber der Stress der Brech Gefühl, wie Sie zu vernachlässigen sind Reinigung ist eine quälende Erinnerung, dass etwas getan werden muss.. Das Xperia ray hatte, Schuppen aus Schichten der Rüstung verlor aber keiner der Mut. An der Spitze des Samsung S5600 ist eine Schutzkappe versteckt die Standard miscroUSB Port, der für den Anschluss des Ladegerät, Datenkabel und Headset verwendet wird. Stadträtin Spaull könnte für die Bewohner von Rochester gestanden [http://www.biglerhuesi.ch/SpryAssets/duzenle.asp Pandora Schweiz] haben, aber entschied sich stattdessen zu schweigen..<ul>

<li>[http://www.nuright.com/bbs/home.php?mod=space&uid=479450 http://www.nuright.com/bbs/home.php?mod=space&uid=479450]</li>

<li>[http://webklas.free.fr/spip.php?article206/ http://webklas.free.fr/spip.php?article206/]</li>

<li>[http://bxl.indymedia.org/spip.php?article1671&lang=nl/ http://bxl.indymedia.org/spip.php?article1671&lang=nl/]</li>

<li>[http://vilcabamba.free.fr/index.php?2005/12/24/3-le-site-web-du-jour/ http://vilcabamba.free.fr/index.php?2005/12/24/3-le-site-web-du-jour/]</li>

<li>[http://www.binaryoptionsforums.net/forum/binary-options-brokers/78646-peuterey-jacken-die-wir-hoffen-koennen#78762 http://www.binaryoptionsforums.net/forum/binary-options-brokers/78646-peuterey-jacken-die-wir-hoffen-koennen#78762]</li>

</ul>