HOME

TheInfoList



OR:

In mathematics, the conjugate gradient method is an
algorithm In mathematics and computer science, an algorithm () is a finite sequence of rigorous instructions, typically used to solve a class of specific problems or to perform a computation. Algorithms are used as specifications for performing ...
for the
numerical solution Numerical analysis is the study of algorithms that use numerical approximation (as opposed to symbolic manipulations) for the problems of mathematical analysis (as distinguished from discrete mathematics). It is the study of numerical methods th ...
of particular
systems of linear equations In mathematics, a system of linear equations (or linear system) is a collection of one or more linear equations involving the same variables. For example, :\begin 3x+2y-z=1\\ 2x-2y+4z=-2\\ -x+\fracy-z=0 \end is a system of three equations in t ...
, namely those whose matrix is
positive-definite In mathematics, positive definiteness is a property of any object to which a bilinear form or a sesquilinear form may be naturally associated, which is positive-definite. See, in particular: * Positive-definite bilinear form * Positive-definite ...
. The conjugate gradient method is often implemented as an
iterative algorithm In computational mathematics, an iterative method is a mathematical procedure that uses an initial value to generate a sequence of improving approximate solutions for a class of problems, in which the ''n''-th approximation is derived from the pr ...
, applicable to sparse systems that are too large to be handled by a direct implementation or other direct methods such as the
Cholesky decomposition In linear algebra, the Cholesky decomposition or Cholesky factorization (pronounced ) is a decomposition of a Hermitian, positive-definite matrix into the product of a lower triangular matrix and its conjugate transpose, which is useful for effic ...
. Large sparse systems often arise when numerically solving
partial differential equation In mathematics, a partial differential equation (PDE) is an equation which imposes relations between the various partial derivatives of a multivariable function. The function is often thought of as an "unknown" to be solved for, similarly to ...
s or optimization problems. The conjugate gradient method can also be used to solve unconstrained
optimization Mathematical optimization (alternatively spelled ''optimisation'') or mathematical programming is the selection of a best element, with regard to some criterion, from some set of available alternatives. It is generally divided into two subfi ...
problems such as energy minimization. It is commonly attributed to Magnus Hestenes and Eduard Stiefel, who programmed it on the Z4, and extensively researched it. The biconjugate gradient method provides a generalization to non-symmetric matrices. Various nonlinear conjugate gradient methods seek minima of nonlinear optimization problems.


Description of the problem addressed by conjugate gradients

Suppose we want to solve the
system of linear equations In mathematics, a system of linear equations (or linear system) is a collection of one or more linear equations involving the same variables. For example, :\begin 3x+2y-z=1\\ 2x-2y+4z=-2\\ -x+\fracy-z=0 \end is a system of three equations in t ...
:\mathbf\mathbf = \mathbf for the vector \mathbf, where the known n \times n matrix \mathbf is
symmetric Symmetry (from grc, συμμετρία "agreement in dimensions, due proportion, arrangement") in everyday language refers to a sense of harmonious and beautiful proportion and balance. In mathematics, "symmetry" has a more precise definit ...
(i.e., AT = A),
positive-definite In mathematics, positive definiteness is a property of any object to which a bilinear form or a sesquilinear form may be naturally associated, which is positive-definite. See, in particular: * Positive-definite bilinear form * Positive-definite ...
(i.e. xTAx > 0 for all non-zero vectors \mathbf in R''n''), and
real Real may refer to: Currencies * Brazilian real (R$) * Central American Republic real * Mexican real * Portuguese real * Spanish real * Spanish colonial real Music Albums * ''Real'' (L'Arc-en-Ciel album) (2000) * ''Real'' (Bright album) (2010) ...
, and \mathbf is known as well. We denote the unique solution of this system by \mathbf_*.


Derivation as a direct 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/ Lanczos iteration for
eigenvalue In linear algebra, an eigenvector () or characteristic vector of a linear transformation is a nonzero vector that changes at most by a scalar factor when that linear transformation is applied to it. The corresponding eigenvalue, often denoted ...
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. We say that two non-zero vectors u and v are conjugate (with respect to \mathbf) if : \mathbf^\mathsf \mathbf \mathbf = 0. Since \mathbf is symmetric and positive-definite, the left-hand side defines an
inner product In mathematics, an inner product space (or, rarely, a Hausdorff pre-Hilbert space) is a real vector space or a complex vector space with an operation called an inner product. The inner product of two vectors in the space is a scalar, often ...
: \mathbf^\mathsf \mathbf \mathbf = \langle \mathbf, \mathbf \rangle_\mathbf := \langle \mathbf \mathbf, \mathbf\rangle = \langle \mathbf, \mathbf^\mathsf \mathbf\rangle = \langle \mathbf, \mathbf\mathbf\rangle. Two vectors are conjugate if and only if they are orthogonal with respect to this inner product. Being conjugate is a symmetric relation: if \mathbf is conjugate to \mathbf, then \mathbf is conjugate to \mathbf. Suppose that :P = \ is a set of n mutually conjugate vectors with respect to \mathbf, i.e. \mathbf_i^\mathsf \mathbf \mathbf_j = 0 for all i \neq j. Then P forms a
basis Basis may refer to: Finance and accounting *Adjusted basis, the net cost of an asset after adjusting for various tax-related items *Basis point, 0.01%, often used in the context of interest rates *Basis trading, a trading strategy consisting of ...
for \mathbb^n, and we may express the solution \mathbf_* of \mathbf = \mathbf in this basis: :\mathbf_* = \sum^_ \alpha_i \mathbf_i \Rightarrow \mathbf \mathbf_* = \sum^_ \alpha_i \mathbf \mathbf_i. Left-multiplying by \mathbf_k^\mathsf yields : \mathbf_k^\mathsf \mathbf = \mathbf_k^\mathsf \mathbf \mathbf_* = \sum^_ \alpha_i \mathbf_k^\mathsf \mathbf \mathbf_i = \sum^_ \alpha_i \left \langle \mathbf_k, \mathbf_i \right \rangle_ = \alpha_k \left \langle \mathbf_k, \mathbf_k \right \rangle_ and so :\alpha_k = \frac. This gives the following method for solving the equation : find a sequence of n conjugate directions, and then compute the coefficients \alpha_k.


As an iterative method

If we choose the conjugate vectors \mathbf_k carefully, then we may not need all of them to obtain a good approximation to the solution \mathbf_*. 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 by (we can assume without loss of generality that , otherwise consider the system Az = b − Ax0 instead). Starting with x0 we search for the solution and in each iteration we need a metric to tell us whether we are closer to the solution (that is unknown to us). This metric comes from the fact that the solution is also the unique minimizer of the following
quadratic function In mathematics, a quadratic polynomial is a polynomial of degree two in one or more variables. A quadratic function is the polynomial function defined by a quadratic polynomial. Before 20th century, the distinction was unclear between a polynomia ...
: f(\mathbf) = \tfrac12 \mathbf^\mathsf \mathbf\mathbf - \mathbf^\mathsf \mathbf, \qquad \mathbf\in\mathbf^n \,. The existence of a unique minimizer is apparent as its
Hessian matrix In mathematics, the Hessian matrix or Hessian is a square matrix of second-order partial derivatives of a scalar-valued function, or scalar field. It describes the local curvature of a function of many variables. The Hessian matrix was developed ...
of second derivatives is symmetric positive-definite : \mathbf(f(\mathbf)) = \mathbf \,, and that the minimizer (use D''f''(x)=0) solves the initial problem is obvious from its first derivative : \nabla f(\mathbf) = \mathbf \mathbf - \mathbf \,. 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 . Starting with an initial guess x0, 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''. Note that p0 is also the residual provided by this initial step of the algorithm. Let r''k'' be the residual at the ''k''th step: : \mathbf_k = \mathbf - \mathbf_k. As observed above, \mathbf_k is the negative gradient of f at \mathbf_k, so the
gradient descent In mathematics, gradient descent (also often called steepest descent) is a first-order iterative optimization algorithm for finding a local minimum of a differentiable function. The idea is to take repeated steps in the opposite direction of ...
method would require to move in the direction r''k''. Here, however, we insist that the directions \mathbf_k must be conjugate to each other. A practical way to enforce this is by requiring that the next search direction be built out of the current residual and all previous search directions. The conjugation constraint is an orthonormal-type constraint and hence the algorithm can be viewed as an example of Gram-Schmidt orthonormalization. This gives the following expression: :\mathbf_ = \mathbf_ - \sum_\frac \mathbf_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 : \mathbf_ = \mathbf_k + \alpha_k \mathbf_k with : \alpha_ = \frac = \frac, where the last equality follows from the definition of \mathbf_k . The expression for \alpha_k can be derived if one substitutes the expression for x''k''+1 into ''f'' and minimizing it w.r.t. \alpha_k : \begin f(\mathbf_) &= f(\mathbf_k + \alpha_k \mathbf_k) =: g(\alpha_k) \\ g'(\alpha_k) &\overset 0 \quad \Rightarrow \quad \alpha_ = \frac \,. \end


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 \mathbf_i is orthogonal to \mathbf_j, i.e. \mathbf_i^\mathsf \mathbf_j=0 , for i ≠ j. And \mathbf_iis \mathbf-orthogonal to \mathbf_j , i.e. \mathbf_i^\mathsf \mathbf \mathbf_j=0 , for i \neq j. This can be regarded that as the algorithm progresses, \mathbf_i and \mathbf_i span the same
Krylov subspace In linear algebra, the order-''r'' Krylov subspace generated by an ''n''-by-''n'' matrix ''A'' and a vector ''b'' of dimension ''n'' is the linear subspace spanned by the images of ''b'' under the first ''r'' powers of ''A'' (starting from A^0=I ...
. Where \mathbf_i form the orthogonal basis with respect to the standard inner product, and \mathbf_i form the orthogonal basis with respect to the inner product induced by \mathbf. Therefore, \mathbf_k can be regarded as the projection of \mathbf on the Krylov subspace. The algorithm is detailed below for solving Ax = b where \mathbf is a real, symmetric, positive-definite matrix. The input vector \mathbf_0 can be an approximate initial solution or 0. It is a different formulation of the exact procedure described above. :\begin & \mathbf_0 := \mathbf - \mathbf_0 \\ & \hbox \mathbf_ \text \mathbf_ \text\\ & \mathbf_0 := \mathbf_0 \\ & k := 0 \\ & \text \\ & \qquad \alpha_k := \frac \\ & \qquad \mathbf_ := \mathbf_k + \alpha_k \mathbf_k \\ & \qquad \mathbf_ := \mathbf_k - \alpha_k \mathbf_k \\ & \qquad \hbox \mathbf_ \text \\ & \qquad \beta_k := \frac \\ & \qquad \mathbf_ := \mathbf_ + \beta_k \mathbf_k \\ & \qquad k := k + 1 \\ & \text \\ & \text \mathbf_ \text \end This is the most commonly used algorithm. The same formula for is also used in the Fletcher–Reeves nonlinear conjugate gradient method.


Restarts

We note that \mathbf_ is computed by the
gradient descent In mathematics, gradient descent (also often called steepest descent) is a first-order iterative optimization algorithm for finding a local minimum of a differentiable function. The idea is to take repeated steps in the opposite direction of ...
method applied to \mathbf_. Setting \beta_=0 would similarly make \mathbf_ computed by the
gradient descent In mathematics, gradient descent (also often called steepest descent) is a first-order iterative optimization algorithm for finding a local minimum of a differentiable function. The idea is to take repeated steps in the opposite direction of ...
method from \mathbf_, i.e., can be used as a simple implementation of a restart of the conjugate gradient iterations. Restarts could slow down convergence, but may improve stability if the conjugate gradient method misbehaves, e.g., due to
round-off error A roundoff error, also called rounding error, is the difference between the result produced by a given algorithm using exact arithmetic and the result produced by the same algorithm using finite-precision, rounded arithmetic. Rounding errors are d ...
.


Explicit residual calculation

The formulas \mathbf_ := \mathbf_k + \alpha_k \mathbf_k and \mathbf_k := \mathbf - \mathbf_k, which both hold in exact arithmetic, make the formulas \mathbf_ := \mathbf_k - \alpha_k \mathbf_k and \mathbf_ := \mathbf - \mathbf_ mathematically equivalent. The former is used in the algorithm to avoid an extra multiplication by \mathbf since the vector \mathbf_k is already computed to evaluate \alpha_k. The latter may be more accurate, substituting the explicit calculation \mathbf_ := \mathbf - \mathbf_ for the implicit one by the recursion subject to
round-off error A roundoff error, also called rounding error, is the difference between the result produced by a given algorithm using exact arithmetic and the result produced by the same algorithm using finite-precision, rounded arithmetic. Rounding errors are d ...
accumulation, and is thus recommended for an occasional evaluation. A norm of the residual is typically used for stopping criteria. The norm of the explicit residual \mathbf_ := \mathbf - \mathbf_ provides a guaranteed level of accuracy both in exact arithmetic and in the presence of the rounding errors, where convergence naturally stagnates. In contrast, the implicit residual \mathbf_ := \mathbf_k - \alpha_k \mathbf_k is known to keep getting smaller in amplitude well below the level of rounding errors and thus cannot be used to determine the stagnation of convergence.


Computation of alpha and beta

In the algorithm, is chosen such that \mathbf_ is orthogonal to \mathbf_. The denominator is simplified from :\alpha_k = \frac = \frac since \mathbf_ = \mathbf_-\mathbf_\mathbf_. The is chosen such that \mathbf_ is conjugate to \mathbf_. Initially, is :\beta_k = - \frac using :\mathbf_ = \mathbf_ - \alpha_ \mathbf \mathbf_ and equivalently \mathbf \mathbf_ = \frac (\mathbf_ - \mathbf_), the numerator of is rewritten as : \mathbf_^\mathsf \mathbf \mathbf_k = \frac \mathbf_^\mathsf (\mathbf_k - \mathbf_) = - \frac \mathbf_^\mathsf \mathbf_ because \mathbf_ and \mathbf_ are orthogonal by design. The denominator is rewritten as : \mathbf_k^\mathsf \mathbf \mathbf_k = (\mathbf_k + \beta_ \mathbf_)^\mathsf \mathbf \mathbf_k = \frac \mathbf_k^\mathsf (\mathbf_k - \mathbf_) = \frac \mathbf_k^\mathsf \mathbf_k using that the search directions p''k'' are conjugated and again that the residuals are orthogonal. This gives the in the algorithm after cancelling .


Example code in

MATLAB MATLAB (an abbreviation of "MATrix LABoratory") is a proprietary multi-paradigm programming language and numeric computing environment developed by MathWorks. MATLAB allows matrix manipulations, plotting of functions and data, implementatio ...
/
GNU Octave GNU Octave is a high-level programming language primarily intended for scientific computing and numerical computation. Octave helps in solving linear and nonlinear problems numerically, and for performing other numerical experiments using a langu ...

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

Consider the linear system Ax = b given by :\mathbf \mathbf= \begin 4 & 1 \\ 1 & 3 \end \begin x_1 \\ x_2 \end = \begin 1 \\ 2 \end, we will perform two steps of the conjugate gradient method beginning with the initial guess :\mathbf_0 = \begin 2 \\ 1 \end in order to find an approximate solution to the system.


Solution

For reference, the exact solution is : \mathbf = \begin \frac \\\\ \frac \end \approx \begin 0.0909 \\\\ 0.6364 \end 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 :\mathbf_0 = \begin 1 \\ 2 \end - \begin 4 & 1 \\ 1 & 3 \end \begin 2 \\ 1 \end = \begin-8 \\ -3 \end = \mathbf_0. Since this is the first iteration, we will use the residual vector r0 as our initial search direction p0; the method of selecting p''k'' will change in further iterations. We now compute the scalar using the relationship : \alpha_0 = \frac = \frac =\frac\approx0.2205 We can now compute x1 using the formula :\mathbf_1 = \mathbf_0 + \alpha_0\mathbf_0 = \begin 2 \\ 1 \end + \frac \begin -8 \\ -3 \end \approx \begin 0.2356 \\ 0.3384 \end. 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 :\mathbf_1 = \mathbf_0 - \alpha_0 \mathbf \mathbf_0 = \begin -8 \\ -3 \end - \frac \begin 4 & 1 \\ 1 & 3 \end \begin -8 \\ -3 \end \approx \begin -0.2810 \\ 0.7492 \end. Our next step in the process is to compute the scalar that will eventually be used to determine the next search direction p1. :\beta_0 = \frac \approx \frac = 0.0088. Now, using this scalar , we can compute the next search direction p1 using the relationship :\mathbf_1 = \mathbf_1 + \beta_0 \mathbf_0 \approx \begin -0.2810 \\ 0.7492 \end + 0.0088 \begin -8 \\ -3 \end = \begin -0.3511 \\ 0.7229 \end. We now compute the scalar using our newly acquired p1 using the same method as that used for . : \alpha_1 = \frac \approx \frac = 0.4122. Finally, we find x2 using the same method as that used to find x1. :\mathbf_2 = \mathbf_1 + \alpha_1 \mathbf_1 \approx \begin 0.2356 \\ 0.3384 \end + 0.4122 \begin -0.3511 \\ 0.7229 \end = \begin 0.0909 \\ 0.6364 \end. 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

The conjugate gradient method can theoretically be viewed as a direct method, as in the absence of
round-off error A roundoff error, also called rounding error, is the difference between the result produced by a given algorithm using exact arithmetic and the result produced by the same algorithm using finite-precision, rounded arithmetic. Rounding errors are d ...
it produces the exact solution after a finite number of iterations, which is not larger than the size of the matrix. In practice, the exact solution is never obtained since the conjugate gradient method is unstable with respect to even small perturbations, e.g., most directions are not in practice conjugate, due to a degenerative nature of generating the Krylov subspaces. As an
iterative method In computational mathematics, an iterative method is a mathematical procedure that uses an initial value to generate a sequence of improving approximate solutions for a class of problems, in which the ''n''-th approximation is derived from the pre ...
, the conjugate gradient method monotonically (in the energy norm) improves approximations \mathbf_ to the exact solution and 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 In numerical analysis, the condition number of a function measures how much the output value of the function can change for a small change in the input argument. This is used to measure how sensitive a function is to changes or errors in the inpu ...
\kappa(A) of the system matrix A: the larger \kappa(A) is, the slower the improvement. If \kappa(A) is large,
preconditioning In mathematics, preconditioning is the application of a transformation, called the preconditioner, that conditions a given problem into a form that is more suitable for numerical solving methods. Preconditioning is typically related to reducing ...
is commonly used to replace the original system \mathbf-\mathbf = 0 with \mathbf^(\mathbf-\mathbf) = 0 such that \kappa(\mathbf^\mathbf) is smaller than \kappa(\mathbf), see below.


Convergence theorem

Define a subset of polynomials as : \Pi_k^* := \left\lbrace \ p \in \Pi_k \ : \ p(0)=1 \ \right\rbrace \,, where \Pi_k is the set of
polynomials In mathematics, a polynomial is an expression consisting of indeterminates (also called variables) and coefficients, that involves only the operations of addition, subtraction, multiplication, and positive-integer powers of variables. An exampl ...
of maximal degree k . Let \left( \mathbf_k \right)_k be the iterative approximations of the exact solution \mathbf_* , and define the errors as \mathbf_k := \mathbf_k - \mathbf_* . Now, the rate of convergence can be approximated as : \begin \left\, \mathbf_k \right\, _\mathbf &= \min_ \left\, p(\mathbf) \mathbf_0 \right\, _\mathbf \\ &\leq \min_ \, \max_ , p(\lambda) , \ \left\, \mathbf_0 \right\, _\mathbf \\ &\leq 2 \left( \frac \right)^k \ \left\, \mathbf_0 \right\, _\mathbf \,, \end where \sigma(\mathbf) denotes the
spectrum A spectrum (plural ''spectra'' or ''spectrums'') is a condition that is not limited to a specific set of values but can vary, without gaps, across a continuum. The word was first used scientifically in optics to describe the rainbow of colors ...
, and \kappa(\mathbf) denotes the
condition number In numerical analysis, the condition number of a function measures how much the output value of the function can change for a small change in the input argument. This is used to measure how sensitive a function is to changes or errors in the inpu ...
. Note, the important limit when \kappa(\mathbf) tends to \infty : \frac \approx 1 - \frac \quad \text \quad \kappa(\mathbf) \gg 1 \,. This limit shows a faster convergence rate compared to the iterative methods of Jacobi or Gauss–Seidel which scale as \approx 1 - \frac . No
round-off error A roundoff error, also called rounding error, is the difference between the result produced by a given algorithm using exact arithmetic and the result produced by the same algorithm using finite-precision, rounded arithmetic. Rounding errors are d ...
is assumed in the convergence theorem, but the convergence bound is commonly valid in practice as theoretically explained by Anne Greenbaum.


Practical convergence

If initialized randomly, the first stage of iterations is often the fastest, as the error is eliminated within the Krylov subspace that initially reflects a smaller effective condition number. The second stage of convergence is typically well defined by the theoretical convergence bound with \sqrt, but may be super-linear, depending on a distribution of the spectrum of the matrix A and the spectral distribution of the error. In the last stage, the smallest attainable accuracy is reached and the convergence stalls or the method may even start diverging. In typical scientific computing applications in
double-precision floating-point format Double-precision floating-point format (sometimes called FP64 or float64) is a floating-point number format, usually occupying 64 bits in computer memory; it represents a wide dynamic range of numeric values by using a floating radix point. Flo ...
for matrices of large sizes, the conjugate gradient method uses a stopping criteria with a tolerance that terminates the iterations during the first or second stage.


The preconditioned conjugate gradient method

In most cases,
preconditioning In mathematics, preconditioning is the application of a transformation, called the preconditioner, that conditions a given problem into a form that is more suitable for numerical solving methods. Preconditioning is typically related to reducing ...
is necessary to ensure fast convergence of the conjugate gradient method. If \mathbf^ is symmetric positive-definite and \mathbf^\mathbf has a better condition number \mathbf, a preconditioned conjugate gradient method can be used. It takes the following form: :\mathbf_0 := \mathbf - \mathbf_0 :\mathbf_0 := \mathbf^ \mathbf_0 :\mathbf_0 := \mathbf_0 :k := 0 \, :repeat ::\alpha_k := \frac ::\mathbf_ := \mathbf_k + \alpha_k \mathbf_k ::\mathbf_ := \mathbf_k - \alpha_k \mathbf_k ::if r''k''+1 is sufficiently small then exit loop end if ::\mathbf_ := \mathbf^ \mathbf_ ::\beta_k := \frac ::\mathbf_ := \mathbf_ + \beta_k \mathbf_k ::k := k + 1 \, :end repeat :The result is x''k''+1 The above formulation is equivalent to applying the regular conjugate gradient method to the preconditioned system :\mathbf^\mathbf(\mathbf^)^\mathsf\mathbf=\mathbf^\mathbf where :\mathbf^\mathsf=\mathbf, \qquad \mathbf=\mathbf^\mathsf\mathbf. The Cholesky decomposition of the preconditioner must be used to keep the symmetry (and positive definiteness) of the system. However, this decomposition does not need to be computed, and it is sufficient to know \mathbf^. It can be shown that \mathbf^\mathbf(\mathbf^)^\mathsf has the same spectrum as \mathbf^\mathbf. 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 In mathematics, preconditioning is the application of a transformation, called the preconditioner, that conditions a given problem into a form that is more suitable for numerical solving methods. Preconditioning is typically related to reducing ...
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 Polak–Ribière formula :\beta_k := \frac instead of the Fletcher–Reeves formula :\beta_k := \frac may dramatically improve the convergence in this case. This version of the preconditioned conjugate gradient method can be called flexible, as it allows for variable preconditioning. The flexible version is also shown to be robust even if the preconditioner is not symmetric positive definite (SPD). The implementation of the flexible version requires storing an extra vector. For a fixed SPD preconditioner, \mathbf_^\mathsf \mathbf_=0, so both formulas for are equivalent in exact arithmetic, i.e., without the
round-off error A roundoff error, also called rounding error, is the difference between the result produced by a given algorithm using exact arithmetic and the result produced by the same algorithm using finite-precision, rounded arithmetic. Rounding errors are d ...
. 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.


Vs. the locally optimal steepest descent method

In both the original and the preconditioned conjugate gradient methods one only needs to set \beta_k := 0 in order to make them locally optimal, using the
line search In optimization, the line search strategy is one of two basic iterative approaches to find a local minimum \mathbf^* of an objective function f:\mathbb R^n\to\mathbb R. The other approach is trust region. The line search approach first finds ...
,
steepest descent In mathematics, gradient descent (also often called steepest descent) is a first-order iterative optimization algorithm for finding a local minimum of a differentiable function. The idea is to take repeated steps in the opposite direction of ...
methods. With this substitution, vectors are always the same as vectors , so there is no need to store vectors . Thus, every iteration of these
steepest descent In mathematics, gradient descent (also often called steepest descent) is a first-order iterative optimization algorithm for finding a local minimum of a differentiable function. The idea is to take repeated steps in the opposite direction of ...
methods is a bit cheaper compared to that for the conjugate gradient methods. However, the latter converge faster, unless a (highly) variable and/or non-SPD
preconditioner In mathematics, preconditioning is the application of a transformation, called the preconditioner, that conditions a given problem into a form that is more suitable for numerical solving methods. Preconditioning is typically related to reducing ...
is used, see above.


Conjugate gradient method as optimal feedback controller for double integrator

The conjugate gradient method can also be derived using
optimal control theory Mathematical optimization (alternatively spelled ''optimisation'') or mathematical programming is the selection of a best element, with regard to some criterion, from some set of available alternatives. It is generally divided into two subfi ...
. Ross, I. M., "An Optimal Control Theory for Accelerated Optimization," , 2019. In this approach, the conjugate gradient method falls out as an optimal feedback controller,u = k(x, v):= -\gamma_a \nabla f(x) - \gamma_b v for the double integrator system,\dot x = v, \quad \dot v = u The quantities \gamma_a and \gamma_b are variable feedback gains.


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 In statistics, ordinary least squares (OLS) is a type of linear least squares method for choosing the unknown parameters in a linear regression model (with fixed level-one effects of a linear function of a set of explanatory variables) by the prin ...
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). : ATAx = ATb 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 In numerical analysis and scientific computing, a sparse matrix or sparse array is a matrix in which most of the elements are zero. There is no strict definition regarding the proportion of zero-value elements for a matrix to qualify as sparse b ...
since these operations are usually extremely efficient. However the downside of forming the normal equations is that the
condition number In numerical analysis, the condition number of a function measures how much the output value of the function can change for a small change in the input argument. This is used to measure how sensitive a function is to changes or errors in the inpu ...
κ(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 In mathematics, preconditioning is the application of a transformation, called the preconditioner, that conditions a given problem into a form that is more suitable for numerical solving methods. Preconditioning is typically related to reducing ...
is often an important part of using the CGNR method. Several algorithms have been proposed (e.g., CGLS, LSQR). Th
LSQR
algorithm purportedly has the best numerical stability when A is ill-conditioned, i.e., A has a large
condition number In numerical analysis, the condition number of a function measures how much the output value of the function can change for a small change in the input argument. This is used to measure how sensitive a function is to changes or errors in the inpu ...
.


Conjugate gradient method for complex Hermitian matrices

The conjugate gradient method with a trivial modification is extendable to solving, given complex-valued matrix A and vector b, the system of linear equations \mathbf \mathbf =\mathbf for the complex-valued vector x, where A is
Hermitian {{Short description, none Numerous things are named after the French mathematician Charles Hermite (1822–1901): Hermite * Cubic Hermite spline, a type of third-degree spline * Gauss–Hermite quadrature, an extension of Gaussian quadrature me ...
(i.e., A' = A) and
positive-definite matrix In mathematics, a symmetric matrix M with real entries is positive-definite if the real number z^\textsfMz is positive for every nonzero real column vector z, where z^\textsf is the transpose of More generally, a Hermitian matrix (that is, a ...
, and the symbol ' denotes the
conjugate transpose In mathematics, the conjugate transpose, also known as the Hermitian transpose, of an m \times n complex matrix \boldsymbol is an n \times m matrix obtained by transposing \boldsymbol and applying complex conjugate on each entry (the complex c ...
using the
MATLAB MATLAB (an abbreviation of "MATrix LABoratory") is a proprietary multi-paradigm programming language and numeric computing environment developed by MathWorks. MATLAB allows matrix manipulations, plotting of functions and data, implementatio ...
/
GNU Octave GNU Octave is a high-level programming language primarily intended for scientific computing and numerical computation. Octave helps in solving linear and nonlinear problems numerically, and for performing other numerical experiments using a langu ...
style. The trivial modification is simply substituting the
conjugate transpose In mathematics, the conjugate transpose, also known as the Hermitian transpose, of an m \times n complex matrix \boldsymbol is an n \times m matrix obtained by transposing \boldsymbol and applying complex conjugate on each entry (the complex c ...
for the real
transpose In linear algebra, the transpose of a matrix is an operator which flips a matrix over its diagonal; that is, it switches the row and column indices of the matrix by producing another matrix, often denoted by (among other notations). The tr ...
everywhere. This substitution is backward compatible, since
conjugate transpose In mathematics, the conjugate transpose, also known as the Hermitian transpose, of an m \times n complex matrix \boldsymbol is an n \times m matrix obtained by transposing \boldsymbol and applying complex conjugate on each entry (the complex c ...
turns into real
transpose In linear algebra, the transpose of a matrix is an operator which flips a matrix over its diagonal; that is, it switches the row and column indices of the matrix by producing another matrix, often denoted by (among other notations). The tr ...
on real-valued vectors and matrices. The provided above Example code in MATLAB/GNU Octave thus already works for complex Hermitian matrices needed no modification.


See also

* Biconjugate gradient method (BiCG) *
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 pr ...
* Gaussian belief propagation * Iterative method: Linear systems *
Krylov subspace In linear algebra, the order-''r'' Krylov subspace generated by an ''n''-by-''n'' matrix ''A'' and a vector ''b'' of dimension ''n'' is the linear subspace spanned by the images of ''b'' under the first ''r'' powers of ''A'' (starting from A^0=I ...
* Nonlinear conjugate gradient method *
Preconditioning In mathematics, preconditioning is the application of a transformation, called the preconditioner, that conditions a given problem into a form that is more suitable for numerical solving methods. Preconditioning is typically related to reducing ...
* Sparse matrix–vector multiplication


References


Further reading

* * * *


External links

* {{DEFAULTSORT:Conjugate Gradient Method Numerical linear algebra Gradient methods Articles with example MATLAB/Octave code