Commenting code may lead to category theory

Not that category theory is hazardous to your health, of course!

This topic came up when i was thinking about how to justify the deficiencies of a particular C++ interface i had written. These functions take the name of a file containing some kind of matrix data. They then construct and return the corresponding matrix object, possibly distributing that object over some number of MPI processes. The trick is that a file could contain all kinds of different data -- both different kinds of matrices (sparse or dense, symmetric or not) and different kinds of data (integer, real, or complex). The file begins with a header containing all of this "metadata" that describes its contents, but you have to open up the file and start reading it in order to find that out.

The issue is that the natural interface calls for each of these functions to return an object of a particular type. The type of object returned by the function is therefore fixed, even though the contents of the file are unknown. There are two reasons for this. First, different kinds of matrices (e.g., sparse or dense) require different kinds of metadata as input, before reading the file at all, in order to interpret the file's contents. (For example, large square sparse matrices may be distributed in a 2-D fashion over processors, so they need both row and column distribution information. Large, dense vectors may only be distributed by rows over processors, so they only need row distribution information.) Second, even if we fix what kind of matrix we want in the interface, we still don't know the type of data in the matrix. It could contain integers, real numbers of various precisions, complex numbers, or even more complicated objects (quotients of polynomials with coefficients in a finite field, for example).

Concretely, in C++ terms, i might have a SparseMatrix or DenseMatrix type. Each of these is "templated" (parametrized) on the type of entries in the matrix: SparseMatrix<double> for a sparse matrix with double-precision real values, SparseMatrix<int> for a sparse matrix with 32-bit signed integer values, etc. This lets us implement (for example) SparseMatrix once for entries of any type. It also lets me write the file reader function once for any parametrization of SparseMatrix. All these are good. The problem is that the user has to pick a specific parametrization of SparseMatrix (e.g., SparseMatrix<double>) before calling the function.

Now, there are ways of implementing types that could be any of various types. Really, what you want are "algebraic data types" that are a disjoint union of different types (for example, "int or double or complex"). C++ doesn't have those, strictly speaking, although you can fake them with so-called "variant" types. Other languages, like Haskell, do have algebraic data types.

Suppose, though, that C++ did have algebraic data types. I will represent the disjoint union here nonstandardly with a plus sign: e.g., "int + double + complex<double>." What i might like the file reader to do is to return a sparse matrix, whose entries are int + double + complex<double>. In our notation, this would be called a "SparseMatrix<int + double + complex<double> >." However, C++ makes this undesirable. This is because in C++, the bitwise arrangement of data in memory matters a lot. Implementing a "variant" type in C++ that can hold any of int, double, or complex<double> would require at least as much space as that required by a complex<double> -- twice as much as a double, and twice as much again as an int. Furthermore, i really don't want a sparse matrix for which different entries might have different types -- all of the entries should have the same type. Yet, i'm forced to do this if the function returns a SparseMatrix<int + double + complex<double> >.

The function might instead return a SparseMatrix<int> + SparseMatrix<double> + SparseMatrix<complex<double> >. However, now when i get the return value, i have to guess what type the object is, trying each of the possibilities until i find out. This is easy enough for the three possibilities here. However, suppose that SparseMatrix takes several type parameters (in our case, it takes five). Do i have to guess and check for all possible combinations?

It seems natural that SparseMatrix<int> + SparseMatrix<double> should be the same type as SparseMatrix<int + double>. However, i was explaining this problem to a colleague tonight, and it occurred to me that templated classes in C++ are just functions between types. For example, SparseMatrix<T> defines a function from the entry type T to the specific sparse matrix type SparseMatrix<T>. If f is a function, why should i expect that f(x + y) = f(x) + f(y)? That would make the function f only one step away from linear -- a very special function indeed! Similarly, i should expect that a SparseMatrix<int> + SparseMatrix<double> would be a very different sort of object than a SparseMatrix<int + double>. This makes sense if you think about possible representations of each of these objects in memory: the former might be a pointer to some kind of matrix, whose arrays of matrix entries are either arrays of int or arrays of doubles, while the latter would be a matrix whose arrays of matrix entries are arrays of a variant type.

I'm not a category theorist, but thinking about this problem made me think about type systems and algebras over type systems. So, i suppose it's a little bit of category theory, no? All i have to do is not embarrass myself too much in front of a real category theorist!

Scientists "vs." programmers: It's all about reproducibility

I was troubled by Dr. John Cook's recent post on the different programming styles of scientists and programmers.  It felt like an exaggerated dichotomy to me: the "cowboy" scientist who tosses off throw-away scripts, and the cube-dwelling programmer implementing a serious product. Dr. Cook was trying hard to think of this dichotomy as a cultural distance, based on his experience working with and trying to understand the needs of both groups.  However, I think he missed the point that joins scientists and programmers:  reproducibility. 

As the article mentions, programmers understand that their code needs to produce predictable results.  They set up tests to verify this.  Entire fields of computer science research revolve around proving code correctness.  As a developer of numerical algorithms, I think of my code as an experimental apparatus for testing my hypothesis about the superiority of one algorithm or optimization over another.  I can't say I've always done my best to make my code work that way, but I'm working on it.

I think many scientists haven't yet learned that code is really just experimental apparatus, and the reason for this in part may be that modern scientists are not rewarded for reproducing "known" results.  Programmers often are -- why else did Google+ follow Facebook, which in turn followed MySpace?  Now, a lot of code really is throw-away analysis scripts.  This is probably what Dr. Cook means by the following:  "Programmers need to understand that sometimes a program really only needs to run once, on one set of input, with expert supervision." However, simulations are not like this.  A simulation code is an experimental apparatus.  It need not be permanent, but it should be constructed well enough to make its results trustworthy.  Just as importantly, it should be described and documented sufficiently well that other scientists can replicate its results. 

The real problem is not a cultural difference between scientists and programmers: the real problem is reproducibility.

A mathematical chauvinist?

An afternoon discussion at work among some of our interns revolved around a question that one (A) posed to another (B): "How many molecules of Abraham Lincoln are in a glass of water?"  After sufficient simplifying assumptions, B (a mathematician) was able to solve the problem with a simple computation.  The process resulted in a debate between A and B over the virtues of being able to derive formulae and constants, rather than looking them up with a search engine.  I sided with B, of course! as having and exercising the ability to carry out that process is much more useful than demanding that a search engine do all the work.  Simple units and order-of-magnitude calculations are things my PhD advisor does very rapidly, and when I was his student, I always found myself lagging behind and feeling not so intelligent.  Exercising that skill at least saved me a little embarrassment!

Later that evening, A was explaining to me how one writes applications for a cell phone.  I found myself feeling perfectly happy that I didn't have to do all that work.  Cell phones are things I want to use to make phone calls, and maybe look up directions in an emergency; I struggled to consider why I would want to write my own app.  However, then I realized I was doing the same thing that A had done above: not valuing the process over the answer.  Learning something about writing a cell phone application might serve me well at some future time, even though I can't imagine the value of this time investment now.  I suppose I was being a mathematical chauvinist!

Git: recovering from an incorrect e-mail address

Our big software project uses Git (the version control system) in a CVS mode, with a single central repository to which we push commits (after they pass the check-in test suite).  This can be good, because it forces us to fix small but annoying things -- like setting the correct e-mail address!  I was testing on a different machine than my usual development workstation; I made some commits there, and then pushed to my workstation's local repository (Git shines for stuff like this).  Then, I activated the check-in tests, but the push failed: my e-mail address wasn't set!  (The tests check for that -- handy!)

While I had set my e-mail address for Git correctly on my workstation, I hadn't set it on the other machine.  There were seven commits sitting on my workstation (where I run the check-in test suite) with the wrong e-mail address.  Here's how I started the fix.  First, I ran

$ git rebase -i HEAD~7

which let me fix up the last 7 local commits.  "-i" means "interactive," so Git fired up a text editor and let me decide which commits to change (and how).  I then repeated the following three tasks seven times:

$ git commit --amend --reset-author
(Edit the commit message -- no author e-mail address here, but it gets fixed due to "--reset-author")
$ git rebase --continue

That fixed it!  There's probably an easier way, but I was pleased.

Larval sparse matrices?

John Rose's post on "larval objects in the JVM" suggests a type distinction between "larval" and "adult" objects.  Larval objects are mutable, suitable for in-place updates and for use as scratch computation, but are not safe for sharing by multiple threads in the same shared-memory space.  Adult objects are immutable and safe to share in parallel, but "changing" them requires copying, which can be expensive.  His examples mainly involve small objects, like complex numbers or dates.  However, he does suggest a "large-scale" example of a database, where the default state is immutable ("adult"), but a mutable API can be made available by request.

This pattern of a "larval" initialization stage -- too complex and incremental to handle in a constructor -- followed by a "publish" command (we call it "fillComplete") that transforms the data structure into its "adult" state -- is more or less how sparse matrices are used in practice.  The data for the sparse matrix is processed in "element" form (more or less equation by equation), and then "assembled" into the optimized data structure.  Changing the optimized final data structure can be expensive and may require re-optimization.  Thus, we like to think of it as immutable when optimized.

The problem is, the sparse matrix API contains operations for the "larval" stage, which are either invalid or expensive when the matrix is in its "adult" stage.  Rather than discouraging users from doing the wrong thing, it would be nice for the API to forbid invalid operations syntactically (not by throwing an exception and frustrating new users with a mysterious error message), and also make clear which operations are allowed but undesirable due to their cost.  Making it syntactically harder to do suboptimal things is a good way to intimidate nonexpert programmers from doing them!

John mentions the library solution to this problem: "builder" objects or functions, like Java's StringBuilder.  They take as input some number of updates, applied incrementally.  When the user signals that all updates are done, they produce a "final" immutable object.  In the case of sparse matrices, that would mean users can only call mutating methods on the "SparseMatrixBuilder" object.  To mutate a "finished" sparse matrix, they would have to request a "builder."


The problem with this approach is the ability to request a builder, after the original completion of the immutable object.  It means that immutable objects can change state.  It also means that the immutable object (which is now only immutable in terms of its API) and the mutable "view" coinhabit the same scope.  Now the state of the "immutable" object depends on the semantics of view updates.  Do they happen as soon as the update methods are called?  (Simultaneity doesn't mean much in the context of multiple threads.)  Do they get accumulated and applied in batches?  Do they involve an asynchronous call to some remote compute device, like a GPU?  Programmers shouldn't be depending on these sorts of semantics.  It's better to avoid the mess and forbid changing immutable objects back into mutable ones, or return to the original solution and provide only a single object type with immutable and mutable states.


What's missing here is a "rebinding syntax" that prevents programmers from accessing the supposedly immutable object while it's in a mutable state.  This is a good use case for the "WITH-" Lisp idiom, or for Python's context guards:  the mutable view created by the "with" isn't allowed to escape its scope.

Using ScaLAPACK: not so painful!

I finally got an excuse this week to try out ScaLAPACK, a library for distributed-memory parallel dense linear algebra computations.  I've mainly been working on iterative methods for sparse problems, so it's not common for me to want to factor some giant dense matrix.  However, part of those iterative methods involves an orthogonalization method for small groups of dense vectors.  Comparing against ScaLAPACK's QR factorization PDGEQRF should make a good sanity check for accuracy and performance.  (PDGEQRF should be slower, since it uses an algorithm that requires more communication, and also since it likely doesn't use multicore parallelism as effectively for the panel factorization.)

It took me two and a half days to get a prototype accuracy test and benchmark up and running, which is a lot less than I had thought.  The Intel MKL Link Line Advisor was a big help in linking the right libraries.  I've tried building ScaLAPACK from source before -- while that was some four years ago, it was a real pain and I was happy not to have to try that.  The last few hours of those two-and-a-half days were spent trying to figure out why the matrix distribution wasn't working.  The routine in question, PDGEQRF, either refused to accept the input data (it helpfully identified the first argument of the routine which was off in some way, but didn't say why it was wrong), or only factored part of the matrix.  Staring at the ScaLAPACK Users' Guide -- in particular, at the diagrams illustrating global matrix layout -- helped me decode the mysterious parameters describing how the matrix is distributed among processors.

Using ScaLAPACK turned out not to be so difficult as I had thought.  The Users' Guide is well-written, with many helpful diagrams.  The routines have great documentation (and helpfully repeat documentation for certain library-standard interfaces, like the matrix descriptor array).  Having a pre-built library with the linker flags ready for copy-and-paste into my build configuration system was a huge help.  I imagine much of the pain people have reported came from trying to build everything themselves.  Note also that Matlab now has support for distributed-memory parallel computation, and distributed-memory data layouts.  It's based on ScaLAPACK, which means the MathWorks' developers have gone through a lot of this pain for you ;-)

A final interesting point:  the ScaLAPACK Users' Guide was written when High Performance Fortran (HPF) was under serious consideration as the way to write distributed-memory parallel code.  ScaLAPACK doesn't use HPF, but it does use HPF as an alternate notation to describe different ways to distribute matrices over processors.  HPF proved unpopular because it was too hard to build an efficient compiler for it, but parts of HPF made their way into modern versions of the Fortran standard. 
 

The user experience of lazy automatic tuning

I've been doing some performance experiments with an MPI-parallel linear algebra kernel.  I've implemented the kernel in two different ways:  a "reduce"  (not quite MPI_Reduce(), though it has a similar communication pattern) and "broadcast," and a "butterfly."  The latter is more general (it can be used to compute things other than what I'm benchmarking).  However, it might be slow in certain cases (e.g., if the network can't handle too many messages at once, or if there are many MPI ranks per node and the network card serializes on processing all those messages going in and out of the node).  Each run involves testing two different implementation alternatives, and four different scalar data types (float, double, complex float, and complex double).

When I ran the benchmarks, I noticed that the very first run (butterfly implementation strategy, "float" data type) of each invocation of the benchmark is 10 - 100x slower than for other data types with the same implementation strategy. When I rebuild the executable using a different MPI library, this problem went away.  What was going on?  I noticed that for the "slow" MPI library, the minimum benchmark run time was reasonable, but the maximum run time was a lot more.  After reading through some messages in the e-mail archives of that MPI library, I found out that the library delays some setup costs until execution time.  That means the first few calls of MPI collectives might be slow, but overall they will be fast.  The other MPI library doesn't seem to be doing this.

I was reminded of some troubles that the MathWorks had with integrating FFTW into Matlab, a few years ago.  FFTW does automatic performance tuning as a function of problem size.  The tuning phase takes time, but once it's finished, running any problem of that size will likely be much faster.  However, users didn't realize this.  They ran an FFT once, noticed that it was much slower than before (because FFTW was busy doing its thing), and then complained.  One can imagine all sorts of fixes for this, but the point is that something happened to performance which was not adequately communicated to users.

Now, MPI users are likely more sophisticated programmers than most Matlab users.  They likely have heard of automatic run-time performance tuning, and maybe even use it themselves.  I'm fairly familiar with such things and have been programming in MPI for a while.  Yet, this phenomenon still puzzled me until I poked around through the e-mail list archives.  What I was missing was an obvious pointer to documentation describing performance phenomena such as this one.  The performance of this MPI library is not transparent.  This is not a bad thing in this case, because the opaqueness is hiding a performance tuning process that will make my code faster overall.  However, I want to know this right away, so I don't have to think about what I could be doing wrong. 

See, I'm not just calling MPI directly.  I've wrapped up some other wrapper over MPI.  I didn't know if my wrapper was broken, the wrapper underneath was broken, the MPI library was broken, the cluster job execution software was broken, or if nothing was wrong and this was just expected behavior.  That's why I want to know right away whether I should expect significantly nonuniform performance behavior from some library.