http://bbs.sciencenet.cn/blog-588243-573819.html
共轭梯度法（Conjugate gradient method）
（cited from:wiki）
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. （看到了吧，系数矩阵A必须是对称的和正定的）The conjugate gradient method is an iterative method, so it can be applied to sparse systems that are too large to be handled by direct methods such as the Cholesky decomposition. Such systems often arise when numerically solving partial differential equations.
----------------------
正定矩阵：
在线性代数里，正定矩阵（即“正数-确定-矩阵”）是埃尔米特矩阵的一种，有时会简称为正定阵。
一个n × n的实对称矩阵 M 是正定的当且仅当对于所有的非零实系数向量z，都有 zTMz > 0。其中zT 表示z的转置。
对于复数的情况，定义则为：一个n × n的埃尔米特矩阵 M 是正定的当且仅当对于每个非零的复向量z，都有z*Mz > 0。其中z* 表示z的共轭转置。由于 M 是埃尔米特矩阵，经计算可知，对于任意的复向量z，z*Mz必然是实数，从而可以与0比较大小。因此这个定义是自洽的。
---------------------
The conjugate gradient method can also be used to solve unconstrained optimization problems such as energy minimization.
The biconjugate gradient method provides a generalization to non-symmetric matrices. Various nonlinear conjugate gradient methods seek minima of nonlinear equations.
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).
Description of the method
Suppose we want to solve the following system of linear equations
Ax = b
where the n-by-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.
We denote the unique solution of this system by x*.
The conjugate gradient method as a direct method
We say that two non-zero vectors u and v are conjugate (with respect to A) if
Since A is symmetric and positive definite, the left-hand side defines an inner product
-
-
<Au,v>=(Au)T*v=(uT)*(AT)*v=<u,ATv>
So, 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. (Note: This notion of conjugate is not related to the notion of complex conjugate.)
Suppose that {pk} is a sequence of n mutually conjugate directions. Then the pk form a basis of Rn, so we can expand the solution x* of Ax = b in this basis:
The coefficients are given by
-
-
(because
are mutually conjugate)
-
This gives the following method for solving the equation Ax = b. We first find a sequence of n conjugate directions and then we compute the coefficients αk.
The conjugate gradient method as an iterative method
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 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 x0 we 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; so if f(x) becomes smaller in an iteration it means that we are closer to x*.
This suggests taking the first basis vector p1 to be the negative of the gradient of f at x = x0. This gradient equals Ax0−b. Since x0 = 0, this means we take p1 = b. 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:
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 the next search direction is 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:
(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:
with
--------------------
在线性代数中，如果内积空间上的一组矢量能够张成一个子空间，那么这一组矢量就称为这个子空间的一个基。
Gram－Schmidt正交化提供了一种方法，能够通过这一子空间上的一个基得出子空间的一个正交基，并可进一步求出对应的标准正交基。
--------------------
The resulting algorithm
The above algorithm gives the most straightforward explanation towards the conjugate gradient method. However, it requires storage of all the previous searching directions and residue vectors, and many matrix vector multiplications, thus could be computationally expensive. In practice, one modifies slightly the condition obtaining the last residue vector, not to minimize the metric following the search direction, but instead to make it orthogonal to the previous residue. Minimization of the metric along the search direction will be obtained automatically in this case. One can then result in an algorithm which only requires storage of the last two residue vectors and the last search direction, and only one matrix vector multiplication.Note that the algorithm described below is equivalent to the previously discussed straightforward procedure.
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.
-
-
-
-
repeat
-
-
-
- ifr k+1 is sufficiently small then exit loop end if
-
-
-
-
- end repeat
-
The result is
x
k+1
This is the most commonly used algorithm.
Numerical example
To illustrate the conjugate gradient method, we will complete a simple example.
Considering the linear system Ax = b given by
we will perform two steps of the conjugate gradient method beginning with the initial guess
in order to find an approximate solution to the system.
Solution
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
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
We can now compute x1 using the formula
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
Our next step in the process is to compute the scalar β0 that will eventually be used to determine the next search direction p1.
Now, using this scalar β0, we can compute the next search direction p1 using the relationship
We now compute the scalar α1 using our newly-acquired p1 using the same method as that used for α0.
Finally, we find x2 using the same method as that used to find x1.
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 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 
If κ(A) is large, preconditioning is used to replace the original system 
按照上面的算法，其MATLAB程序，如下：（需要的人可自行拷取，如有错误，敬请告知！）
% % % %-------------------------------------------------------------------------------------
function x = myconjgrad(A,b,tol)
% MYCONJGRAD my Conjugate Gradient Method.
% x = MYCONJGRAD(A,b,tol) attemps to solve the system of linear equations A*x=b
% for x. The N-by-N coefficient matrix A must be symmetric,positive definite and the right
% hand side column vector B must have the same length as A.
% input parameter tol specifies the tolerance of the method.
%
% Example:
%{
n = 600;
m = 800;
A = randn(n,m);
A = A * A';
b = randn(n,1);
tic, x = myconjgrad(A,b,1e-4); toc
norm(A*x-b)
%}
%
x = b; % assuming x0 is b
r = b - A*x; % the residual vector r0
p = r; % the first search direction
for k = 1:numel(b);
z = A*p;
alpha = (r'*r)/(p'*z); % the scalar
x = x + alpha*p; % x1 = x0 + alpha0*p0
s = r'*r;
r = r - alpha*z; % the residual vector r
if( norm(r) < tol )
return;
end
B = (r'*r)/s;
p = r + B*p; % y is the search direction p
end
% % % %-------------------------------------------------------------------------------------
从程序中可看到，关键的变量有系数αk,由之可得新的计算点x位置，新的残差向量rk，以及β，由它可以求得新的搜索方向pk。
程序编写时，参考了该地址的程序。两者结果不一样，在所实验的条件下，本文结果较优。但是不明白参考版本的计算思路。或许有其他优点？有知道的请告知。