Skip to content

Comment on Don’t invert that matrixparent

Comments

You can solve amazingly large problems with a sparse matrix with out understanding what you are doing, as long as you use direct (factorization) methods. If you type A\b with some gigantic sparse A, matlab will apply reordering heuristics that often just work, and are actually incredibly hard to improve upon by manual fiddling. When you get to really, really big problems, you need to use iterative methods (like conjugate gradient), and for those you have to have a very firm understanding of everything. The gap between direct and indirect methods is huge.

Matlab's backslash is the GPL'd Umfpack (relicensed from the author). The most important thing to know about iterative solvers is that the Krylov method is a minor component, preconditioning is the hard part.

You need to have a basic understanding of the different Krylov subspace methods to be able to use them. If you are using CG on a matrix that isn't SPD, then you are doing something very wrong. Same goes for the other methods -- they all have their strengths, weaknesses and gotchas, and it's pretty important to know what they are if you are going to use them! Preconditioning can greatly reduce the number of iterations you have to do once you have chosen a method, but it's not particularly helpful if you are using the wrong method.

With CG, the preconditioner also needs to be SPD (or have suitable factors). You can get by pretty well using only GMRES and spend all your time thinking about preconditioners. Once you have some good candidate preconditioners, you can try other Krylov methods, but the gains are rarely very large. I have frequently found SPD systems that are more efficiently solved with GMRES and a non-symmetric preconditioner, and even more so for symmetric indefinite systems where popular wisdom says to use MINRES, CR, or SQMR.

For amusement, it's really worth knowing about

http://www.comlab.ox.ac.uk/people/nick.trefethen/publication...

which shows that each of the popular nonsymmetric iterations can beat all others by a huge factor.

But this is all justification that you should use a library like PETSc for your linear solvers. Then all aspects of the solve become run-time (typically command-line) options including an algebra for composing preconditioners and other components, with a plugin architecture for every component.

If you are using CG on a matrix that isn't SPD, then you are doing something very wrong.

Actually, this isn't so weird. There are a number of domain decomposition schemes that use CG, but keep all iterates in a benign space on which the operator is SPD, even though the real system is indefinite (porous media and Stokes/mixed elasticity). See work by Axel Klawonn, and recent stuff from Xuemin Tu.

When I working on l1-minimize problem, I wrote a very simple cg solver here and it can solve for a sparse matrix of A (Ax=B where x and B are dense) (c with cxcore from opencv):

http://github.com/liuliu/l1qc/tree/master/src/

AboutSource Built by g1lg1l

Hackerly is an independent reader for Hacker News, built on the public HN API. Not affiliated with Y Combinator.