Communication-Avoiding Algorithms and the Latency Floor

A Krylov solver on a large machine typically spends most of its wall-clock time waiting rather than computing. The inner products that dominate a conjugate gradient iteration are global reductions, and a global reduction has a latency that grows logarithmically with the process count while the arithmetic per rank shrinks linearly. At some process count the arithmetic becomes negligible and the iteration time is set entirely by how many times the ranks must agree on a scalar.

Communication-avoiding algorithms attack that ratio directly. They reorganize the linear algebra so that a block of s iterations needs one synchronization instead of s of them, at the cost of more arithmetic and more local memory. The idea generalizes beyond Krylov methods to any computation in which a sequence of dependent steps each requires a global exchange.

Why the inner product dominates

A conjugate gradient iteration performs one sparse matrix-vector product, two inner products, and three vector updates. The matrix-vector product costs a halo exchange plus local arithmetic, and its cost per rank is proportional to the number of nonzeros per rank. The inner products each require a global reduction, whose cost has a latency term proportional to the logarithm of the process count and a bandwidth term proportional to the vector length.

Those terms scale in opposite directions. As P grows at fixed problem size, the local arithmetic falls as 1/P and the halo volume falls more slowly, while the reduction latency grows as log P. At some point the reductions are the entire iteration. The allreduce cost model gives the floor: two reductions of a single scalar per iteration cost on the order of 4 log2(P) * alpha, which at a thousand ranks and a microsecond of latency is tens of microseconds per iteration before any arithmetic happens.

This is not a matter of a slow network. It is a consequence of the dependency structure: iteration k+1 cannot start until the scalars from iteration k are known, and knowing them requires every rank to contribute and every rank to receive the sum.

Blocking the recurrence

The standard fix is to run s iterations at once and delay the synchronizations. The mechanism is a change of basis. Instead of storing the Krylov basis vectors one at a time, the algorithm builds a block of s vectors using the matrix powers A^j r_0, computes the Gram matrix of that block with a single block inner product, and then solves a small s x s system to recover the coefficients that the individual iterations would have produced.

The consequence for communication is that the s inner products become one reduction over s scalars, and the s matrix-vector products can be fused into a kernel that reads the matrix once and multiplies by s vectors. The matrix read is the expensive part on a sparse matrix, so the fusion saves memory traffic as well as messages.

The price is explicit. The block of s basis vectors must be stored, which costs s vectors of memory per rank. The Gram matrix computation does s^2 inner products’ worth of arithmetic, which is s times more flops than the sequential version’s s inner products. The small system must be solved, and it becomes ill-conditioned as s grows, because the columns of the block approach the dominant eigenvector. Choosing s is therefore a numerical decision as much as a performance one.

#include <math.h>

/* One step of a blocked recurrence: given the current block of s
   vectors, orthonormalize against the previous block with a single
   reduction per column pair rather than one per vector dot product. */
void block_orthogonalize(double *W, const double *V, int s, long n)
{
    double h[16][16];            /* s <= 16 in this instantiation */
    if (s > 16) return;

    /* h = V^T W, computed as one block inner product. */
    for (int i = 0; i < s; i++)
        for (int j = 0; j < s; j++) {
            double acc = 0.0;
            for (long k = 0; k < n; k++)
                acc += V[i * n + k] * W[j * n + k];
            h[i][j] = acc;
        }

    /* W -= V h, a block rank-s update; one pass over memory. */
    for (int i = 0; i < s; i++)
        for (int j = 0; j < s; j++)
            for (long k = 0; k < n; k++)
                W[j * n + k] -= V[i * n + k] * h[i][j];
}

In a distributed setting the inner-product loop becomes a single collective over s^2 scalars, and the update is local. The reduction count per unit of progress drops by a factor of s.

A different exponent, not a smaller constant

The interesting claim about communication-avoiding algorithms is not that they reduce the constant in front of the latency term. It is that they change how the latency term scales with the problem size.

For a dense matrix multiplication distributed over P processors, the classical algorithm communicates a number of words proportional to n^2 / P^(1/2) and takes n^3 / P flops, so the communication-to-computation ratio grows as P^(1/2) / n. A blocked algorithm that uses fewer messages, at the cost of more flops, reduces the message count by a factor of P^(1/2) and keeps the word count at the same order. The exponent in the scalability analysis changes, which is qualitatively different from tuning a constant.

For sparse iterative methods the gain is in the number of sequential steps rather than in the total bytes. The blocked method still reads the matrix s times’ worth of data per s iterations, or once if the kernel fuses, but it performs one reduction per s iterations instead of s. The iteration count to convergence may also rise slightly, because the block version solves a modified problem.

The isoefficiency consequence is that a solver whose per-iteration reduction cost grows as log P can be restored to a per-iteration cost that is roughly independent of P, provided the block size grows with P and the problem size grows enough to keep the arithmetic substantial.

Numerical cost of the change

Every communication-avoiding variant trades stability for latency.

The block basis [r, Ar, A^2 r, ...] becomes progressively worse conditioned because repeated multiplication by A amplifies the dominant eigencomponent. Orthogonalizing that basis explicitly loses accuracy as the condition number of the block grows. The standard remedy is to use a different polynomial basis, such as a Newton or Chebyshev basis, whose columns are better conditioned, and to compute the small-system solution in that basis and convert back. The conversion itself must be done stably, which usually means a QR factorization of a small matrix rather than an explicit inverse.

Pipelined variants take a different approach: they reformulate the recurrence so that the inner products for iteration k+1 can be computed before the result of iteration k is known, by expressing them in terms of quantities already available. The reformulation is exact in exact arithmetic and changes the rounding, so the computed residual can drift from the true residual faster than in the classical recurrence. The practical consequence is that the algorithm must occasionally compute the true residual explicitly, which costs a matrix-vector product and a reduction, reintroducing exactly the synchronization the method was designed to remove.

The deeper numerical issue is that these methods change the answer. A classical conjugate gradient and a s-step communication-avoiding version applied to the same system will produce iterates that differ in the last digits and converge in a different number of iterations. The difference is usually small and occasionally not, particularly for ill-conditioned systems where the block basis is nearly dependent. Any claim about the method’s performance must therefore include the convergence behaviour, not just the time per iteration.

What else the idea applies to

The same restructuring applies wherever a dependence chain forces synchronization.

A sparse matrix-vector product on a graph with a high diameter can be restructured to compute several products at once by reading the matrix once and multiplying by a block of vectors, which raises arithmetic intensity without changing the sparsity. The technique is the basis of block and multi-vector Krylov methods and of the “sparse matrix, dense block” kernels used in spectral methods. The limit is memory: s vectors of length n must be resident, and the block must fit alongside the matrix.

A dense triangular solve has a strictly sequential dependency along the diagonal, and blocking it means computing a block of the solution at once and then correcting, which is the standard blocked forward substitution and converts a sequence of rank-one updates into a sequence of matrix-matrix products. The same transformation underlies the recursive formulations used in high-performance libraries.

A global scan, such as the prefix sums used in particle tracing and in sorting, has a dependency chain that appears to require log P steps, and a blelloch-style tree does exactly that. A communication-avoiding scan can compute a block of local scans and then combine them, reducing the number of global steps per element at the cost of a local pass over the data.

The common pattern is that the algorithm reads more data or does more arithmetic in order to shorten the critical path through the machine. On a network where a message costs a microsecond and a flop costs a fraction of a nanosecond, that trade is usually favourable, and its favourability increases with the process count. On a single node it is usually not, which is why the same code often has two paths and chooses between them by measuring the message latency. The interconnect determines whether the trade pays, and the answer differs between a torus-connected machine and a full-bisection tree.