[Rcpp-devel] C++ design for least squares solvers

[Douglas Bates]

A fundamental operation in the work that I do is to solve least
squares problems in various configurations.  The simplest case is to
solve a least squares problem with a dense model matrix and a numeric
response, but there are many other cases based on properties of the
model matrices.

For the model matrices I can have one or two matrices and the matrices
can be sparse or dense.  For the types of problems I may want to
incorporate weights and/or a penalty term.  These properties are
somewhat orthogonal.

The underlying code for the solution is very different for the
different classes but the types of operations (solve the system for a
given response, evaluate the fitted values, ...) are common.

Coming from the C/R programming world I initially approached this by
defining a base class and defining classes that specialized this
class.  The functions to do the actual operations were virtual
functions but that approach seems to get very complicated, perhaps
because I haven't thought it through completely yet.  It seems that
templates would be a more effective way of approaching this but I am
still feeling my way around the corners of the template idiom.

So I am trying to decide how best to express such ideas in C++ code.
Ultimately I want to combine various forms of least squares solvers
with the idea of a generalized linear model family so that I can write
an iteratively reweighted least squares algorithm at a high level and
have the classes take care of the fussy details.

I realize that this may not be enough information to offer an opinion.
 If desired I can make things more explicit by describing exactly
which operations should be incorporated in a least squares solver and
some of the possible combinations.