An Introduction to the Conjugate Gradient Method without the Agonizing Pain Jonathan Richard Shewchuk Reading Group Presention By David Cline
Linear System Unknown vector (what we want to find) Known vector Square matrix
Positive Definite Matrix [ x1 x2 … xn] > 0 * Also, all eigenvalues of the matrix are positive
Quadtratic form • An expression of the form
Why do we care? The gradient of the quadratic form is our original system if A is symmetric:
Visual representation f(x) f’(x) f(x)
Solution • the solution to the system, x, is the global minimum of f. … if A is symmetric, • And since A is positive definite, x is the global minimum of f
Definitions • Error • Residual Whenever you read ‘residual’, Think ‘the direction of steepest Descent’.
Method of steepest descent • Start with arbitrary point, x(0) • move in direction opposite gradient of f, r(0) • reach minimum in that direction at distance alpha • repeat
Steepest descent does well: Steepest descent converges in one Iteration if the error term is an Eigenvector. Steepest descent converges in one Iteration if the all the eigenvalues Are equal.
Steepest descent does poorly If the error term is a mix of large and small eigenvectors, steepest descent will move back and forth along toward the solution, but take many iterations to converge. The worst case convergence is related to the ratio of the largest and smallest eigenvalues of A, called the “condition number”:
Convergence of steepest descent: # iterations “energy norm” at iteration i “energy norm” at iteration 0
How can we speed up or guarantee convergence? • Use the eigenvectors as directions. • terminates in n iterations.
Method of conjugate directions • Instead of eigenvectors, which are too hard to compute, use directions that are “conjugate” or “A-orthogonal”:
How to find conjugate directions? • Gram-Shmidt Conjugation: Start with n linearly independent vectors u0…un-1 • For each vector, subract those parts that are not A-orthogonal to the other processed vectors:
Problem • Gram-Schmidt conjugation is slow and we have to store all of the vectors we have created.
Conjugate Gradient Method • Apply the method of conjugate directions, but use the residuals for the u values: ui = r(i)
How does this help us? • It turns out that the residual ri is A-orthogonal to all of the previous residuals, except ri-1, so we simply make it A-orthogonal to ri-1, and we are set.
Simplifying further k=i-1
Putting it all together Start with steepest descent Compute distance to bottom Of parabola Slide down to bottom of parabola Compute steepest descent At next location Remove part of vector that Is not A-orthogonal to di
Starting and stopping • Start either with a rough estimate of the solution, or the zero vector. • Stop when the norm of the residual is small enough.
Diagonal preconditioning • Just use the diagonal of A as M. A diagonal matrix is easy to invert, but of course it isn’t the best method out there.
CG on the normal equations If A is not symmetric, or positive-definite, or not square, we can’t use CG directly to solve However, we can use it to solve the system is always symmetric, positive definite and square. The problem that we solve with this is the least-squares fit but the condition number increases. Also note that we never actually have to form Instead we multiply by AT and then by A.