0
$\begingroup$

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.

| cite | improve this question | |
$\endgroup$
  • 2
    $\begingroup$ Did you compute the residual for your solution? $\endgroup$ – nicoguaro♦ Dec 3 '17 at 23:44
  • $\begingroup$ I will do that. But I read that even if the b-A*x is small, it doesn't necessarily mean x is accurate? $\endgroup$ – David Dec 4 '17 at 0:23
  • 2
    $\begingroup$ 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. $\endgroup$ – Brian Borchers Dec 4 '17 at 4:03
  • $\begingroup$ @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?" $\endgroup$ – David Dec 4 '17 at 5:02
  • 1
    $\begingroup$ 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). $\endgroup$ – origimbo Dec 4 '17 at 18:29
3
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – David Dec 5 '17 at 5:45
  • $\begingroup$ 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. $\endgroup$ – Brian Borchers Dec 5 '17 at 6:18
  • $\begingroup$ thanks. I will see if I can figiure out how to implement Tikhonov regularization for my problem. $\endgroup$ – David Dec 6 '17 at 0:52
  • $\begingroup$ 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. $\endgroup$ – Brian Borchers Dec 6 '17 at 1:07
  • $\begingroup$ 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? $\endgroup$ – David Dec 6 '17 at 18:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.