Tuesday, December 20, 2011

Numerically stable standard deviation calculation and code perforation

How would you implement a class which has an append(double x) function to collect values, and a get_std_dev() function which returns the standard deviation of the values collected?

Sounds like an easy problem, here is my guess at what your code would look like:

This looks like an standard implementation, but it contains some subtle bugs. These bugs do not show up in regular usage, but they lie patiently in wait. Remember the tricky binary search implementation bug which escaped detection in the JDK for about a decade and is now staple fodder of all technical interviews?
I recently came across Knuth's numerically stable standard deviation calculation algorithm on the NS3 developer's mailing list. That algorithm looks like this:

 (Implementation borrowed from Mitch's version for NS3.)

To see the difference between the two implementations, try this:

For the binary search algorithm, the fix was relatively simple: change (high + low) / 2 with low + (high - low) / 2; this prevented the potential integer overflow. The changes required in this case are similar, but for floating point overflow. They are more subtle but less disruptive since getting an inf is more natural than getting negative array index values. This might not have been so in Knuth's time.

This is a humbling lesson again, a reminder that even simple code snippets can be difficult to get absolutely right and our world works on tech which is be too dizzying to be understood by one mind. Each engineer plays his part in the affair, does his best, and hopes that others do the same. And, amazingly, it all works. Well, most of the time anyway.

Interestingly, there is a down-side to doing the best job without having the complete picture in mind: one will implement a more complex algorithm to ensure correctness where perhaps speed was the primary requirement. This happens so often that Prof. Rinard at CSAIL-MIT is investigating code-perforation as a means to hasten execution. In layman terms, code-perforation works by executing fewer iterations of each loop in the code and seeing how much it effects the speed and accuracy of the results. Here, the line between correctness and accuracy is fairly thin and the natural instinct of the developers is to make each function, each loop, each segment as accurate as possible. In such settings, code perforation manages to relax those loops whose good hard and accurate work is undone by a later loop which produces less accurate estimates. This works well in applications such as video encoding/decoding and Monte Carlo simulations.

I think that errors like this numeric instability of Standard Deviation would be much harder to catch (or notice?) after perforation. However, that is stuff for another post.

No comments: