overly-ambitious-isqo and the design of numerical codes

I have finally released the overly-ambitious-isqo project on github!

I wanted to call out two particular design concerns I had.

rich language

My first goal was to try very hard to build up the C language to very succinctly express the main algorithm in src/isqo_functor.cpp in extremely rich language. It seems like numerical code is typically implemented with loops like for (int i=0; i<N; i) and method calls like deltal(xhat, mu). I have found it much easier to reason and think deeply about codes like for (int primal_index=0; primal_index < num_primal; primal_index) and method calls like linear_model_reduction(penalty_iterate, penalty_parameter).

It’s a simple thing, but if you want to calculate the linear model reduction of the penalty iterate, spotting the bug between these lines:

linear_model_reduction(penalty_iterate, penalty_parameter)
linear_model_reduction(feasibility_iterate, penalty_parameter)

is much easier to spot than the same bug in the equivalent code

deltal(hatx, mu)
deltal(barx, mu)

since the later requires remembering that hat represents “penalty” and bar represents “feasibility.” This kind of code is not AT ALL common in numerical codes, and I think that is a big huge bummer, because the code is incredibly difficult to understand, and it doesn’t need to be.

The other advantage of building the interfaces this way is that it is trivial to later add performance “tricks” like caching the last linear model reduction for a particular iterate, since the entire class hierarchy is set up so that such a mechanism could be used in common for all of the computable functions we have added. Pretty neat, right?

strong types

The other big goal was to try to have extremely strong types in the code. Typical numerical codes have functions which accepts 20 pointers to double arrays, each of which is critical for the computation.

This seems, to say the least, completely insane to me! There are so many errors that the compiler could catch, simply by defining them as a structure and passing pointers/references to those instead! And it enforces a very clear standard for otherwise arbitrary things. The most obvious convention is simply the order in which they’re called–why should a programmer have to remember or look up 20 items (all of identical types, mind you!) when really they’re just trying to pass in one particular structure? It’s not just that, though. Another example: Some drafts of our document included a matrix equality constraint like

 A*x = b

but some external libraries we linked to specified the same constraint as

 A*x + c = 0.

The original code for the algorithm used the former convention, which was an endless source of bugs and countless, “multiply every element by -1 here” type code. There’s no reason for it! It’s much simpler to define a type

• EqualityConstraint - a constraint of the type A*x=b
• Matrix A - stores values of matrix
• Value b - stores values of the “right-hand side” of the equality constraint

Don’t you agree?