3
$\begingroup$

I'm trying to solve

$Ax = b$

given a vector $b$ and a large, symmetric positive definite, sparse, banded matrix $A$ that has a very poor condition number.

I know several iterative methods that could be used for this task (but haven't implemented them all):
a) (Preconditioned) Conjugate Gradient.
b) Simple Gradient Descent.
c) SSOR / Gauss-Seidel.
d) Multigrid flavors of the above.

Also, I know that there are sparse direct solvers but I have been reluctant to try because they seem very complex and hard to implement.

So what's the method of choice here?

|cite|improve this question
$\endgroup$
  • $\begingroup$ Yes. May I just repost there? $\endgroup$ – Rafael Jan 19 '12 at 19:14
  • $\begingroup$ Yes, but provide links between the questions. $\endgroup$ – Thomas Klimpel Jan 19 '12 at 20:08
  • $\begingroup$ @Rafael: if you repost it there, the moderators can do some voodoo and link the two questions, as long as you let us know the url of the version over at scicomp. $\endgroup$ – Willie Wong Jan 20 '12 at 13:12
  • $\begingroup$ Thanks. I did some research on scicomp and noticed that there are very similar questions already. Therefore didn't want to post. $\endgroup$ – Rafael Jan 20 '12 at 14:14
  • $\begingroup$ @Rafael What sort of physical problem does this come from? What causes the ill-conditioning? How would you characterize the near-null space? How large are you target problems? $\endgroup$ – Jed Brown May 20 '12 at 21:04
2
$\begingroup$

If your bandwidth is small, a direct method will be by far the fastest, and if implemented with iterative refinement can provide the solution to almost working accuracy unless the condition is so poor that no method provides any accuracy.

As your matrix is positive definite, it is easy to implement a direct banded method - because of poor condition, increase all pivots smaller than 13-8 times the largest diagonal element to this value before continuing the factorization to increase numerical stability. I would use the resulting factorization as the preconditioner of a conjugate gradient method, then you get the best of both worlds.

|cite|improve this answer
$\endgroup$
  • 2
    $\begingroup$ Note that bandwidth is not the critical factor for sparse direct solver performance, it is the size of vertex separators, and a low-bandwidth ordering will produce much more fill than necessary. $\endgroup$ – Jed Brown May 20 '12 at 18:08
  • $\begingroup$ @JedBrown: If the matrix is banded (as the question implies) and full between the bands (as I assumed in LAPACK manner) the optimal ordering is the original band, and there is no fill in. $\endgroup$ – Arnold Neumaier May 21 '12 at 11:34
  • $\begingroup$ This depends how close the band is to the diagonal. A query reveals that essentially all use of the term "banded" on this site involve PDEs on structured grids which are decided not banded in the LAPACK definition. Due to this, and years of the same on the PETSc mailing lists, I disregard claims of "banded" structure until the questioner is precise about what they mean. $\endgroup$ – Jed Brown May 21 '12 at 11:42
1
$\begingroup$

There are several other better and more complicated iterative methods (gmres, bicgstab, etc.), for example here is a list of methods that Matlab supports. If the matrix is very large, direct methods could do very poorly, depending on the sparsity structure.

But why do you want to implement the method yourself in the first place? If the system is large and poorly conditioned it is unlikely that you will be able to do better yourself than to use a highly-optimized method already included in a mature library.

|cite|improve this answer
$\endgroup$
  • 2
    $\begingroup$ The original poster said $A$ is spd (symmetric positive definite), so gmres and bicgstab will perform worse than (Preconditioned) Conjugate Gradient. $\endgroup$ – Thomas Klimpel Jan 18 '12 at 23:53
  • $\begingroup$ You're right, I agree. $\endgroup$ – Kirill Jan 19 '12 at 0:07
  • $\begingroup$ I don't necessarily want to implement solvers myself. However, many solver libraries I saw also seemed very complex to integrate so I wanted to get a better understanding of the problem before I give a shot at a particular solver/lib, taking lots of work to integrate. $\endgroup$ – Rafael Jan 19 '12 at 19:12
  • $\begingroup$ @Rafael What environment are you working in? Libraries are a good place to learn about methods and make it trivial to experiment with many classes of methods. $\endgroup$ – Jed Brown May 20 '12 at 19:24
0
$\begingroup$

Best method would be MG preconditioned CG. I believe that is the ultimate solver for this case. If you know how the system is generated, i.e the geometry is known, you can apply geometric multigrid. If not, algebraic multigrid should be used. Surely, for small problems - direct methods could be favored yet even in that case MG+PCG is really fast -if properly implemented, of course.

|cite|improve this answer
$\endgroup$

Your Answer

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