# How to verify solution to pre-conditioned linear systems solver?

I am solving Ax=b. A has a very large condition number (> O(10^10))

I am using the conjugate gradients method with point jacobi pre-conditioning. I obtained a solution 'x' that "looks" reasonable. What is the best way to verify that my solution is correct using this method with preconditioning?

Note that I cannot obtain an exact solution to the original equation.

My matrix is symmetric, and I believe it is negative definite/semidefinite as at least one of the its eigenvalues is 0.

• Did you compute the residual for your solution? Dec 3 '17 at 23:44
• I will do that. But I read that even if the b-A*x is small, it doesn't necessarily mean x is accurate? Dec 4 '17 at 0:23
• Your system is very ill-conditioned, so there are lots of solutions that are basically equivalent in terms of having a small residual. Unless you regularize the problem, the answer is that there isn't just one "correct" solution- this is a fundamental property of the ill-conditioned system and not of the algorithm you've used to solve the system. Dec 4 '17 at 4:03
• @BrianBorchers Can "regularization" be done by pre-conditioning alone? If pre-conditioned, are there still "lots of solutions that are basically equivalent in terms of having a small residual?" Dec 4 '17 at 5:02
• Regularisation isn't frequently directly implemented in numerical toolkits, since the choice of regularisation term will depend on the underlying physics the equation is meant to capture. For example, you could be looking for "smooth" solutions (in which case it might look like a derivative operator) or ones with small norm (in which case you might use a scalar multiple of the identity matrix). Dec 4 '17 at 18:29

You've started with a singular linear system of equations $Ax=b$. As a practical matter, it's unlikely that $b$ lies exactly in the range of $A$, so at best you can find a least squares solution that minimizes

$\min \| Ax - b \|_{2}$

Because the system is singular, the null space of $A$ is non-empty, and there will be an infinite number of solutions to this least squares problem. In particular, if $x_{LS}$ is some particular least squares solution and $v$ is any vector in the null space of $A$, then any $x$ of the form $x=x_{LS}+\alpha v$ will be a least squares solution.

Which of these infinitely many least squares solutions is the "correct" one?

Perhaps you don't care and will accept any least squares solution. In that case, any solution that satisfies the normal equations $A^{T}Ax=A^{T}b$ will be correct.

Perhaps you want to simultaneously minimize the norm of $x$ to get a minimum norm least squares solution? You can regularize the least squares problem by minimizing

$\min \| Ax-b \|_{2}^{2} + \lambda^{2} \| x \|_{2}^{2}$

where $\lambda$ is a small parameter. This least squares problem is strictly convex, so it has a unique optimal solution.

What happens in floating point arithmetic with limited precision? In that case, it's likely that $A$ will have a very large but still finite condition number.

You can throw a preconditioner $M$ at $A$, in hopes that $M^{-1}A$ will be better conditioned than $A$. This might work in the sense at $M^{-1}A$ is numerically better conditioned than $A$, but since we know that $A$ is theoretically singular, this is a mirage! $M^{-1}A$ should still be singular.

So what is the solution to the numerically well-conditioned system $M^{-1}Ax=M^{-1}b$ doing? This is similar to regularization in the sense that $M^{-1}A$ is better conditioned than $A$, but you don't know much about the direction in which the preconditioner has pushed the solution. If $M^{-1}Ax-M^{-1}b$ is small, and $M^{-1}$ is well behaved, then perhaps $Ax-b$ is small and you've got a least squares solution. On the other hand, if $M^{-1}$ is very badly scaled, it may be that the residual in the original problem $Ax-b$ is quite large.

There has been some research on preconditioners that are designed to achieve a particular regularization. See for example:

Calvetti, Daniela. 2007. “Preconditioned Iterative Methods for Linear Discrete Ill-Posed Problems from a Bayesian Inversion Perspective.” Journal of Computational and Applied Mathematics, Special Issue: Applied Computational Inverse Problems, 198 (2):378–95. https://doi.org/10.1016/j.cam.2005.10.038.

• I just want to verify my understanding if regularization. We can use regularization to obtained a "preferred" solution, but the "preferred" solution is still not unique? Dec 5 '17 at 5:45
• The regularized solution will be unique and furthermore, the regularized problem will be numerically well conditioned so that it can reasonably be solved in floating point. Dec 5 '17 at 6:18
• thanks. I will see if I can figiure out how to implement Tikhonov regularization for my problem. Dec 6 '17 at 0:52
• In terms of PETSC, you can use the LSQR implementation that's already in PETSC. The to combine the $\| Ax -b \|_{2}^{2}$ and the $\lambda \|x \|_{2}^{2}$ into one least squares problem involving the matrix $[A; \sqrt{\lambda} I]$. Don't explicitly form this matrix but rather provide a routine that does the matrix vector multiplication. Dec 6 '17 at 1:07
• You mentioned that the Regularization parameter is small. If I choose between a “range” of small values, should I expect to obtain the same regularized LS solution? Dec 6 '17 at 18:14

If you don't care too much about which $$b$$ you're working with (say just for linear algebra), then you can use the Method of Manufactured Solutions (MMS) in Linear Algebra much like we differential equations geeks would for our problems to start with a known exact solution, $$x^*$$, derive a $$b$$, and then apply your solver to see how decent an $$x$$ you can solve for. By picking $$b=Ax^*$$ for any $$x^*$$ you make up, you can see how well your solver performs on that $$b$$ and ask yourself as many questions as you like about how your proposed solution $$x$$ and $$x^*$$ compare. If you require that $$b$$ or $$x^*$$ have certain properties, then you may have more trouble. This way you can compute $$\|{x-x^*}\|$$ in almost any norm you like. If you're committed to your $$A$$ and $$b$$, then the problem is somewhat harder, and I'd suggest using the MMS in the original problem that led to $$A$$ and $$b$$ if you can control them (which is what we PDE geeks do when we can). I can't tell for sure if you're trying to pick the best solver for your underlying problem or if you're trying to design a new linear algebra solver for a class of poorly conditioned problems. Either way, you can use MMS to give you exact solutions to work with. For your PDE, the MMS $$x^*$$s might not be highly physical, but the $$b$$ that is generated for you will be as representative as possible for what your system, $$A$$, can create for $$x^*$$.

I recommend starting with simple $$x^*$$s and working up to complex ones. In the end, for error analysis, people tend to recommend using functions that excite as many of the possible modes of the physical solution as possible and have no singularities (unless you expect them), like transcendentals and other high-order functions, but starting with polynomials on simple domains is probably best when it comes to picking functions for arbitrary, smooth PDEs.

The original paper on MMS can be found at the original MMS paper link (PDF warning, if you care).