Nobody has solved the fault resilience problem

Lessons from Supercomputing 2011:

1. Nobody has solved the programming models program.  This has never stopped anyone from writing unmaintainable but functional scientific codes, however.

2. Nobody has solved the fault resilience problem.  This actually terrifies everyone who has thought about it, because we have never had to deal with this problem before.  Furthermore, most of us are ignorant of the huge body of existing research and practice on fault-resilient computing in highly fault-prone environments with a huge penalty of failure (space-based computing, sometimes with human payloads).  Most of us balk at the price of fault resilience (triple modular redundancy) that this experience recommends.  Nevertheless, we're being asked by policymakers to do computations to help them make high-consequence decisions (asteroid-earth collisions? climate change? whether The Bomb might go off by mistake if unfortunate techs drop it on their big toe?). 

The #2 situation is comparable to that of floating-point arithmetic in the awful, awful days before industry and academia converged on a standard.  The lack of a sensible standard required immense floating-point expertise in order to have some confidence in the result of a computation.  These days, the lack of reliability of data and computations might require the same level of expertise in order to prove correctness of existing algorithms.  Around six years ago, Prof. William Kahan (UC Berkeley) suggested that perhaps only one PhD expert in floating-point computation graduates per year.  While his standards are exceptionally high, I have no doubt that the lack of training in numerical analysis is the bottleneck for making use of computers whose data and computations are likely to experience undetected errors.

Fastest Krylov method for symmetric indefinite problems?

One of my engineer correspondents asked what Krylov subspace method would be fastest for solving a symmetric indefinite problem: GMRES, a least-squares solver such as LSQR, or BiCG.  I've thought about this a lot for symmetric positive definite problems, but not so much for the indefinite case.  Here are my thoughts; i welcome correction or further suggestions.

BiCG reduces to a redundantly computed CG if the matrix is symmetric positive definite.  Thus, CG should always be faster than BiCG in that case.  I haven't thought about this question before for the case of a symmetric indefinite problem, but my guess is that MINRES would be faster than BiCG (or, say, QMR) in most cases.  MINRES promises monotonic convergence of the residual norm; BiCG doesn't.  Their storage requirements are comparable.

MINRES is a special case of GMRES for symmetric indefinite problems.  MINRES should be faster than GMRES, though MINRES is a bit less robust.

LSQR works a little bit like BiCG (it builds up two Krylov bases).  My guess is that MINRES would be faster than LSQR.

If the linear system is singular or very ill conditioned, one may have more success with LSQR than with GMRES, more success with GMRES than with MINRES, and more success with MINRES than SYMMLQ.  LSQR is able to solve rank-deficient problems in a least-squares sense.  GMRES and MINRES (depending on the implementation) may be able to solve rank-deficient problems in a least-squares sense.

On the other hand, "LSQR with preconditioning" means something rather different than users might think.

Alas, I want to write more, but time constrains me...

The ghosts of departed quantities

Yea, and Cauchy did gaze upon the paper, saying: "Behold, I have reduced the infinite to the finite; the ghosts of departed quantities I have constrained with sigils and brought again to life."

My friend and old colleague John Cassel responds:  "Which is why we refer to Cauchy today as having invented 'the calculus of the undead.'"

Newton:  "BLAAAAARGHARARRR FLUUUUXIONSSSSSS..." 

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!