1

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.

David
  • 251
  • 2
  • 8
  • 2
    Did you compute the residual for your solution? – nicoguaro Dec 03 '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? – David Dec 04 '17 at 00:23
  • 2
    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. – Brian Borchers Dec 04 '17 at 04: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?" – David Dec 04 '17 at 05:02
  • Preconditioning doesn't change the condition number of the underlying problem. It will push the solution that you get in some not particularly well-controlled direction. Regularization strategies (such as Tikhonov regularization) will push the solution towards some understandably preferred solution among the many equally good ones (e.g. by finding a minimum norm approximate least squares solution.) – Brian Borchers Dec 04 '17 at 05:25
  • I'll add that there has been research done on using specially designed preconditioners to achieve particular regularization effects. See for example this paper by Daniela Calvetti: https://www.sciencedirect.com/science/article/pii/S0377042705007211 – Brian Borchers Dec 04 '17 at 05:50
  • @BrianBorchers Thank you for the paper. Read it all. I have honestly never heard of regularization. Based on your comments, it seems regularization appear to be the superior technique to obtain an accurate solution for an originally ill-conditioned problem. But pre-conditioning seems to be more commonly used. Why is that? – David Dec 04 '17 at 06:17
  • I am using the PETSc library to solve my Ax=b system. It does not seem the Tikhonov Regularization is a capability within this library. – David Dec 04 '17 at 06:18
  • 1
    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). – origimbo Dec 04 '17 at 18:29
  • @origimbo I am solving the linear elasticity steady state equations. https://scicomp.stackexchange.com/questions/21476/which-preconditioning-for-large-linear-elasticity-problem – David Dec 04 '17 at 18:33
  • My only experience of that equation set is solving it in a "nice" parameter/boundary condition space, but if you have Neumann boundary conditions, it's probable you want to look at a term penalising rigid body motions. – origimbo Dec 04 '17 at 18:55
  • I am using finite volumes to solve this equation. I also have a thermal stress term implemented into the original equation. My boundary conditions are currently only Dirichlet, and I will be working on Neumann BCs later on. Have you had experience with adding this thermal stress term? My thermal stress term is factored into the forcing vector 'b' and does not impact the Jacobia matrix A. – David Dec 04 '17 at 18:58
  • My current dirichlet BC is very sample. It uses symmetry, i.e., all boundary cells' values equal the adjacent interior cell value. What drives the change is the thermal stress (temperature is changing). – David Dec 04 '17 at 19:00

2 Answers2

5

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.

Brian Borchers
  • 18,719
  • 1
  • 39
  • 70
  • 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? – David Dec 05 '17 at 05: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. – Brian Borchers Dec 05 '17 at 06:18
  • thanks. I will see if I can figiure out how to implement Tikhonov regularization for my problem. – David Dec 06 '17 at 00: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. – Brian Borchers Dec 06 '17 at 01: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? – David Dec 06 '17 at 18:14
  • No- the choice of the regularization parameter will effect the solution that you get. You have to manage the trade-off between the least squares objective and the regularization penalty. There is a huge literature on ways to do this. – Brian Borchers Dec 06 '17 at 19:32
  • I'll review the literature. I was speaking with another research group who is solving the same equations as me using the finite volume method, and they claim that their A matrix is well-conditioned. It puzzles me as to how his system is well-conditioned.I am now suspecting that it could possibly be with how I discretized my gradients (I used Discrete Green Gauss theorem). in the steady state linear elasticity equations. Perhaps this discretization technique led to the ill-conditioned A matrix. – David Dec 06 '17 at 22:10
2

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).

Bill Barth
  • 10,905
  • 1
  • 21
  • 39