Don’t invert that matrix

http://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/

by KonaB • 17 years ago

90 38 17 years ago

jacobolus• 17 years ago

Does anyone trying to solve problems involving matrices with millions of entries not have some idea of the computational complexity of various operations they can do on them, given the nature of their matrices, etc.? I’d expect that either (a) someone has a small enough problem that the computational time is irrelevant and they should just arrive at the solution in whatever manner is conceptually clearest, given their current mathematical knowledge, or (b) the problem is going to take a lot of computation time, in which case it’s worth exploring to learn exactly what the fastest method will be.

lliiffee• 17 years ago

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.

jedbrown6• 17 years ago

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.

liuliu• 17 years ago

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/

Retric• 17 years ago

Speed is not the only measurement, accuracy can also be a problem.

roundsquare2• 17 years ago

Not necessarily. The specific mathematical/programming techniques can be very complex. It is probably not reasonable to learn them just because you are dealing with large matrices. Of course, some knowledge of the computational complexity is definitely helpful.

memetichazard• 17 years ago

I'm a bit rusty on my linear algebra. Without inverting the matrix, what's the algorithm for solving Ax=b? Factorization is mentioned but I'm not sure what it means in this case.

ajross3• 17 years ago

LU decomposition is the most common choice, which is a similar operation to matrix inversion in that you're iterating row (or column) operations on both sides of an equation. Numerical Recipes has a great section on this, IIRC.

But that kind of obscures the fact that, while algorithms like LU decomposition are faster than a pure matrix inversion, they're faster by a constant factor of about 2 or so. It's not a fundamentally different algorithm, just a different set of tricks to solve a bunch of simultaneous equations. If you aren't hugely performance constrained and already have a matrix inversion routine, you'd never be bothered to re-write your equation solver just for this speedup.

eru2• 17 years ago

As far as I know LU decomposition is more than a constant factor faster than inversion. And it is numerically more stable.

JabavuAdams• 17 years ago

Gaussian elimination. Remember where you subtract multiples of one row from another, in order to get zeros below the diagonal? That way you can back-solve a system of equations by reading off solutions starting with the bottom-right corner.

The result of carrying out the elimination steps is a triangular matrix with all zeros below the main diagonal. Call it U, for "upper".

The elimination steps themselves can be encoded as a matrix, rather than describing them in words. For instance, "subtract 2 times row 1 from row 2" is the matrix:

[ 1 0 0]

[-2 1 0]

[ 0 0 1]

So, to solve the system, you can perform elimination on A x = b, which gives you two matrices (U, and the elimination steps L). In some sense, you've factored A in to L and U.

Example:

Using Strang's terminology, let's call the matrix that "eliminates (i.e. sets to 0) the (i,j) entry" Eij

So, to solve a 3x3 system A x = b, you want to do elimination as: (E32 E31 E21) A = U

I.e.

1) Take A and subtract some multiple of row 2 from row 1, to put 0 in (2,1)

2) Take the result and subtract some multiple of row 3 from row 1 to put 0 in (3,1)

3) Take the result and subtract some multiple of row 3 from row 2 to put 0 in (3,2)

Now the goal of elimination is to get an uppper-triangular matrix U. I.e. if you were solving some system of equations in x, y, z, then because U is triangular, you could just read off the answer to z. Then plug that into the equation for row 2 to solve for y, then plug that into the equation for row 1 to solve for x. That's back-substitution.

BUT, let's return to our operations. We have:

E32 E31 E21 A = U

(E32 E31 E21) A = U

inverse(E32 E31 E21) (E32 E31 E21) A = inverse(E32 E31 E21) U

I A = inverse(E32 E31 E21) U

Now, it just so happens that the inverse of (E32 E31 E21) is lower-triangular. I.e. all entries above the diagonal are zero, so let's call it L.

I A = L U

A = L U

LU-decomposition. Or "factorization", since we've factored A into two matrices: L, and U.

whyenot2• 17 years ago

Yes, but it's a little more complicated than that. In most cases you cannot trust the solution of straight Gaussian elimination. You need to take additional measures to avoid excessive error growth, such as using partial pivoting.

see http://en.wikipedia.org/wiki/Pivoting and http://en.wikipedia.org/wiki/Gaussian_elimination

SandB0x• 17 years ago

As others have said, Gaussian elimination, or if the matrix is not invertible (square and singular or not square) then the Moore-Penrose pseudoinverse can be used to find the least squares solution.

jedbrown6• 17 years ago

Except that the normal equations are rarely a good way to compute the action of the pseudoinverse (though explicitly computing it would be even worse). QR is the workhorse for under- and over-determined dense systems, with SVD as a fallback for especially nasty ones.

imurray2• 17 years ago

A Matlab/Octave example of inv being slower and less accurate:

  >> A = randn(100); x = randn(100,1); b = A*x;
  >> tic; x1 = inv(A)*b; toc
  Elapsed time is 0.142381 seconds.
  >> tic; x2 = A\b; toc     
  Elapsed time is 0.000736 seconds.
  >> max(abs(x-x1))
  ans =
     7.1054e-14
  >> max(abs(x-x2))
  ans =
     4.8406e-14

Edit: Actually this is a warning about timing stuff. The inv function hadn't "warmed up". On a second run both times drop, and the difference is much smaller:

  >> tic; x1 = inv(A)*b; toc
  Elapsed time is 0.000737 seconds.
  >> tic; x2 = A\b; toc
  Elapsed time is 0.000471 seconds.

ibsulon• 17 years ago

What happens when you run it 1000 times, skipping the first?

imurray2• 17 years ago

If I use this:

http://www.mathworks.com/matlabcentral/fileexchange/18798-ti...

which warms up the function, averages over repetitions and tries to account for various overheads, the ratio states about the same (1:1.5 ish). I would normally use timeit(), but for a quick post, I was being sloppy.

ewjordan• 17 years ago

Here's another piece of advice that covers a lot more than matrix inversion: don't solve math problems yourself in code if you can in any way avoid it.

Use a library instead, and use the most specific methods that are available: if there are specialized routines to solve Ax=b, use them, don't use the matrix inversion methods (which are probably also there), or even the LU decomposition ones, even if Real Math tells you that the result is equivalent - in numerical computation, mathematical equivalence is an entirely different entity than what "real" mathematicians concern themselves with.

Seriously. I don't care if you're a math PhD, I don't even care if you know a decent amount about numerical methods. Whatever code you're going to put together to do the job will be inferior to what a mature library can do, simply because the library has been put through many iterations and has had a lot more time to get patched up to handle all those special cases (usually related to finite precision or discretization) that Real Math doesn't have to worry about. And it will probably be vastly more optimized than anything you'll put together, because it's spent years of heavy use in some critical research environments.

The only time you should be writing your own math code (beyond simple arithmetic) is if no mature libraries are available to you. In which case you should allocate at the very least twice the amount of time it will take to research and implement the methods, and probably a lot more, because you're going to have a lot of fiddling to do.

DrJokepu• 17 years ago

Another reason would be (self-)education. There's no better way to gain a deeper understanding of a subject or learn the ins-and-outs of optimization than writing some code and then possibly trying to improve it. If you always blindly depend on third party libraries, you will remain ignorant and dumb.

Then, you can go on and use third party libraries in production while feeling satisfied that you know what's going on inside the library.

gaius• 17 years ago

If you always blindly depend on third party libraries, you will remain ignorant and dumb.

Depends if your goal is to become an aeronautical engineer or an assembly language hacker, really. Neither one's "better" than the other, but the the former won't benefit much from understanding exactly how to optimize on SPARC vs x64... Nor will the latter need to worry about his compiler's stalling speed :-)

nitrogen• 17 years ago

What sort of keywords would one enter into a search engine to find such a library?

gjm11• 17 years ago

<<<linear algebra library>>> finds several good ones with Google.

[EDITED to add: That's if what you want is specifically linear algebra. <<<numerical library>>> does pretty well if your needs are more general.]

metaguri• 17 years ago

The other important reason to use a library is that if someone improves the algorithm someday, you can benefit for nothing more than installing and testing the new version of the library.

I remember this also being the case for Constraint Satisfaction Problem (CSP) algorithms in AI. If you can formulate your problem as a CSP and feed it to a library, you automatically benefit from future improvements in the field of CSP solvers.

brettnak• 17 years ago

Are people interested enough in this for a blog post? I studied numerical analysis and would be happy to post some of the basics if people are interested.

nimrody• 17 years ago

What is usually missing from these articles is even a simplified error analysis.

We can all count the floating-point operations, and today caches sometimes influence matrix algorithms more than raw flop counts. The problem gets more interesting when the matrix becomes ill conditioned and different algorithms will produce very different results.

indigoviolet• 17 years ago

Yes! Please.

profquail• 17 years ago

I did numerical analysis too...I don't think most CS departments cover it though, so chances are a lot of people would be interested in it.

sown• 17 years ago

In the same vein as "Statistics Hacks" I'm sure those same fans will be fond of "Numerical Hacks". :)

roundsquare2• 17 years ago

I'd be interested, but probably if you went into less common subjects. "Don't invert that matrix" is probably the first thing I, or anyone else, learned in numerical analysis.

Side note: My professor also mentioned it at least 3 times in each lecture... so it is important.

Don’t invert that matrix

Related Stories