Live data from Hacker News

Don’t invert that matrix (2010)

johndcook.com

41–48 of 48 posts

Re: Don’t invert that matrix (2010)

#41
Sparse matrix factorization is a black art and is best done with some existing out-of-the-box solution.

The problem is that they're all these FORTRAN packages that one person wrote, and then everybody, in the name of human progress, reused and reused ever since -- akin to bootstrapping human intelligence.

Nobody ever wants to write the Lanczos bidiagonalization algorithm with partial reorthogonalization from scratch, and when I considered writing my own C++ implementation to bring human progress to the 21st century, my academic advisor told me "just don't."

EDIT: The real reason is that FORTRAN is less apt to copy and more apt to reference than C++.

Re: Don’t invert that matrix (2010)

#42
post #9

Earlier quoted context omitted.

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

This. Or you know how to calculate the closed-form inverse a priori like with projection matrices (e.g. gluUnProject).

Re: Don’t invert that matrix (2010)

#43
post #31

Earlier quoted context omitted.

In games you're probably applying many small inverse matrices to many small vectors. In partial differential equations, you solve a single large linear system where the unknown vector represents a field defined at every point on a spatial mesh. In my research, specifically, that spatial mesh can have as many as 10^7 grid points. For linear systems of this size, any exact solution at all is impractical - we exclusivel…

And not only are the matrices small, they're probably known (or assumed) to be well-conditioned. Inverting even a small hilbert matrix will probably go boom.

Not just well conditioned, but special orthogonal.

Re: Don’t invert that matrix (2010)

#44
post #27

It's funny because the whole discussion only exists because there is no division operator in "matrix math". For scalars we can simply write a / b (instead of a^-1 * b) and everyone knows how to calculate it. But for some reason there's no such thing when dealing with matrices. Well it turns out there is something like division for matrices and it's called factorization (followed by back-substitution) and there are ac…

It's a little more complicated than: """It's funny because the whole discussion only exists because there is no division operator in "matrix math"."""

Matrices in general (insert caveats on size here) form a ring, but not a division ring (https://en.wikipedia.org/wiki/Division_ring). This means that not all (square) matrices will have a multiplicative inverse. This is why we often restrict to things like GL(n), the group of non-singular square matrices of size n.

Furthermore when you say "scalars", they have multiplicative inverses because they're elements of a field (and sometimes the inverse might be a pain like if you're in a number field or something).

This isn't even getting into the Moore-Penrose pseudoinverse (https://en.wikipedia.org/wiki/Moore–Penrose_pseudoinverse) or linear operators or so on, but I'll let someone else shout about algebra.

Re: Don’t invert that matrix (2010)

#45

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 statist…

[deleted]

Re: Don’t invert that matrix (2010)

#46

Earlier quoted context omitted.

You calculate the factorization (QR, LU, or similar), but you never calculate the inverse of the factorization: A x = b A => QR # QR decomposition Q R x = b R x = Q' b # Q is inverted by transposing x = backsubstitute(R, Q' b) Q was the transpose of it's inverse, and so you don't get any rounding error when you transpose it. R didn't need to be inverted because it was in a form where back substitution lets you apply…

MATLAB's linsolve uses your strategy! ( http://in.mathworks.com/help/matlab/ref/linsolve.html )

Linear models in R are also solved by QR factorization.

Re: Don’t invert that matrix (2010)

#47
post #3
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.

Another relevant difference is that precision is not nearly as critical in games as in most numerical applications. Errors of a few percent will likely go unnoticed. Graphics implementations therefore usually only use single precision.

Re: Don’t invert that matrix (2010)

#48
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 t…

I avoid floating point wherever possible (especially after reading "What Every Programmer Needs To Know About Floating Point") and just use large integers. Not sure if the particular math in question here can be modified to use large integers...
Post reply on HN