Comment by albertzeyer | Hacker News Reader

Comment by albertzeyer | Hacker News Reader

Parent Full thread

I wonder if it was all handwritten. Esp. functions like invertMatrix or premultiplyMatrix. Or what he used to generate the code.

I once used WolframAlpha to get me some complicated formulas and wrote a custom WA-plaintext-output-to-C translator for that output to use it in my code. Here is the script: https://github.com/albertz/helpers/blob/master/wolframalpha_...

Example input: http://www.wolframalpha.com/input/?i=solve+a*x_1%5E3+%2B+b*x...

Output: s->a = ((x1 + x2 - 2.) / pow(x2 - x1, 3.)); s->b = ((- (((x1 + x2 - 2.) * pow(x1, 2.)) / pow(x2 - x1, 3.)) - ((4. * x2 * (x1 + x2 - 2.) * x1) / pow(x2 - x1, 3.)) + ((6. * (x1 + x2 - 2.) * x1) / pow(x2 - x1, 3.)) - ((7. * pow(x2, 2.) * (x1 + x2 - 2.)) / pow(x2 - x1, 3.)) + ((6. * x2 * (x1 + x2 - 2.)) / pow(x2 - x1, 3.)) - 1.) / (4. * x2 - 4.)); s->c = (1. / 2.) * ((((x1 + x2 - 2.) * pow(x1, 2.)) / pow(x2 - x1, 3.)) + ((4. * x2 * (x1 + x2 - 2.) * x1) / pow(x2 - x1, 3.)) - ((6. * (x1 + x2 - 2.) * x1) / pow(x2 - x1, 3.)) + ((pow(x2, 2.) * (x1 + x2 - 2.)) / pow(x2 - x1, 3.)) - ((6. * x2 * (x1 + x2 - 2.)) / pow(x2 - x1, 3.)) - ((4. * (- (((x1 + x2 - 2.) * pow(x1, 2.)) / pow(x2 - x1, 3.)) - ((4. * x2 * (x1 + x2 - 2.) * x1) / pow(x2 - x1, 3.)) + ((6. * (x1 + x2 - 2.) * x1) / pow(x2 - x1, 3.)) - ((7. * pow(x2, 2.) * (x1 + x2 - 2.)) / pow(x2 - x1, 3.)) + ((6. * x2 * (x1 + x2 - 2.)) / pow(x2 - x1, 3.)) - 1.)) / (4. * x2 - 4.)) + 1.); s->d = (1. / 4.) * ((((x1 + x2 - 2.) * pow(x1, 3.)) / pow(x2 - x1, 3.)) - ((4. * x2 * (x1 + x2 - 2.) * pow(x1, 2.)) / pow(x2 - x1, 3.)) - (((x1 + x2 - 2.) * pow(x1, 2.)) / pow(x2 - x1, 3.)) - ((pow(x2, 2.) * (x1 + x2 - 2.) * x1) / pow(x2 - x1, 3.)) + ((2. * x2 * (x1 + x2 - 2.) * x1) / pow(x2 - x1, 3.)) + ((6. * (x1 + x2 - 2.) * x1) / pow(x2 - x1, 3.)) + x1 - ((pow(x2, 2.) * (x1 + x2 - 2.)) / pow(x2 - x1, 3.)) + ((6. * x2 * (x1 + x2 - 2.)) / pow(x2 - x1, 3.)) + ((4. * (- (((x1 + x2 - 2.) * pow(x1, 2.)) / pow(x2 - x1, 3.)) - ((4. * x2 * (x1 + x2 - 2.) * x1) / pow(x2 - x1, 3.)) + ((6. * (x1 + x2 - 2.) * x1) / pow(x2 - x1, 3.)) - ((7. * pow(x2, 2.) * (x1 + x2 - 2.)) / pow(x2 - x1, 3.)) + ((6. * x2 * (x1 + x2 - 2.)) / pow(x2 - x1, 3.)) - 1.)) / (4. * x2 - 4.)) + 1.);

(https://github.com/albertz/music-player/blob/master/ffmpeg_p... SmoothClipCalc::setX)

It's worth pointing out that you should rarely, if ever, use exact formulas since there are a lot of problems with numerical stability amongst other things. Often it is faster to just numerically solve the system.

For example, in your case you are solving Ax=b which you can easily solve by a simple LU decomposition. Plus if your A changes, your exact formula is wrong, so using an LU decomposition just makes sense. Plus your code will look way cleaner as well. And way easier to debug.

http://en.wikipedia.org/wiki/LU_decomposition

His system of equations is not linear so vanilla Ax=b won't work.

It is actually linear. a,b,c,d are the unknowns and x1 and x2 are the inputs.

A bit confusing since normally x1,x2 denote the unknowns.

Not if you are doing econometrics/statistics.

Unless you like econometrics: then X are data, betas and gammas and sigmas are unknown (oh my!).

Jargon is fun, unless I want to communicate across fields :/.

Well today I learned. I never knew that. Stupid mathematics and its inconsistent notation.

It makes sense when you look at historical development, but there is no reason to learn separate notations now without learning all sides :).

He's solving for abcd, which is linear in powers of x. I'm not going to argue for what's the correct solution, but it's certainly not that mess.

In case that the exact formulas for this solution would be numerically unstable, why would the LU decomposition improve the numerical stability? Also, why should it be faster?

Well, to be fair, in that example, if the compiler does those `pow(x2 - x1, 3.)` calls multiple times, this would not be optimal, but otherwise, it should be ok.

I did this because I will probably never ever change the matrix A here and I wanted to make that code very fast. Otherwise, it's of course a good idea to use some more generic solution.

I'm not implying yours are not unstable, it's just there is a chance it might be. Typically exact formulas exhibit such properties.

LU decomposition is numerically stable, for the most part. You can run into problems with floating point arithmetic if you're not careful. This just has to do with the way you solve the problem. By breaking up the matrix into an upper and lower triangular part and solving the resulting triangular system, you often avoid things like powers or square roots and just reduce everything to basic addition and multiplication which tend to have more stable numerical properties.

As for speed, well in this small case of N=4 there is probably very little speed increase since LU is O(N^{3}), although this can be improved depending on symmetries of the matrix. Maybe for N=5 I might write out the exact formulas but beyond that I would just a LU solver since the amount of code for an LU solver is fairly compact.

But, as you correctly point out, if you are never changing the matrix A then you are probably safe with writing it out this way.

I just come from a numerical fluids background so seeing exact formulas makes me uneasy. And in my work since the number of grid points N tends to be variable so general formulas are not available. Plus I've noticed a lot of programmers tend to not know numerical algorithms. In fairness, I don't blame them. All the numerical computations courses I've taken and TAed are very boring and don't make you do anything fun and numerical analysis has this reputation of being dry. If you actually made them write stuff like a Navier-Stokes solver, you would get them interested.

By the way if you want to learn more about numerical linear algebra, which is the cornerstone of most scientific and high performance computing, I personally enjoy Trefethen and Bau's book. Although it's aimed at a mathematical audience and assumes such.

Wait, no, LU decomposition is not stable unless you're doing pivoting. And if you're doing pivoting, then it very much stops looking like an "exact formula", because it'll get filled up by ifs and variable renamings, etc.

There's plenty of cases in graphics and simulation where solving even a 4x4 matrix needs careful numeric consideration (see this paper about tetrahedral mesh simplification, section III.D, page 7: http://www.sci.utah.edu/~hvo/papers/tetstream.pdf)

If you care about accuracy at all, do not use the straightforward "exact formulas" like Cramer's rule. That's just asking for trouble.

And for numerical linear algebra, I'd go start with Strang and chapter 9: http://math.mit.edu/linearalgebra/

You are quite correct about stability but I didn't want to get into technical details. And indeed you bring up correctly that such numerical issues mustn’t be treated lightly. I don't do computer graphics and the matrices I deal with are well conditioned such that problems of stability that arise due to linear algebra routines are negligible. In my simulations, instability arises due to the wrong time-stepping scheme, which is a completely different issue.

And good recommendation by Strang. I believe he has some fantastic MIT OpenCourse lectures on numerical linear algebra as well.

The conventional wisdom (in scientific computing) that you seldom need to compute an inverse of a matrix does not apply to computer graphics. 4x4 matrix inverse is one of the most useful functions in computer graphics and is used a lot.

It is most often used as a hand rolled 4x4 matrix inversion function, not a generic NxN matrix inversion.

> I wonder if it was all handwritten. Esp. functions like invertMatrix or premultiplyMatrix. Or what he used to generate the code.

Yes, it probably is hand written. Or more likely, copied from some well known resource such as MESA.

Computer graphics deals extensively with 4x4 matrices, so there are hand written implementations of elementary operations such as matrix inverse, all unrolled and precomputed so that there's no loops or anything. (NOTE: this isn't really basic loop unrolling, these are not simple loops to begin with).

http://stackoverflow.com/questions/1148309/inverting-a-4x4-m... http://stackoverflow.com/questions/2624422/efficient-4x4-mat...

For cases where the most basic simplifications are enough (0*x=0 etc.) there is this one: http://glat.info/flatorize/

Reply on news.ycombinator.com