Skip to content

Finish your derivations, please

October 21, 2010

Every time you ship a product with half-assed math in it, God kills a fluffy kitten by feeding it to an ill-tempered panda bear (don’t ask me why – I think it’s bizarre and oddly specific, but I have no say in the matter). There’s tons of ways to make your math way more complicated (and expensive) than it needs to be, but most of the time it’s the same few common mistakes repeated ad nauseam. Here’s a small checklist of things to look out for that will make your life easier and your code better:

  • Symmetries. If your problem has some obvious symmetry, you can usually exploit it in the solution. If it has radial symmetry around some point, move that point to the origin. If there’s some coordinate system where the constraints (or the problem statement) get a lot simpler, try solving the problem in that coordinate system. This isn’t guaranteed to win you anything, but if you haven’t checked, you should – if symmetry leads to a solution, the solutions are usually very nice, clean and efficient.
  • Geometry. If your problem is geometrical, draw a picture first, even if you know how to solve it algebraically. Approaches that use the geometry of the problem can often make use of symmetries that aren’t obvious when you write it in equations. More importantly, in geometric derivations, most of the quantities you compute actually have geometrical meaning (points, distances, ratios of lengths, etc). Very useful when debugging to get a quick sanity check. In contrast, intermediate values in the middle of algebraic manipulations rarely have any meaning within the context of the problem – you have to treat the solver essentially as a black box.
  • Angles. Avoid them. They’re rarely what you actually want, they tend to introduce a lot of trigonometric functions into the code, you suddenly need to worry about parametrization artifacts (e.g. dealing with wraparound), and code using angles is generally harder to read/understand/debug than equivalent code using vectors (and slower, too).
  • Absolute Angles. Particularly, never never ever use angles relative to some arbitrary absolute coordinate system. They will wrap around to negative at some point, and suddenly something breaks somewhere and nobody knows why. And if you’re about to introduce some arbitrary coordinate system just to determine some angles, stop and think very hard if that’s really a good idea. (If, upon reflection, you’re still undecided, this website has your answer).
  • Did I mention angles? There’s one particular case of angle-mania that really pisses me off: Using inverse trigonometric functions immediately followed by sin / cos. atan2 / sin / cos: World’s most expensive 2D vector normalize. Using acos on the result of a dot product just to get the corresponding sin/tan? Time to brush up your trigonometric identities. A particularly bad offender can be found here – the relevant section from the simplified shader is this:
    float alpha = max( acos( dot( v, n ) ), acos( dot( l, n ) ) );
    float beta  = min( acos( dot( v, n ) ), acos( dot( l, n ) ) );
    C = sin(alpha) * tan(beta);

    Ouch! If you use some trig identities and the fact that acos is monotonically decreasing over its domain, this reduces to:

    float vdotn = dot(v, n);
    float ldotn = dot(l, n);
    C = sqrt((1.0 - vdotn*vdotn) * (1.0 - ldotn*ldotn))
      / max(vdotn, ldotn);

    ..and suddenly there’s no need to use a lookup texture anymore (and by the way, this has way higher accuracy too). Come on, people! You don’t need to derive it by hand (although that’s not hard either), you don’t need to buy some formula collection, it’s all on Wikipedia – spend the two minutes and look it up!

  • Elementary linear algebra. If you build some matrix by concatenating several transforms and it’s a performance bottleneck, don’t get all SIMD on it, do the obvious thing first: do the matrix multiply symbolically and generate the result directly instead of doing the matrix multiplies every time. (But state clearly in comments which transforms you multiplied together in what order to get your result or suffer the righteous wrath of the next person to touch that code). Don’t invert matrices that you know are orthogonal, just use the transpose! Don’t use 4×4 matrices everywhere when all your transforms (except for the projection matrix) are affine. It’s not rocket science.
  • Unnecessary numerical differentiation. Numerical differentiation is numerically unstable, and notoriously difficult to get robust. It’s also often completely unnecessary. If you’re dealing with analytically defined functions, compute the derivative directly – no robustness issues, and it’s usually faster too (…but remember the chain rule if you warp your parameter on the way in).

Short version: Don’t just stick with the first implementation that works (usually barely). Once you have a solution, at least spend 5 minutes looking over it and check if you missed any obvious simplifications. If you won’t do it by yourself, think of the kittens!

About these ads

From → Maths

9 Comments
  1. Won permalink

    Amen! Particularly about avoiding trig.

  2. Justin Paver permalink

    Nice post!

    > Don’t invert matrices that you know are orthogonal, just use the transpose!

    I believe the basis vectors should have unit length too for inverse==transpose, so really the transpose should only be used if the matrices are orthonormal.

    > World’s most expensive 2D vector normalize.

    No, *this* is the World’s most expensive 2D vector normalize :)

    void vec2dNormalize(vec2d& outNormal, const vec2d &inNormal)
    {
    while (1)
    {
    outNormal.x = frand(-1.0f, 1.0f);
    outNormal.y = frand(-1.0f, 1.0f);
    // yeah I said it. no epsilon thresholds
    bUnitLen = (outNormal.x*outNormal.x + outNormal.y*outNormal.y) == 1.0f;
    bParallel = (-outNormal.y*inNormal.x + outNormal.x*inNormal.y) == 0.0f;
    if (bUnitLen && bParallel)
    return;
    }
    }

    Possibly so expensive it’ll never finish. Sadly, I’ve actually seen someone sort lists like this too. Randomly permuting the list, and then testing whether it’s sorted until they get a sorted list :(

    • “I believe the basis vectors should have unit length too for inverse==transpose, so really the transpose should only be used if the matrices are orthonormal.”
      The terminology is a bit unfortunate there. Matrices are called orthogonal if and only if the corresponding basis vectors are orthonormal (the usual definition is that Q^T Q = Q Q^T I). There’s no special name for matrices with orthogonal but non-normalized basis vectors.

  3. Per Vognsen permalink

    Numerics checklist:

    – Special cases are ugly but necessary in most numerical code. Be suspicious of too beautiful code that looks like a one-to-one translation from a paper or textbook. Most likely it is full of pitfalls.
    – As a rule, prefer self-correcting iterative methods over direct methods. Even a venerable direct method like Gaussian elimination needs iterative refinement for numerical robustness in the face of ill-conditioned problems. This also applies to simple closed forms like the cubic formula. The chief virtue of self-correction is that it makes intermediate round-off error almost a non-issue. But you must still think carefully about your stopping condition (error vs residual conditioning, circle of uncertainty in calculating error and residual deltas, etc).

  4. royalestel permalink

    I hate to look unsmart, but I’m working on this:
    Ouch! If you use some trig identities and the fact that acos is monotonically decreasing over its domain, this reduces to:

    Think you could expand that a little?

    • Arc cosine is monotonically decreasing: see the diagram here (blue curve). Note how the curve strictly goes down as you go further to the right. This is what “monotonically decreasing” means, visually.

      In formulas, this means that if x acos(y): x is further to the left, so its arc cosine value is higher, because the arc cosine is steadily decreasing as you go to the right!

      Hence alpha = max( acos( dot( v, n ) ), acos( dot( l, n ) ) ) is the same as alpha = acos( min( dot( v, n ), dot( l, n ) ) – max turns into min because smaller values of x imply larger values of acos(x). Similar for beta.

      Now, the only thing we do with alpha is compute its sine. By the Pythagorean theorem, sin(x)^2 + cos(x)^2 = 1. So sin(alpha) ^ 2 + cos(alpha) ^ 2 = 1, which in turn means that sin(alpha)^2 = 1 - cos(alpha)^2, so sin(alpha) = +- sqrt(1 - cos(alpha)^2). (Positive sign turns out to be correct in this case). So we never need to compute the acos for alpha in the first place: we already have cos_alpha, which is just cos_alpha = min( dot( v, n ), dot( l, n ) ).

      tan(beta) is similar: tan(beta) = sin(beta) / cos(beta) = sqrt(1 - cos_beta^2) / cos_beta. So C = sin_alpha * tan_beta = sqrt(1 - cos_alpha^2) * sqrt(1 - cos_beta^2) / cos_beta = sqrt((1 - cos_alpha^2) * (1 - cos_beta^2)) / cos_beta.

      Now note that the square root term contains the same expressions for both cos_alpha and cos_beta and just multiplies them together – so it doesn’t matter which of them is the min of the two dot products and which is the max; we calculate the same value for both dot products and multiply them together. Introducing some more names:

      vdotn = dot( v, n )
      ldotn = dot( l, n )
      cos_alpha = min( vdotn, ldotn )
      cos_beta = max( vdotn, ldotn )

      so

      C = sqrt((1 - cos_alpha^2) * (1 - cos_beta^2)) / cos_beta = sqrt( ( 1 - vdotn^2 ) * ( 1 - vdotl^2 ) ) / max( vdotn, ldotn ).

      • royalestel permalink

        Thank you so much for taking the time to spell this out for me! I get it, yet at the same time I know it’s just convoluted enough that I wouldn’t arrive at the solution. After reading your blog and Emil Persson’s “Low Level Thinking” presentation at GDC’13, it made me look at some of the shader code we have here which is ripe for optimization.

        The fact that this is a pet peeve of yours and his, means that you both see this as obviously important and fairly trivial to do.

        My thinking is instead that the both of you had lots of practice optimizing code using trig identities and many other methods, that are not even a consideration when teaching programming.

        Your combined works imply is that there is a strong need for practice optimization problems, at the algorithmic level and low level for professionals.

        In the long term, University course work ought to involve lots of these sorts of optimizations, particularly involving dot products and trig identities.

        I’ll keep plugging away at my code optimizations (I know I will get rid of that atan somehow), but thanks so much for showing me how a pro does it.

Trackbacks & Pingbacks

  1. (Linear) Algebra Toolbox 1 « The ryg blog
  2. Frustum culling: turning the crank « The ryg blog

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 207 other followers

%d bloggers like this: