## Linear algebra in C++ is no fun

One of the reasons I'm in the UK this week is to work with David Balding on HapCluster, a method developed in his group that I have extended in various directions.

There's a few extensions we want to add, and maybe (finally) get a paper written about all the extensions.  The latter probably won't happen ... even with lots of extensions is is hard getting a paper accepted that just adds features here and there, and it is always necessary to compare with a lot of other tools (when really no single tool can beat all the others unless you find just the right "sweet spot" in parameter space where you are lucky, and that takes months of fiddling with parameters and I just can't be bothered).

Anyway, one of the extensions we are looking at is mapping of quantitative traits.  Ed Waldron already did this in his thesis, but I never got around to implementing it in my tool, 'cause I didn't fully understand the mathematics.  I still don't, but I think that I get enough to implement the formulas now at least.

The only problem is that it is a lot of matrix expressions that I need to implement, and that is just hell to do in C++.  At least if you, as I do, use GSL and CBLAS.

A simple expression such as

becomes

gsl_matrix *
calc_Vstar(const gsl_matrix *B, const gsl_matrix *invV)
{
gsl_matrix *copy = gsl_matrix_alloc(invV->size1, invV->size2);
gsl_matrix_memcpy(copy, invV);
gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, B, B, 1.0, copy);
gsl_matrix *Vstar = invert_matrix(copy);
gsl_matrix_free(copy);
return Vstar;
}

with a pre-computed .

Computing

we can re-use the computation from before, but it still ends up as this long function:

gsl_matrix *
calc_mstar(const gsl_matrix *Vstar, const gsl_matrix *invV,
const gsl_matrix *B, const gsl_matrix *m,
const gsl_matrix *y)
{
assert(Vstar->size1 == Vstar->size2);
int p = Vstar->size1;

// x := V^-1 m
gsl_matrix *x = gsl_matrix_alloc(p, 1);
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, invV, m, 0.0, x);

// z := B'Y
gsl_matrix *z = gsl_matrix_alloc(p, 1);
gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, B, y, 0.0, z);

// x := V^-1 m + B'Y = x + z
for (int i = 0; i < p; ++i)
gsl_matrix_set(x, i, 0, gsl_matrix_get(x, i, 0) + gsl_matrix_get(z, i, 0));

// m* := (V^-1 + B'B)^-1 (V^-1 m + B'Y) where V* = (V^-1 + B'B)^-1
gsl_matrix *mstar = gsl_matrix_alloc(p, 1);
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, Vstar, x, 0.0, mstar);

gsl_matrix_free(x);
gsl_matrix_free(z);

return mstar;
}

You have to break down the expressions into tiny little chunks, like you are writing assembly programs for arithmetic.  It is very frustrating.

It is low-level, of course, because the direct translation of matrix expressions to code is pretty inefficient, with lots of temporary matrices being created and deleted.  Still, that is pretty much how my code looks anyway, I am just manually doing the direct translation, 'cause I don't want to spend too much time trying to optimise the code before I know that it actually works (and is too slow to be worth optimising).

Template expressions, like in Blitz++, should be able to handle this, but I have never managed to learn that library.  I don't really deal with matrix expressions that often, so I haven't bothered.  If I end up having to do much more matrix hacking, I will.  Life is just too short to manually coding it up...

--
118-134=-16

Tags: ,

### 8 Responses to “Linear algebra in C++ is no fun”

1. mw Says:

That's one of the places where operator overloading actually makes sense. Never used a library that didn't use overloading for matrix (and vector) operations.

2. Thomas Mailund Says:

Damn, I could also have checked Boost... http://www.boost.org/doc/libs/1_35_0/libs/numeric/ublas/doc/overview.htm

Damn it! Always check Boost when coding C++!

3. Thomas Mailund Says:

mw: yes, you want overloading here! It is just a bit tricky, 'cause evaluation order and such can matter a great deal, and obvious things like just changing your indices in loops instead of first transposing a matrix and then multiplying just doesn't work out of the box with simple classes and overloading.

You really want template expressions to combine convenient notation with computational efficiency.

uBLAS, that I just linked to, gets most of this right, I just don't like that * seems to be component wise and not matrix multiplication. I don't want to write C = prod(A,B), I want to write C = A*B, damn it! I am more likely to use matrix multiplication than component wise multiplication, so the former and not the latter should use the operator.

But I can probably get used to that, if the alternative means implementing my own library :)

I have to do something about quadratic forms, though. I don't see any easy way of computing v' A v using this library... that still looks like two multiplication operations, but it shouldn't really be.

4. Me Says:

You really should have a look at eigen2!
http://eigen.tuxfamily.org/index.php?title=Main_Page

5. TomD Says:

I encountered the same issues when coding against the Intel MKL BLAS from C++, so I wrote a convenience wrapper here: http://wiki.github.com/tdoris/cppmkl which you may find useful. It doesn't use operator overloading for performance reasons (cf. boost ublas noalias issues), but it takes most of the hard work out of calling blas functions and provides an aligned allocator for speed+precision. If you're already familiar with the cblas then it might be what you need.

6. Thomas Mailund Says:

Thanks for the links, guys, I'll have a look!

I haven't used much linear algebra in C++ before (usually I have used R, Octave/Matlab or SciPy for that), but for this project I need it to get my MCMC fast enough.

Now, of course, I'm getting punished for lack of experience and not knowing which libraries are available, so the links are much appreciated!

7. amix Says:

The current languages have a really, really bad support for mathematical expressions. One of the new languages that challenges this is Fortress:
http://projectfortress.sun.com/Projects/Community

One of the designers of Fortress is Guy Steele (Lisp*, designer of Java). Apart from an unicode syntax that can be used for math expressions it also offers implicit parallelism and transactions.

Another thing is Mathematica - it's becoming a pretty powerful programming environment with Mathematica 7 - - and I could bet that if you do lots of math, then Mathematica is much faster and more correct than the average C++/Python library. For example, I found this post quite interesting:
* http://blog.wolfram.com/2007/09/25/arithmetic-is-hard-to-get-right/ (12 years spent on doing arithmetic right :))

8. Anton Says:

Cool!