The matrix and the table are probably the most common data structures in computing: matrices hold the data of scientific simulations and machine learning models, and tables hold the records of databases. On a single compute node, their main operations, the matrix multiplication and the database join, are well studied, with fast algorithms and libraries (e.g. BLAS for matrices). However, when the data is spread across a network of compute nodes (i.e. a distributed memory environment), each node reads only its own memory, and any other data must arrive as a message, which is much slower.

The data is spread because it is large: a dense matrix of 100,000 × 100,000 numbers in double precision takes 80 GB, and the product of two such matrices takes \(10^{15}\) multiply-adds. On 100 nodes, the three matrices of the product take only 2.4 GB per node. However, each element of the product needs a whole row of one matrix and a whole column of the other, which are spread across all nodes. How the data is split then decides how much of it each node receives, and whether it arrives while the node computes.

This post describes how to split the data of both operations so that each node works mostly alone, for dense matrices, sparse matrices and database tables. The algorithms use the terminology of MPI: they run on \(p\) processes, called ranks (usually one per compute node), or servers in the database part.

Part 1: dense matrices

The product \(C = AB\) of two \(n \times n\) matrices is given by:

\[C_{ij} = \sum_{k=0}^{n-1} A_{ik} \, B_{kj},\]

where each element takes \(n\) multiply-adds, i.e. \(n^3\) operations for the whole matrix, while the matrix-vector product \(y = Ax\) takes \(n^2\) operations.

Warm-up: the matrix-vector product

Each rank holds \(n/p\) consecutive rows of \(A\), and the matching entries of \(x\) and \(y\). To compute \(y_i\), a rank needs row \(i\) of \(A\), which it holds, and all of \(x\), which it does not. Therefore, the algorithm runs in two steps, illustrated below on 4 ranks:

  1. an all-gather (MPI_Allgather) sends the entries of \(x\) of each rank to all other ranks;
  2. each rank multiplies its rows of \(A\) by \(x\).

Matrix-vector product on 4 ranks, colored by rank. Left: the initial layout. Right: rank 1 after the all-gather, ready to compute \(y_2\) and \(y_3\).

Each rank performs \(n^2/p\) operations and receives about \(n\) entries of \(x\). Note that \(y\) ends up split like \(x\), ready for the next product of an iterative solver, such as the conjugate gradient method.

The matrix-matrix product: shifting tiles

The same idea fails for \(C = AB\): each rank would need all of \(B\), and an all-gather would store the \(n^2\) elements of \(B\) on every rank, which defeats the purpose of splitting the data. Instead, each rank holds only one piece of each matrix. At every step, it multiplies the pieces that it holds, and then shifts some of them to a neighbor, cyclically (i.e. the first rank sends to the last). The shifted pieces belong to one matrix, or to two.

Shifting one matrix: the ring algorithm

The ranks form a ring, and rank \(i\) holds row block \(i\) of \(A\), \(B\) and \(C\). Its block of \(A\) is cut into \(p\) tiles \(A_{i,k}\), one for each row block \(B_k\) of \(B\), so that row block \(i\) of \(C\) is computed as:

\[C_i = \sum_{k=0}^{p-1} A_{i,k} \, B_k,\]

where each tile \(A_{i,k}\) has \(n/p \times n/p\) elements, and each block \(B_k\) has \(n/p \times n\) elements. \(A\) and \(C\) stay in place, and only \(B\) moves. At each step, every rank multiplies the tile of \(A\) that matches its current block of \(B\), adds the result to \(C_i\), and passes the block of \(B\) to the previous rank. After \(p\) steps, \(C_i\) is complete, as shown below on 3 ranks:

The ring algorithm on 3 ranks. Blocks of \(B\) shift up, and rank 0 sends to rank 2. At each step, a rank multiplies the tile of \(A\) that has the color of its block of \(B\).

The output \(C\) is split by rows, like \(A\). In the mirror case, \(B\) stays in place, the column blocks of \(A\) move, and \(C\) is split by columns. Either way, each rank stores \(3n^2/p\) elements, but it receives \(p-1\) blocks of \(n^2/p\) elements, i.e. \((p-1)\,n^2/p\), almost \(n^2\), however many ranks are added.

Shifting two matrices: Cannon’s algorithm

To send less data, the matrices are cut in both dimensions, which is the idea of Cannon’s algorithm, first described by Lynn Elliot Cannon in 1969. The ranks form a \(q \times q\) process grid (with \(q = \sqrt{p}\)), and each matrix is cut into \(q \times q\) tiles of \(n/q \times n/q\) elements. Rank \((i,j)\) starts with \(A_{i,j}\) and \(B_{i,j}\), and computes:

\[C_{i,j} = \sum_{k=0}^{q-1} A_{i,k} \, B_{k,j}.\]

As an example, rank \((1,2)\) of a \(3 \times 3\) grid computes \(C_{1,2} = A_{1,0} B_{0,2} + A_{1,1} B_{1,2} + A_{1,2} B_{2,2}\). It needs the tiles of row 1 of \(A\) and of column 2 of \(B\), but only one matching pair at a time. Cannon’s algorithm delivers these pairs in three steps:

  1. row \(i\) of \(A\) is shifted left by \(i\), and column \(j\) of \(B\) up by \(j\) (the alignment), so that rank \((i,j)\) holds the matching pair \(A_{i,k}\) and \(B_{k,j}\), with \(k = (i+j) \bmod q\);
  2. every rank adds \(A_{i,k} B_{k,j}\) to \(C_{i,j}\);
  3. \(A\) is shifted left by one and \(B\) up by one, so that rank \((i,j)\) holds the next matching pair, \(A_{i,k+1}\) and \(B_{k+1,j}\) (indices modulo \(q\)).

Steps 2 and 3 are repeated until each rank has done \(q\) multiplications, as illustrated below:

Cannon’s algorithm on a \(3 \times 3\) grid of ranks. Tiles of \(A\) (top) shift left, and tiles of \(B\) (bottom) shift up. The thick frame marks rank \((1,2)\), and colors show its three pairs: blue for \(k=0\), green for \(k=1\) and red for \(k=2\).

Looking at the picture above, rank \((1,2)\) starts with \(A_{1,2}\) and \(B_{1,2}\), which do not match. The alignment shifts row 1 of \(A\) left by one and column 2 of \(B\) up by two: \(A_{1,0}\) arrives from rank \((1,0)\) and \(B_{0,2}\) from rank \((0,2)\), the pair \(k = (1+2) \bmod 3 = 0\) (blue). The next two shifts bring the pairs \(k = 1\) (green) and \(k = 2\) (red), and the three products add up to \(C_{1,2}\).

In brief, \(A\) moves along the grid rows, \(B\) along the grid columns, and \(C\) stays, so that tile \(C_{i,j}\) ends on rank \((i,j)\), in the same layout as the inputs. Note that, on a 2D grid, two of the three matrices must move.

Each rank still stores \(3n^2/p\) elements, but receives at most \(2q\) tiles (one pair in the alignment, and one pair at each of the \(q-1\) shifts), or \(2n^2/\sqrt{p}\) elements. Unlike in the ring algorithm, this volume shrinks as ranks are added, and it is the smaller of the two from 9 ranks on. On 4 ranks, a rank receives up to \(n^2\) elements with Cannon’s algorithm and \(0.75\,n^2\) with the ring; on 9 ranks, \(0.67\,n^2\) and \(0.89\,n^2\); and on 64 ranks, \(0.25\,n^2\) and \(0.98\,n^2\).

In MPI, the process grid is a communicator with a Cartesian topology, created by MPI_Cart_create with periodic dimensions, so that the shifts wrap around. For each shift, MPI_Cart_shift returns the ranks to send to and to receive from along a row or a column (at a distance of \(i\) or \(j\) for the alignment, and of one afterwards), and MPI_Sendrecv_replace sends the current tile and receives the next one into the same buffer.

Cannon’s algorithm needs a square grid of ranks. To find more about the topic, SUMMA (van de Geijn and Watts, 1995) replaces the shifts with broadcasts along the grid rows and columns, and works on rectangular grids. The 2.5D algorithm (Solomonik and Demmel, 2011) stores extra copies of the matrices to send even less data.

Overlapping computation and communication

The tiles of the next step do not depend on the current product, so a rank can receive them while it multiplies:

  1. the rank starts a non-blocking send of the current tiles, and a non-blocking receive of the next ones into a second buffer (MPI_Isend, MPI_Irecv);
  2. it multiplies the current tiles;
  3. it waits for both messages (MPI_Waitall), and swaps the buffers.

The following picture compares the timelines with and without this overlap:

Timeline of one rank in Cannon’s algorithm on a \(3 \times 3\) grid. Top: shifts and multiplications take turns. Bottom: each shift runs during a multiplication.

This technique, called double buffering, costs one extra pair of tiles, and works for the ring algorithm too. A shift is hidden whenever the multiplication takes longer. In practice, a multiplication costs \((n/q)^3\) operations and a shift moves \(2(n/q)^2\) elements, i.e. \(n/(2q)\) operations for each element received, so large tiles hide the communication easily. As an example, with \(n = 10{,}000\) on a \(10 \times 10\) grid, the tiles have 1,000 × 1,000 elements, and a rank performs 500 multiply-adds for every element that it receives.

Part 2: sparse matrices

Many matrices are mostly zeros, e.g. in finite element simulations, where each point of a mesh interacts only with a few neighbors. As an example, on a mesh of 1,000 × 1,000 points where each point interacts with its 4 nearest neighbors, the matrix has \(10^{12}\) entries (8 TB in double precision), of which only about 5 million (0.0005%) are nonzero. Such matrices are called sparse, and a sparse format stores only their nonzeros.

The compressed sparse row format

The compressed sparse row format (CSR, also called the Yale format), in use since the mid-1960s, stores a matrix in three arrays:

  • values: the nonzeros, row by row;
  • col_idx: the column of each nonzero;
  • row_ptr: where each row starts, i.e. row \(i\) is at positions row_ptr[i] to row_ptr[i+1]-1 of the other two arrays.

An example is displayed below:

A sparse matrix (left) and its CSR arrays (right), colored by row. Row 3 is empty, so row_ptr[3] and row_ptr[4] are both 6.

For \(n\) rows and \(nnz\) nonzeros, CSR stores \(2 \, nnz + n + 1\) numbers: 22 in the example above, instead of 30 for the dense layout, and about 11 million for the matrix of the mesh, instead of \(10^{12}\). Reading a row is easy, but reading a column requires a search through all rows, which makes the distributed transpose hard. The product \(y = Ax\) visits each row and each nonzero once, in \(O(n + nnz)\) operations:

for (int i = 0; i < n; i++) {
  y[i] = 0;
  for (int k = row_ptr[i]; k < row_ptr[i+1]; k++)
    y[i] += values[k] * x[col_idx[k]];
}

The sparse matrix-vector product

As before, \(A\) is split by rows for the sparse matrix-vector product (SpMV). An all-gather of \(x\) would work, but a rank needs \(x_j\) only for the columns \(j\) where its rows have a nonzero, which is usually a small set. So only those entries are sent:

  1. once, in a setup step, each rank tells the owners of \(x\) which of their entries it needs;
  2. before every product, an exchange (MPI_Alltoallv) sends these entries to the ranks that need them. The received copies are called ghost entries, or the halo of the rank;
  3. each rank runs the CSR loop over its own and ghost entries.

The following picture illustrates the exchange on 3 ranks:

Sparse matrix-vector product on 3 ranks. Left: nonzeros colored by rank. A dot marks a nonzero whose entry of \(x\) is on another rank. Right: the entries of \(x\) on each rank after the exchange. Dashed boxes are ghost entries.

In the picture above, rank 0 needs \(x_2\) and \(x_5\) from the other ranks, and ranks 1 and 2 need one entry each: 4 entries cross the network, instead of the 12 of an all-gather, which sends each rank the 4 entries of the other two ranks. The communication now grows with the number of ghost entries, and not with \(n\).

For a matrix from a mesh, the ghost entries are only the entries at the border of the region of each rank. As an example, the mesh of 1,000 × 1,000 points on 100 ranks gives each rank at most 2,000 ghost entries for strips of 10 × 1,000 points, and at most 400 for squares of 100 × 100 points, while an all-gather would send each rank 990,000 entries. In practice, a graph partitioning tool assigns the rows to the ranks, so that the regions have equal sizes and short borders.

The sparse matrix-matrix product

In the sparse matrix-matrix product (SpGEMM, for generalized sparse matrix-matrix multiplication), both \(A\) and \(B\) are sparse, and row \(i\) of \(C = AB\) is given by a weighted sum of rows of \(B\):

\[C_{i,:} = \sum_{k \, : \, A_{ik} \neq 0} A_{ik} \, B_{k,:},\]

where \(B_{k,:}\) is row \(k\) of \(B\). So a rank needs row \(k\) of \(B\) whenever its rows of \(A\) have a nonzero in column \(k\), i.e. the same indices as the ghost entries, but whole rows instead of single entries. With the matrix of the picture above as \(A\), rank 0 needs rows 2 and 5 of \(B\) from the other ranks, rank 1 needs row 0, and rank 2 row 1. This row-by-row formulation (Gustavson, 1978) suits CSR, which stores both matrices by rows, and leads to the following algorithm:

  1. each rank requests the remote rows of \(B\) it needs;
  2. the owners send them in CSR, with two MPI_Alltoallv calls: first the row sizes, so that the receiver can allocate memory, then the indices and values;
  3. each rank builds each of its rows of \(C\) in a sparse accumulator, a dense array of \(n\) values with the list of the positions that hold a nonzero, so that the compression of the row into CSR visits only its nonzeros.

The number of nonzeros of \(C\) is not known in advance, so step 3 often runs twice: a symbolic pass counts the nonzeros and allocates memory, and a numeric pass computes the values. The output \(C\) is split by rows, like \(A\).

To find more about the topic, Çatalyürek and Aykanat assign the rows of a sparse matrix to the ranks with a hypergraph instead of a graph, as the hypergraph counts the ghost entries of the matrix-vector product exactly. Buluç and Gilbert note that the row split of the sparse matrix-matrix product does not scale to thousands of ranks, and split both matrices on a 2D grid of ranks instead, as SUMMA does for dense matrices.

Part 3: database joins

In 1970, Edgar F. Codd proposed the relational model, where the data is stored in tables, with one record per row. A key is a column that identifies a record, and links tables. In the example below, Person has one row per person, with key id, and Messages has one row per message, whose id is the key of its author (a foreign key).

Inner and outer joins

A join pairs the rows of two tables whose keys match. When the condition is the equality of the keys, the join is called an equi-join, as in this SQL query:

SELECT Person.id, name, age, msg, year
FROM Person INNER JOIN Messages ON Person.id = Messages.id;

There are four kinds of inner and outer joins:

  • inner join: only the matching pairs. In the example below, Kim wrote two messages, so Kim appears twice, while Mary, Joe and the message with key 9 have no match, and are dropped;
  • left outer join: also the unmatched rows of the left table (Person), with NULL in the other columns;
  • right outer join: also the unmatched rows of the right table (Messages);
  • full outer join: also the unmatched rows of both tables.

The following picture shows the two tables and three of these joins:

Top: the tables Person (blue) and Messages (green). Bottom: three joins on id. Dark grey cells are NULL.

On one node, there are three classic algorithms to join tables of \(N\) and \(M\) rows (their costs exclude the size of the output):

  • nested loop: compares all pairs of rows, in \(O(N \cdot M)\);
  • sort-merge: sorts both tables by key and scans them together, in \(O(N \log N + M \log M)\), or \(O(N + M)\) if already sorted;
  • hash join: hashes the keys of the smaller table and looks up the rows of the other, in \(O(N + M)\) on average.

In practice, the hash join is the usual choice for large tables, and the sort-merge join is preferred when both tables are already sorted by key. The nested loop suits small tables, or join conditions other than an equality of keys. For outer joins, each row also records whether it found a match, and the rows without one are returned with NULL. However, all these algorithms assume that rows with the same key are on the same node, and the following sections show how to ensure this for three layouts.

Semi-join: each table on a different server

Here, Person is on server 1, and Messages is on server 2.

Inner join. Sending a whole table works, but a semi-join (Bernstein and Chiu, 1981) sends less:

  1. server 1 sends its keys to server 2;
  2. server 2 returns only its rows with one of these keys;
  3. server 1 joins them with its table, locally.

Steps 1 and 2 follow the pattern of the setup step of the sparse matrix-vector product, with keys in place of indices: the server that needs data sends the list of what it needs, and only that data comes back.

Outer joins. In step 3, server 1 also returns its unmatched persons, with NULL (the left outer join). Server 2 knows which of its messages had no match, as their keys were not in the list of step 1, so it returns them with NULL in a step 4. The inner join and step 4 give the right outer join, and steps 3 and 4 the full outer join, as displayed below:

A semi-join between two servers: the keys (1), the matching rows (2), the left outer join on server 1 (3), and the unmatched messages of server 2 (4). The first three rows of (3) are the inner join.

For a long list of keys, a Bloom filter is a smaller substitute. It lets a few wrong rows (false positives) through in step 2, and server 1 treats them as unmatched messages.

Shuffle and broadcast joins: tables split by rows

Each server holds some rows of each table, which is called horizontal partitioning, or sharding. Matching rows may then sit on different servers: below, John is on server 2, but his message “Bye” is on server 1. A local join misses this pair, and a local left outer join is even wrong, as it returns John with NULL.

Inner join. Each row is placed by its key, with the same rule for both tables (e.g. a hash of the key). Below, even keys go to server 1, and odd keys to server 2. Matching rows then share a server, and every server joins locally. If the tables are not placed this way, a shuffle first sends each row to its server, with an all-to-all exchange (MPI_Alltoallv), as illustrated below:

Two tables split by rows across 2 servers, colored by target server: blue for server 1 (even keys), and orange for server 2 (odd keys). Left: before the shuffle. Right: after it.

Outer joins. After a shuffle, outer joins are local too, as a row without a match on its server has no match anywhere. Above, server 2 returns Mary and Joe for the left outer join, and “Hey” for the right outer join. However, a shuffle may cause a load imbalance when a few keys are frequent (data skew), as all rows with the same key go to the same server.

Broadcast join. If Person is small, it can be copied to every server instead, while Messages stays in place, so that a frequent key does not overload any server. Each message then sits with all persons, so the inner and right outer joins are local.

The left outer join is not: John matches on server 1 but not on server 2, so no server can tell alone that a person has no match. Therefore, each server flags its local matches, and a logical-or all-reduce (MPI_Allreduce with MPI_LOR) combines the flags. Server 1 then returns the persons with a final flag of 0, with NULL, as presented below:

A broadcast join, for a left outer join. After the all-reduce of the match flags, only Mary and Joe have no match.

Join index: tables split by columns

Each server holds some columns, for all rows. This is called vertical partitioning, as in column-oriented databases, and the rows are aligned by position, i.e. row \(r\) on every server belongs to the same record. Below, server 1 holds the keys, server 2 the other columns of Person, and server 3 the other columns of Messages.

Inner join. Only the keys decide which rows match, so the keys are matched first, and the other columns move last:

  1. server 1 matches the keys into a join index (Valduriez, 1987): a list of pairs of row numbers, one from each table;
  2. server 1 sends the Person row numbers to server 2, and the Messages row numbers to server 3;
  3. servers 2 and 3 gather their rows in that order.

In practice, only row numbers cross the network, and the result stays split by columns. Building the output rows only at the end, and only for the rows that match, is called late materialization.

Outer joins. The join index also holds pairs with one side missing, such as Mary (row 1) and Joe (row 3) in the left outer join, and the servers write NULL where a row number is missing, as illustrated below:

A left outer join of tables split by columns across 3 servers. Column # is the implicit row number. Server 1 builds the join index (1) and sends row numbers (2), and servers 2 and 3 gather their rows (3). The first three rows are the inner join.

To find more about the topic, Blanas et al. compare the shuffle join (called repartition join there), the broadcast join and the semi-join on a cluster of 100 nodes running MapReduce, and say when to use each. The survey of Abadi et al. on column-oriented databases details late materialization (section 4.4) and the joins of column stores (section 4.5).

Summary

To sum up, both problems follow one rule: the data that is used together should be kept together. When the layout does not fit the operation, a few communication steps move only the data that is needed, overlapped with the computation where possible.

In practice, the ring algorithm works for any number of ranks, while Cannon’s algorithm sends less data from 9 ranks on, but needs a square grid. A semi-join pays off when few rows match, a shuffle when both tables are large, a broadcast join when one table is small or a few keys are frequent, and a join index when the tables are stored by columns. The following tables compare the methods:

Matrix algorithm What moves Output split by Received per rank
matrix-vector all of \(x\) rows \(\approx n\)
ring row blocks of \(B\) rows \(\approx n^2\)
Cannon tiles of \(A\) and \(B\) tiles \(2n^2/\sqrt{p}\)
sparse matrix-vector ghost entries of \(x\) rows the ghost entries
sparse matrix-matrix rows of \(B\) at the ghost indices rows these rows
Join layout Inner join Outer joins Sent over the network
each table on a server semi-join each server adds its unmatched rows keys, matching rows
split by rows shuffle by key, or broadcast a small table local after a shuffle, all-reduce of flags after a broadcast moved rows, or the small table
split by columns join index pairs with a missing side row numbers