How can a CG solver solve a non positive definite sparse matrix

Question

I am using the CUSP CG solver and I ran it on a couple of sparse matrices from the University of Florida sparse matrix collection. The solver was able to solve non positive definite sparse matrices. My understanding is that CG solvers can't solve non positive definite sparse matrices? So is how was CUSP's CG solver able to do it?

Answer 1

I highly recommend the following read:

• J.R. Shewchuk, "An Introduction to the Conjugate Gradient Method Without the Agonizing Pain"

In short, if the matrix is non-positive definite, there is no guarantee that CG will fail. It might be able to solve it (for some RHSs and certain tolerances), but it is just not supposed to and, probably, converges slower than the appropriate methods.

You definitely should look into the other iterative solvers (say, BiCGStab or something else, depending on your problem properties).

NB: I am not familiar with CUPS, so I would also check if there is a function inside that checks the property of the matrix and switches to some other algorithm (unlikely, but possible) – especially if the solver log is ambiguous.

Answer 2

In addition to Anton's answer I would like to provide a special case where CG is able to solve a system even if it is non-positive definite.

When simulating incompressible fluid, and if the fluid domain is surrounded with pure Dirichlet BCs (i.e. only boundary velocities are specified, and we don't have any free surface or boundary pressure), then we have a nullspace in the system. Thus the matrix is non-positive definite.

That nullspace is all ones, $\mathbf{e}=(1, 1, \cdots, 1)^T$. Which is not surprising because we only have pressure gradient $\nabla p$ in the equation. If $p$ is a solution to $\nabla p = 0$, then so is $p + c*\mathbf{e}$, where $c$ is any constant.

Since the nullspace is known, we can project out that nullspace whenever we do the matrix multiplication in CG: $$ \begin{align} \mathbf{k} &= \mathbf{A} \mathbf{x} \\ \mathbf{k}_{\perp} &= \mathbf{k}-(\mathbf{k}\cdot \mathbf{e}) \mathbf{e}/N \end{align} $$ ($N$ is the width or height of the square matrix $\mathbf{A}$)

Then with just this additional projection, the non-positive definite system can still be solved by CG.