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.