Hacker News Board Hacker News, browsed like a classic message board
Live data from Hacker News

Don’t invert that matrix (2010)

johndcook.com

Post reply on HN
Print view View original

11–20 of 48 posts

Re: Don’t invert that matrix (2010)

#11

Post by tclmeelmo »

jacobolus wrote: post #8
mwfunk wrote: post #4

I get that inverting a matrix can be computationally intensive, but what's the alternative? "What if you have to solve Ax = b for a lot of different b‘s? Surely then it’s worthwhile to find A-1. No. The first time you solve Ax = b, you factor A and save that factorization. Then when you solve for the next b, the answer comes much faster." I don't get it, although I'm not a numerics expert so maybe I'm missing somethi…

For example, you factorize A into LU, so that you have LUx = b. Now you can let Ux = y, solve the system Ly = b for y, then solve Ux = y for x. Solving each triangular system is fast: O(n^2). If you need to do this with a different b, you can do it again, and you already have A factorized, which was the expensive part, i.e. O(n^3). (Depending on the type of problem, you might want to factorize the original matrix in…

IIRC, the HP-15C used LU decomposition in its matrix math routines.

Re: Don’t invert that matrix (2010)

#12

Post by [deleted] »

jacobolus wrote: post #8
mwfunk wrote: post #4

I get that inverting a matrix can be computationally intensive, but what's the alternative? "What if you have to solve Ax = b for a lot of different b‘s? Surely then it’s worthwhile to find A-1. No. The first time you solve Ax = b, you factor A and save that factorization. Then when you solve for the next b, the answer comes much faster." I don't get it, although I'm not a numerics expert so maybe I'm missing somethi…

For example, you factorize A into LU, so that you have LUx = b. Now you can let Ux = y, solve the system Ly = b for y, then solve Ux = y for x. Solving each triangular system is fast: O(n^2). If you need to do this with a different b, you can do it again, and you already have A factorized, which was the expensive part, i.e. O(n^3). (Depending on the type of problem, you might want to factorize the original matrix in…

[deleted]

Re: Don’t invert that matrix (2010)

#13

Post by GFK_of_xmaspast »

vvanders wrote: post #9
rcxdude wrote: post #3

Earlier quoted context omitted.

More talking about applications where the matrices are large, like statistics, machine learning (but I repeat myself), signal processing, etc. In games pretty much all matrices are around 4x4, instead of millions by millions.

Ah, that's the context I was missing, thanks! In games we can usually just transpose the inner 3x3 R+S, invert scale and handle translation similarly to get a semi-quick inverse.

This is the difference between matrices that are SO(3) (which have trivial inverses) and matrices that are just GL(n) (which: don't).

Re: Don’t invert that matrix (2010)

#14

Post by [deleted] »

vvanders wrote: post #9
rcxdude wrote: post #3

Earlier quoted context omitted.

More talking about applications where the matrices are large, like statistics, machine learning (but I repeat myself), signal processing, etc. In games pretty much all matrices are around 4x4, instead of millions by millions.

Ah, that's the context I was missing, thanks! In games we can usually just transpose the inner 3x3 R+S, invert scale and handle translation similarly to get a semi-quick inverse.

[deleted]

Re: Don’t invert that matrix (2010)

#15

Post by Someone »

pmarreck wrote: post #7

I am not a mathematician, but simply based on the fact that beautiful/simpler answers to questions seem to get closer to the ultimate truth, the current state of large-dimension matrix inversion is far from the ideal possible solution, because it is ugly as hell.

It's not the current state; it is inherent in the problem, so the situation will never get better.

Applied math is like a writing a program for customers; the result is more important than the beauty of the road towards it.

Re: Don’t invert that matrix (2010)

#16

Post by ThatGeoGuy »

pmarreck wrote: post #7

I am not a mathematician, but simply based on the fact that beautiful/simpler answers to questions seem to get closer to the ultimate truth, the current state of large-dimension matrix inversion is far from the ideal possible solution, because it is ugly as hell.

Arguably LU / LDLT factorizations before solving the direct system are the more elegant / beautiful / simpler answers if you take the algorithm at face value. However, large-dimension matrix inversion has problems outside of just beauty. Part of the problem is that we don't have fast, perfect representations of real or complex numbers. At best we have double precision floating points, which means that when you take the adjoint over the determinant, you're looking at a lot of error.

That's not to say that the implementations of these have to be ugly or that our current BLAS / LAPACK / whatever can't be improved. There's a lot of room for improvement in this space, but the fact that they're 1. fast and 2. well tested and robust means that improving these sorts of interfaces come with baggage. People complain about the wheel being reinvented every time somebody invents a new Javascript framework, but could you imagine if this same problem applied to BLAS and other high performance math libraries? (Arguably, at least then we'd have more choices in math space)

Re: Don’t invert that matrix (2010)

#17

Post by andy_wrote »

In R there is a little syntactical nudge in this direction in that the matrix inversion function is actually called `solve`. `solve(A, B)` gives you the solution to AX = B, and B's default argument is the identity, so `solve(A)` just gives you the inverse, and any instances of `solve(A) %*% B` in your code should stand out as red flags. As rcxdude mentions this is something that you tend to care more about in statistical applications, and that's what R is for.

Re: Don’t invert that matrix (2010)

#18

Post by theoh »

rcxdude wrote: post #3
vvanders wrote: post #2

Is there more context around what exactly he's getting at? In games we use matrix inverses all the time(and since they're usually a specialized affine transform there's cheap shortcuts to calculate them) for stiff like hit detection or animation.

More talking about applications where the matrices are large, like statistics, machine learning (but I repeat myself), signal processing, etc. In games pretty much all matrices are around 4x4, instead of millions by millions.

There's an neat technique for fitting RBF representations to dense point clouds. Instead of the impractical matrix inversion, the trick is to solve the linear system with a fast multipole method, similar to the Barnes-Hut algorithm for the n-body problem.

Carr, J.C. et al, "Reconstruction and Representation of 3D Objects with Radial Basis Functions", Siggraph 2001

http://dl.acm.org/citation.cfm?id=383266

Re: Don’t invert that matrix (2010)

#19

Post by mturmon »

askee wrote: post #5
vvanders wrote: post #2

Is there more context around what exactly he's getting at? In games we use matrix inverses all the time(and since they're usually a specialized affine transform there's cheap shortcuts to calculate them) for stiff like hit detection or animation.

In general there are no shortcuts, and when the matrix A is not small (e.g. m x n with m,n << 100), you don't want to invert A, especially if you are only interested in x (from A x = b). Instead, use numerical minimization schemes like conjugated gradients and variants thereof or Quasi-Newton methods (BFGS). Combined with preconditioning as well as regularization or denoising this usually yields good results. Compres…

This comment seems to be a diversion from the question asked. The main point of OP is that solving A x = b is commonly done with factorizations of A, not by explicit computation of inv(A).

Resorting to a numerical minimization (like MINRES) would be unusual unless the dimension is much bigger than hundreds. Certainly it's not relevant in a graphics computation with dimension of 3, 4, or 6. My laptop inverts a 4Kx4K general matrix in 1.5s. It solves a 4K linear system in 0.5s.

Re: Don’t invert that matrix (2010)

#20

Post by IngoBlechschmid »

mwfunk wrote: post #4

I get that inverting a matrix can be computationally intensive, but what's the alternative? "What if you have to solve Ax = b for a lot of different b‘s? Surely then it’s worthwhile to find A-1. No. The first time you solve Ax = b, you factor A and save that factorization. Then when you solve for the next b, the answer comes much faster." I don't get it, although I'm not a numerics expert so maybe I'm missing somethi…

"How else would you solve Ax = b for x, if not by computing the inverse of A (assuming A is a 3x3 or 4x4)?"

For instance, there a fast iterative solvers. They start with an initial guess for the solution (which can be an arbitrary vector) and successively generate more accurate approximations.

These iterative solvers don't need access to the inverse of A. Furthermore, many of those solvers don't even need access to A: They only need access to a subroutine which, given a vector u, returns Au.

Many matrices arising in applications (for instance one often reduces the problem of solving partial differential equations to the problem of solving a (very huge) linear system of equations) are sparse (have only few non-zero entries). For these it can make sense to write specialized code which specifically calculates Au.

If you want to read more about this, then a good starting point is https://en.wikipedia.org/wiki/Kaczmarz_method.

Post reply on HN