I'm having problems understanding why the conjugate gradient method breaks down for singular matrices. I've read a good introduction to intuitively understanding the CG method through visualizing the quadratic form, but I'm not really understanding mathematically why this seems to be. Does anyone have a good mathematical explanation or any sources that explain why conjugate gradient iteration doesn't work for singular matrices?
Thanks in advance.
It's important to distinguish between CG as an exact method, as originally conceived by Hestenes and Stiefel and as an iterative Krylov-type method. As an iterative method, CG works fine for singular matrices provided you don't do too many iterations. For ill-conditioned systems loss of A-orthogonality of the search vectors causes the method to slow down (and produce spurious eigenvalues if you're calculating them, rather than a solution). Also, as Hestenes later pointed out (1975), CG will converge to the pseudo-inverse solution if you start with an initial approximation of zero.
Like any Krylov subspace method CG will compute the components of the solution starting with the largest eigenvalues. In some problems the largest few eigenvalues are sufficient to get a reasonably accurate solution. In this case CG is great. If you really need the small eigenvalue contribution, then you need to re-orthogonalize the search vectors, or partially re-orthogonalize them. This is not hard. In any case actual singularity is no real problem for the iterative method. Stability issues were solved many years ago.
You can even use CG to compute an SVD of a sparse rectangular matrix. Look at SVDPACK by Michael Berry.
Let my reply here to your comment:
My question was geared more towards any mathematical explanation of why singular matrices don't work. My understanding, from the previously referenced article (pg 5), is that singular matrices have multiple solutions rather than a single solution.
Unfortunately, that's not so. If $A$ is a singular matrix, then for most $b$, the equation $Ax=b$ has no solution. It has a solution if and only if $b$ is in the column space of $A$ (or orthogonal to the kernel of $A$, if $A$ is symmetric), in which case there are infinitely many solutions. This is basic linear algebra.
Let's say $A$ is symmetric positive semidefinite, so we can at least entertain the thought of conjugate gradients. Conjugate gradients can be interpreted as maximizing the quadratic function $$f(x) = b^Tx-\frac12 x^TAx.$$ If $A$ is singular, then this function typically has no maximum value, as I will demonstrate. Since $A$ is singular, there are vectors $v$ such that $Av = 0$. If there is any such vector $v$ that is not orthogonal to $b$, consider substituting $x=tv$ into $f(x)$: $$f(tv) = t(b^Tv) - \frac12t^2(v^TAv) = t(b^Tv) = ct,$$ where $c$ is the nonzero real number $b^Tv$. Obviously, this is unbounded as a function of $t$, so we see that there is no maximum value if $b$ is not orthogonal to the kernel of $A$.
A search for "conjugate gradient singular matrix" took me to this question. While the answer is obviously given by the responses, the question can be refined: Can CG still give a working algorithm if the matrix is singular, but behaves as a symmetric positive definite form on a (large) subspace?
A standard example is given by the finite element discretization of the Neumann problem on a simply connected domain. The constant functions are both the kernel and the cokernel of the Laplacian. On functions with vanishing mean, the Laplacian is still a positive definite symmetric operator, and we would like to leverage this structure.
This is non-trivial and best our numerical method is derived from a fully analytic setting, because this might provide us the convergence analysis as well. --- This appraoch is for example elaborated in
On the Finite Element Solution of the Pure Neumann Problem
Pavel Bochev and R. B. Lehoucq
SIAM Review
Vol. 47, No. 1 (Mar., 2005), pp. 50-66
Published by: Society for Industrial and Applied Mathematics
Article Stable URL: http://www.jstor.org/stable/20453601
Apart from this canon standard example for the Laplacian, system matrices with larger kernel appear in numerical methods for the de-Rham-complex, in particular if the domain is topologically non-trivial (Finite Element Exterior Calculus, Discrete Exterior Calculus). Singular system solves are still no standard material for education in computational science. As far as I may dare to give an estimate, there is still much room for a better theory building.
If the matrix is singular then then the images of some vectors may be zero and that could mean some of the steps may end up with division by zero.
The conjugate gradient method becomes unstable when the matrix A is singular, i.e., if you compute the output of the k-th iteration for a small change in the initial data, large deviations will occur. A similar large deviation occurs due the computer roundoff error. Even for non singular matrices A that are close to to be singular, for example when one of the eigenvalues of A is smaller by orders of magnitudes than the rest, instability problems can occur: The convergence rate becomes slower, and even convergence can be lost due to roundoff errors.
