conjugate gradient iteration - MathOverflow most recent 30 from http://mathoverflow.net 2013-05-26T04:06:21Z http://mathoverflow.net/feeds/question/7903 http://www.creativecommons.org/licenses/by-nc/2.5/rdf http://mathoverflow.net/questions/7903/conjugate-gradient-iteration conjugate gradient iteration john 2009-12-05T21:51:52Z 2011-09-17T22:17:33Z <p>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 <a href="http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf" rel="nofollow">CG method</a> 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?</p> <p>Thanks in advance.</p> http://mathoverflow.net/questions/7903/conjugate-gradient-iteration/7940#7940 Answer by Kristal Cantwell for conjugate gradient iteration Kristal Cantwell 2009-12-06T01:40:19Z 2009-12-06T01:40:19Z <p>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. </p> http://mathoverflow.net/questions/7903/conjugate-gradient-iteration/7943#7943 Answer by Darsh Ranjan for conjugate gradient iteration Darsh Ranjan 2009-12-06T02:20:12Z 2009-12-06T02:20:12Z <p>Let my reply here to your comment:</p> <blockquote> <p>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.</p> </blockquote> <p>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. </p> <p>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$. </p> http://mathoverflow.net/questions/7903/conjugate-gradient-iteration/7990#7990 Answer by David Bar Moshe for conjugate gradient iteration David Bar Moshe 2009-12-06T11:37:06Z 2009-12-06T11:37:06Z <p>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. </p> http://mathoverflow.net/questions/7903/conjugate-gradient-iteration/75704#75704 Answer by Martin for conjugate gradient iteration Martin 2011-09-17T22:17:33Z 2011-09-17T22:17:33Z <p>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?</p> <p>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.</p> <p>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</p> <pre><code>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 </code></pre> <p>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.</p>