This paper is best motivated by the following observation (on p. 2). Currently, the fastest machine in the TOP500 list is rated at 16.3 peta floating-point operations per second (PFLOPS), which is achieved on matrices of dimension 12,681,215. Twenty years ago, the speed was 170.6 giga FLOPS, and the corresponding matrix dimension was 42,000. This seems like tremendous progress, a factor of ten in speed every four years. However, the 1995 machine took less than five minutes to demonstrate this speed, whereas the 2015 machine took more than 23 hours. This is because the usual linear algebra algorithms have critical paths of length Ω (N) when dealing with N × N matrices, independent of the number of processors used.

Algorithms that have shorter (polynomial in log N) critical paths are known, but they perform asymptotically more work in all but the classical algorithms. This paper proposes algorithms for inverting N × N symmetric positive definite matrices where the critical path is polynomial in log N, but the work involved is at most a constant multiple of that of the classic algorithms.

The key idea is to use Strassen’s inversion [1] (not to be confused with Strassen’s multiplication). However, if we use this all the way, we get a critical path of length Ω(N). Hence, the authors stop at some point n and switch to Pan’s adaptation of Newton’s iteration [2], using O(log N) iterations (which suffices if the condition number of the matrix is polynomial in N). This does give the desired result: algorithm NeSt requires work (1+o(1))N³ and critical path length O(log^2.5+εN).

However, is it stable? This is a key question for nonstandard inversion algorithms. The authors in fact propose doing a single further Newton iteration after each Strassen step. This algorithms NeStS multiplies the work by five, but under appropriate conditions can be proved to be stable; in fact, it is so under much wider conditions in practice.

Modeling for a hypothetical (but pretty realistic) exascale machine shows that the basic NeST algorithm is 10 times faster than LU decomposition of N<1,000,000. It is a pity that the authors didn’t also model NeStS, as this seems theoretically recommended.

The authors also perform practical tests on a 4×8-core NUMA machine, measuring their algorithms against a recent (inversion has apparently improved recently) Intel MKL implementation. Their NeSt uses MKL for the addition and multiplication subroutines. On a single core, NeSt takes 1.27 times as long as MKL, better than theory would predict. On 32 cores, with N=2¹⁴, MKL seems to take 50 percent longer than NeSt.

There are then some rather surprising results on accuracy, performed on N=2¹³ matrices (always symmetric positive definite). Up to condition number 2⁶, the methods are similar. Then, pure Newton and, surprisingly, NeStS are worse, while NeSt tracks MKL up to condition number 2¹⁴, and then is worse (but still much better than NeStS or Newton).

Finally, the authors use the theoretical trick of reducing inversion of arbitrary matrices to symmetric positive definite ones. As expected, the accuracy behavior is dire (as this trick squares the condition number).