Communication-Optimal Parallel Algorithm for Strassen's Matrix Multiplication

Introduction

Matrix multiplication is one of the most fundamental algorithmic problems in numerical linear algebra, distributed computing, scientific computing, and high-performance computing. Parallelization of matrix multiplication has been extensively studied (e.g., ). It has been addressed using many theoretical approaches, algorithmic tools, and software engineering methods in order to optimize performance and obtain faster and more efficient parallel algorithms and implementations.

We obtain a new parallel algorithm based on Strassen’s fast matrix multiplication.Our actual implementation uses the Winograd variant ; see Appendix A for details. It is more efficient than any other parallel matrix multiplication algorithm of which we are aware, including those that are based on classical ($\Theta(n^{3})$) multiplication, and those that are based on Strassen’s and other Strassen-like matrix multiplications. We compare the efficiency of the new algorithm with previous algorithms, and provide both asymptotic analysis (Sections 3 and 4) and benchmarking data (Section 5).

To design efficient parallel algorithms, it is necessary not only to load balance the computation, but also to minimize the time spent communicating between processors. The inter-processor communication costs are in many cases significantly higher than the computational costs. Moreover, hardware trends predict that more problems will become communication-bound in the future . Even matrix multiplication becomes communication-bound when run on sufficiently many processors. Given the importance of communication costs, it is preferable to match the performance of an algorithm to a communication lower bound, obtaining a communication-optimal algorithm.

2 Communication costs of matrix multiplication

We consider a distributed-memory parallel machine model as described in Section 2.1. The communication costs are measured as a function of the number of processors $P$, the local memory size $M$ in words, and the matrix dimension $n$. Irony, Toledo, and Tiskin proved that in the distributed memory parallel model, the bandwidth cost of classical $n$-by-$n$ matrix multiplication is bounded by $\Omega\left(\frac{n^{3}}{PM^{1/2}}\right)$ words. Using their technique one can also deduce a memory-independent bandwidth cost bound of $\Omega\left(\frac{n^{2}}{P^{2/3}}\right)$ and generalize it to other classes of algorithms . For a shared-memory model similar bounds were shown in . Until recently, parallel classical matrix multiplication algorithms (e.g., “2D” , and “3D” ) have minimized communication only for specific $M$ values. The first algorithm that minimizes the communication costs for the entire range of $M$ has recently been obtained by Solomonik and Demmel . See Section 4.1 for more details.

None of these lower bounding techniques and parallelization approaches generalize to fast matrix multiplication, such as . A communication cost lower bound for fast matrix multiplication algorithms has only recently been obtained : Strassen’s algorithm run on a distributed-memory parallel machine has bandwidth cost $\Omega\left(\left(\frac{n}{M^{1/2}}\right)^{\omega_{0}}\cdot\frac{M}{P}\right)$ and latency cost $\Omega\left(\left(\frac{n}{M^{1/2}}\right)^{\omega_{0}}\cdot\frac{1}{P}\right)$, where $\omega_{0}=\log_{2}7$ (see Section 2.4). These bounds generalize to other, but not all, fast matrix multiplication algorithms, with $\omega_{0}$ being the exponent of the computational complexity.

In the sequential case,See for a discussion of the sequential memory model. the lower bounds are attained by the natural recursive implementation which is thus optimal. However, a parallel communication-optimal Strassen-based algorithm was not previously known. Previous parallel algorithms that use Strassen (e.g., ), decrease the computational costs at the expense of higher communication costs. The factors by which these algorithms exceed the lower bounds are typically small powers of $P$ and $M$, as discussed in Section 4. However both $P$ and $M$ can be large (e.g. on a modern supercomputer, one may have $P\sim 10^{5}$ and $M\sim 10^{9}$).

3 Parallelizing Strassen’s matrix multiplication in a communication efficient way

The main impetus for this work was the observation of the asymptotic gap between the communication costs of existing parallel Strassen-based algorithms and the communication lower bounds. Because of the attainability of the lower bounds in the sequential case, we hypothesized that the gap could be closed by finding a new algorithm rather than by tightening the lower bounds.

We made three observations from the lower bound results of that lead to the new algorithm. First, the lower bounds for Strassen are lower than those for classical matrix multiplication. This implies that in order to obtain an optimal Strassen-based algorithm, the communication pattern for an optimal algorithm cannot be that of a classical algorithm but must reflect the properties of Strassen’s algorithm. Second, the factor $M^{\omega_{0}/2-1}$ that appears in the denominator of the communication cost lower bound implies that an optimal algorithm must use as much local memory as possible. That is, there is a tradeoff between memory usage and communication (the same is true in the classical case). Third, the proof of the lower bounds shows that in order to minimize communication costs relative to computation, it is necessary to perform each sub-matrix multiplication of size $\Theta(\sqrt{M})\times\Theta(\sqrt{M})$ on a single processor.

With these observations and assisted by techniques from previous approaches to parallelizing Strassen, we developed a new parallel algorithm which achieves perfect load balance, minimizes communication costs, and in particular performs asymptotically less computation and communication than is possible using classical matrix multiplication.

4 Our contributions and paper organization

Our main contribution is a new algorithm we call Communication-Avoiding Parallel Strassen, or CAPS.

CAPS asymptotically minimizes computational and bandwidth costs over all parallel Strassen-based algorithms. It also minimizes latency cost up to a logarithmic factor in the number of processors.

CAPS performs asymptotically better than any previous previous classical or Strassen-based parallel algorithm. It also runs faster in practice. The algorithm and its computational and communication cost analyses are presented in Section 3. There we show it matches the communication lower bounds.

We provide a review and analysis of previous algorithms in Section 4. We also consider two natural combinations of previously known algorithms (Sections 4.4 and 4.5). One of these new algorithms that we call “2.5D-Strassen” performs better than all previous algorithms, but is still not optimal, and performs worse than CAPS.

We discuss our implementations of the new algorithms and compare their performance with previous ones in Section 5 to show that our new CAPS algorithm outperforms previous algorithms not just asymptotically, but also in practice. Benchmarking our implementation on a Cray XT4, we obtain speedups over classical and Strassen-based algorithms ranging from $24\%$ to $184\%$ for a fixed matrix dimension $n=94080$, where the number of nodes ranges from $49$ to $7203$.

In Section 6 we show that our parallelization method applies to other fast matrix multiplication algorithms. It also applies to classical recursive matrix multiplication, thus obtaining a new optimal classical algorithm that matches the 2.5D algorithm of Solomonik and Demmel . In Section 6, we also discuss numerical stability, hardware scaling, and future work.

Preliminaries

We model communication of distributed-memory parallel architectures as follows. We assume the machine has $P$ processors, each with local memory of size $M$ words, which are connected via a network. Processors communicate via messages, and we assume that a message of $w$ words can be communicated in time $\alpha+\beta w$. The bandwidth cost of the algorithm is given by the word count and denoted by $BW(\cdot)$, and the latency cost is given by the message count and denoted by $L(\cdot)$. Similarly the computational cost is given by the number of floating point operations and denoted by $F(\cdot)$. We call the time per floating point operation $\gamma$.

We count the number of words, messages and floating point operations along the critical path as defined in . That is, two messages that are communicated between separate pairs of processors simultaneously are counted only once, as are two floating point operations performed in parallel on different processors. This metric is closely related to the total running time of the algorithm, which we model as

We assume that (1) the architecture is homogeneous (that is, $\gamma$ is the same on all processors and $\alpha$ and $\beta$ are the same between each pair of processors), (2) processors can send/receive only one message to/from one processor at a time and they cannot overlap computation with communication (this latter assumption can be dropped, affecting the running time by a factor of at most two), and (3) there is no communication resource contention among processors. That is, we assume that there is a link in the network between each pair of processors. Thus lower bounds derived in this model are valid for any network, but attainability of the lower bounds depends on the details of the network.

2 Strassen’s algorithm

Strassen showed that $2\times 2$ matrix multiplication can be performed using $7$ multiplications and $18$ additions, instead of the classical algorithm that does $8$ multiplications and $4$ additions . By recursive application this yields an algorithm with multiplies two $n\times n$ matrices $O(n^{\omega_{0}})$ flops, where $\omega_{0}=\log_{2}7\approx 2.81$. Winograd improved the algorithm to use $7$ multiplications and $15$ additions in the base case, thus decreasing the hidden constant in the $O$ notation . We review the Strassen-Winograd algorithm in Appendix A.

3 Previous work on parallel Strassen

In this section we breifly describe previous efforts to parallelize Strassen. More details, including communication analyses, are in Section 4. A summary appears in Table 1.

Luo and Drake explored Strassen-based parallel algorithms that use the communication patterns known for classical matrix multiplication. They considered using a classical 2D parallel algorithm and using Strassen locally, which corresponds to what we call the “2D-Strassen” approach (see Section 4.2). They also consider using Strassen at the highest level and performing a classical parallel algorithm for each subproblem generated, which corresponds to what we call the “Strassen-2D” approach. The size of the subproblems depends on the number of Strassen steps taken (see Section 4.3). Luo and Drake also analyzed the communication costs for these two approaches.

Soon after, Grayson, Shah, and van de Geijn improved on the Strassen-2D approach of by using a better classical parallel matrix multiplication algorithm and running on a more communication-efficient machine. They obtained better performance results compared to a purely classical algorithm for up to three levels of Strassen’s recursion.

Other parallel approaches have used more complex parallel schemes and communication patterns. However, they restrict attention to only one or two steps of Strassen and obtain modest performance improvements over classical algorithms.

4 Strassen lower bounds

For Strassen-based algorithms, the bandwidth cost lower bound has been proved using expansion arguments on the computation graph, and the latency cost lower bound is an immediate corollary.

(Memory-dependent lower bound) Consider a Strassen-based algorithm running on $P$ processors each with local memory size $M$. Let $BW(n,P,M)$ be the bandwith cost and $L(n,P,M)$ be the latency cost of the algorithm. Assume that no intermediate values are computed twice. Then

A memory-independent lower bound has recently been proved using the same expansion approach:

(Memory-independent lower bound) Consider a Strassen-based algorithm running on $P$ processors. Let $BW(n,P)$ be the bandwith cost and $L(n,P)$ be the latency cost of the algorithm. Assume that no intermediate values are computed twice. Assume only one copy of the input data is stored at the start of the algorithm and the computation is load-balanced in an asymptotic sense. Then

and the latency cost is $L(n,P)=\Omega(1)$.

Note that when $M=O(n^{2}/P^{2/\omega_{0}})$, the memory-dependent lower bound is dominant, and when $M=\Omega(n^{2}/P^{2/\omega_{0}})$, the memory-independent lower bound is dominant.

Communication-Avoiding Parallel Strassen

In this section we present the CAPS algorithm, and prove it is communication-optimal. See Algorithm 1 for a concise presentation and Algorithm 2 for a more detailed description.

CAPS has computational cost $\Theta\left(\frac{n^{\omega_{0}}}{P}\right)$, bandwidth cost $\Theta\left(\max\left\{\frac{n^{\omega_{0}}}{PM^{\omega_{0}/2-1}},\frac{n^{2}}{P^{2/\omega_{0}}}\right\}\right)$, and latency cost $\Theta\left(\max\left\{\frac{n^{\omega_{0}}}{PM^{\omega_{0}/2}}\log P,\log P\right\}\right)$.

By Theorems 2.1 and 2.2, we see that CAPS has optimal computational and bandwidth costs, and that its latency cost is at most $\log P$ away from optimal. Thus Theorem 1.1 follows. We prove Theorem 3.1 in Section 3.5.

Consider the recursion tree of Strassen’s sequential algorithm. CAPS traverses it in parallel as follows. At each level of the tree, the algorithm proceeds in one of two ways. A “breadth-first-step” (BFS) divides the 7 subproblems among the processors, so that $\frac{1}{7}$ of the processors work on each subproblem independently and in parallel. A “depth-first-step” (DFS) uses all the processors on each subproblem, solving each one in sequence. See Figure 1.

In short, a BFS step requires more memory but reduces communication costs while a DFS step requires little extra memory but is less communication-efficient. In order to minimize communication costs, the algorithm must choose an ordering of BFS and DFS steps that uses as much memory as possible.

Let $k=\log_{7}P$ and $s\geq k$ be the number of distributed Strassen steps the algorithm will take. In this section, we assume that $n$ is a multiple of $2^{s}7^{\lceil k/2\rceil}$. If $k$ is even, the restriction simplifies to $n$ being a multiple of $2^{s}\sqrt{P}$. Since $P$ is a power of 7, it is sometimes convenient to think of the processors as numbered in base 7. CAPS performs $s$ steps of Strassen’s algorithm and finishes the calculation with local matrix multiplication. The algorithm can easily be generalized to other values of $n$ by padding or dynamic peeling.

It is possible to use a more complicated scheme that interleave BFS and DFS steps to reduce communication. We show that the LM scheme is optimal up to a constant factor, and hence no more than a constant factor improvement can be attained from interleaving.

2 Data layout

We require that the data layout of the matrices satisfies the following two properties:

At each of the $s$ Strassen recursion steps, the data layouts of the four sub-matrices of each of $A$, $B$, and $C$ must match so that the weighted additions of these sub-matrices can be performed locally. This technique follows and allows communication-free DFS steps.

Each of these submatrices must be equally distributed among the $P$ processors for load balancing.

There are many data layouts that satisfy these properties, perhaps the simplest being block-cyclic layout with a processor grid of size $7^{\lfloor k/2\rfloor}\times 7^{\lceil k/2\rceil}$ and block size $\frac{n}{2^{s}7^{\lfloor k/2\rfloor}}\times\frac{n}{2^{s}7^{\lceil k/2\rceil}}$. (When $k=\log_{7}P$ is even these expressions simplify to a processor grid of size $\sqrt{P}\times\sqrt{P}$ and block size $\frac{n}{2^{s}\sqrt{P}}$.) See Figure 2.

Any layout that we use is specified by three parameters, $(n,P,s)$, and intermediate stages of the computation use the same layout with smaller values of the parameters. A BFS step reduces a multiplication problem with layout parameters $(n,P,s)$ to seven subproblems with layout parameters $(n/2,P/7,s-1)$. A DFS step reduces a multiplication problem with layout parameters $(n,P,s)$ to seven subproblems with layout parameters $(n/2,P,s-1)$.

Note that if the input data is initially load-balanced but distributed using a different layout, we can rearrange it to the above layout using a total of $O\left(\frac{n^{2}}{P}\right)$ words and $O(n^{2})$ messages. This has no asymptotic effect on the bandwidth cost but significantly increases the latency cost in the worst case.

3 Unlimited Memory scheme

In the UM scheme, we take $k=\log_{7}P$ BFS steps in a row. Since a BFS step reduces the number of processors involved in each subproblem by a factor of 7, after $k$ BFS steps each subproblem is assigned to a single processor, and so is computed locally with no further communication costs. We first describe a BFS step in more detail.

The matrices $A$ and $B$ are initially distributed as described in Section 3.2. In order to take a recursive step, the 14 matrices $S_{1},\dots S_{7},T_{1},\dots,T_{7}$ must be computed. Each processor allocates space for all 14 matrices and performs local additions and subtractions to compute its portion of the matrices. Recall that the submatrices are distributed identically, so this step requires no communication. If the layouts of $A$ and $B$ have parameters $(n,P,s)$, the $S_{i}$ and the $T_{i}$ now have layout parameters $(n/2,P,s-1)$.

The next step is to redistribute these 14 matrices so that the 7 pairs of matrices $(S_{i},T_{i})$ exist on disjoint sets of $P/7$ processors. This requires disjoint sets of $7$ processors performing an all-to-all communication step (each processor must send and receieve a message from each of the other 6). To see this, consider the numbering of the processors base-7. On the $m^{\textrm{th}}$ BFS step, the communication is between the seven processors whose numbers agree on all digits except the $m^{\textrm{th}}$ (counting from the right). After the $m^{\textrm{th}}$ BFS step, the set of processors working on a given subproblem share the same $m$-digit suffix. After the above communication is performed, the layout of $S_{i}$ and $T_{i}$ has parameters $(n/2,P/7,s-1)$, and the sets of processors that own the $T_{i}$ and $S_{i}$ are disjoint for different values of $i$. Note that since each all-to-all only involves seven processors no matter how large $P$ is, this algorithm does not have the scalability issues that typically come from an all-to-all communication pattern.

The extra memory required to take one BFS step is the space to store all $7$ triples $S_{j}$, $T_{j}$, $Q_{j}$. Since each of those matrices is $\frac{1}{4}$ the size of $A$, $B$, and $C$, the extra space required at a given step is $7/4$ the extra space required for the previous step. We assume that no extra memory is required for the local multiplications.If one does not overwrite the input, it is impossible to run Strassen in place; however using a few temporary matrices affects the analysis here by a constant factor only. Thus, the total local memory requirement for taking $k$ BFS steps is given by

3.2 Computational costs

The computation required at a given BFS step is that of the local additions and subtractions associated with computing the $S_{i}$ and $T_{i}$ and updating the output matrix $C$ with the $Q_{i}$. Since Strassen performs 18 additions and subtractions, the computational cost recurrence is

with base case $F_{\text{UM}}(n,1)=c_{s}n^{\omega_{0}}-6n^{2}$, where $c_{s}$ is the constant of Strassen’s algorithm. See Appendix A for more details. The solution to this recurrence is

3.3 Communication costs

Consider the communication costs associated with the UM scheme. Given that the redistribution within a BFS step is performed by an all-to-all communication step among sets of $7$ processors, each processor sends $6$ messages and receives $6$ messages to redistribute $S_{1},\dots,S_{7}$, and the same for $T_{1},\dots,T_{7}$. After the products $Q_{i}=S_{i}T_{i}$ are computed, each processor sends $6$ messages and receive $6$ messages to redistribute $Q_{1},\dots,Q_{7}$. The size of each message varies according to the recursion depth, and is the number of words a processor owns of any $S_{i}$, $T_{i}$, or $Q_{i}$, namely $\frac{n^{2}}{4P}$ words.

As each of the $Q_{i}$ is computed simultaneously on disjoint sets of $P/7$ processors, we obtain a cost recurrence for the entire UM scheme:

with base case $L_{\text{UM}}(n,1)=BW_{\text{UM}}(n,1)=0$. Thus

4 Limited Memory scheme

Consider taking a single DFS step. Rather than allocating space for and computing all 14 matrices $S_{1},T_{1},\dots,S_{7},T_{7}$ at once, the DFS step requires allocation of only one subproblem, and each of the $Q_{i}$ will be computed in sequence.

Consider the $i^{\text{th}}$ subproblem: as before, both $S_{i}$ and $T_{i}$ can be computed locally. After $Q_{i}$ is computed, it is used to update the corresponding quadrants of $C$ and then discarded so that its space in memory (as well as the space for $S_{i}$ and $T_{i}$) can be re-used for the next subproblem. In a DFS step, no redistribution occurs. After $S_{i}$ and $T_{i}$ are computed, all processors participate in the computation of $Q_{i}$.

We assume that some extra memory is available. To be precise, assume the matrices $A$, $B$, and $C$ require only $\frac{1}{3}$ of the available memory:

The extra memory requirement for a DFS step is the space to store one subproblem. Thus, the extra space required at this step is $1/4$ the space required to store $A$, $B$, and $C$. The local memory requirements for the LM scheme is given by

where the last line follows from (3) and (2). Thus, the limited memory scheme does not exceed the available memory.

4.2 Computational costs

As in the UM case, the computation required at a given DFS step is that of the local additions and subtractions associated with computing the $S_{i}$ and $T_{i}$ and updating the output matrix $C$ with the $Q_{i}$. However, since all processors participate in each subproblem and the subproblems are computed in sequence, the recurrence is given by

4.3 Communication costs

Since there are no communication costs associated with a DFS step, the recurrence is simply

Thus the total communication costs are given by

5 Communication optimality

(of Theorem 3.1). In the case that $M\geq\text{Mem}_{\text{UM}}(n,P)=\Omega\left(\frac{n^{2}}{P^{2/\omega_{0}}}\right)$ the UM scheme is possible. Then the communication costs are given by (1) which matches the lower bound of Theorem 2.2. Thus the UM scheme is communication-optimal (up to a logarithmic factor in the latency cost and assuming that the data is initially distributed as described in Section 3.2). For smaller values of $M$, the LM scheme must be used. Then the communication costs are given by (4) and match the lower bound of Theorem 2.1, so the LM scheme is also communication-optimal.

We note that for the LM scheme, since both the computational and communication costs are proportional to $\frac{1}{P}$, we can expect perfect strong scaling: given a fixed problem size, increasing the number of processors by some factor will decrease each cost by the same factor. However, this strong scaling property has a limited range. As $P$ increases, holding everything else constant, the global memory size $PM$ increases as well. The limit of perfect strong scaling is exactly when there is enough memory for the UM scheme. See for details.

Analysis of Other Algorithms

In the section we detail the asymptotic communication costs of other matrix multiplication algorithms, both classical and Strassen-based. These communication costs and the corresponding lower bounds are summarized in Table 1.

Many of the algorithms described in this section are hybrids of two different algorithms. We use the convention that the names of the hybrid algorithms are composed of the names of the two component algorithms, hyphenated. The first name describes the algorithm used at the top level, on the largest problems, and the second describes the algorithm used at the base level on smaller problems.

Classical algorithms must communicate asymptotically more than an optimal Strassen-based algorithm. To compare the lower bounds, it is necessary to consider three cases for the memory size: when the memory-dependent bounds dominate for both classical and Strassen, when the memory-dependent bound dominates for classical, but the memory-independent bound dominates for Strassen, and when the memory-independent bounds dominate for both classical and Strassen. This analysis is detailed in Appendix B. Briefly, the factor by which the classical bandwidth cost exceeds the Strassen bandwidth cost is $P^{a}$ where $a$ ranges from $\frac{2}{\omega_{0}}-\frac{2}{3}\approx 0.046$ to $\frac{3-\omega_{0}}{2}\approx 0.10$ depending on the relative problem size. The same sort of analysis is used throughout Section 4 to compare each algorithm with the Strassen-based lower bounds.

Various parallel classical matrix multiplication algorithms minimize communication relative to the classical lower bounds for certain amounts of local memory $M$. For example, Cannon’s algorithm minimizes communication for $M=O(n^{2}/P)$. Several more practical algorithms exist (such as SUMMA ) which use the same amount of local memory and have the same asymptotic communication costs. We call this class of algorithms “2D” because the communication patterns follow a two-dimensional processor grid.

Another class of algorithms, known as “3D" because the communication pattern maps to a three-dimensional processor grid, uses more local memory and reduces communication relative to 2D algorithms. This class of algorithms minimizes communication relative to the classical lower bounds for $M=\Omega(n^{2}/P^{2/3})$. As shown in , it is not possible to use more memory than $M=\Theta(n^{2}/P^{2/3})$ to reduce communication.

Recently, a more general algorithm has been developed which minimizes communication in all cases. Because it reduces to a 2D and 3D for the extreme values of $M$ but interpolates for the values between, it is known as the “2.5D” algorithm .

2 2D-Strassen

One idea to parallelize Strassen-based algorithms is to use a 2D classical algorithm for the inter-processor communication, and use the fast matrix multiplication algorithm locally . We call such an algorithm “2D-Strassen”. It is straightforward to implement, but cannot attain all the computational speedup from Strassen since it uses a classical algorithm for part of the computation. In particular, it does not use Strassen for the largest matrices, when Strassen provides the greatest reduction in computation. As a result, the computational cost exceeds $\Theta(n^{\omega_{0}}/P)$ by a factor of $P^{(3-\omega_{0})/2}\approx P^{0.10}$. The 2D-Strassen algorithm has the same communication cost as 2D algorithms, and hence does not match the communication costs of CAPS. In comparing the 2D-Strassen bandwidth cost, $\Theta(n^{2}/P^{1/2})$, to the CAPS bandwidth cost in Section 3, note that for the problem to fit in memory we always have $M=\Omega(n^{2}/P)$. The bandwidth cost exceeds that of CAPS by a factor of $P^{a}$, where $a$ ranges from $(3-\omega_{0})/2\approx.10$ to $2/\omega_{0}-1/2\approx.21$, depending on the relative problem size. Similarly, the latency cost, $\Theta(P^{1/2})$, exceeds that of CAPS by a factor of $P^{a}$ where $a$ ranges from $(3-\omega_{0})/2\approx.10$ to $1/2=.5$.

3 Strassen-2D

4 2.5D-Strassen

A natural idea is to replace a 2D classical algorithm in 2D-Strassen with the superior 2.5D classical algorithm to obtain an algorithm we call 2.5D-Strassen. This algorithm uses the 2.5D algorithm for the inter-processor communication, and then uses Strassen for the local computation. When $M=\Theta(n^{2}/P)$, 2.5D-Strassen is exactly the same as 2D-Strassen, but when there is extra memory it both decreases the communication cost and decreases the computational cost since the local matrix multiplications are performed (using Strassen) on larger matrices. To be precise, the computational cost exceeds the lower bound by a factor of $P^{a}$ where $a$ ranges from $1-\frac{\omega_{0}}{3}\approx 0.064$ to $\frac{3-\omega_{0}}{2}\approx 0.10$ depending on the relative problem size. The bandwidth cost exceeds the bandwidth cost of CAPS by a factor of $P^{a}$ where $a$ ranges from $\frac{2}{\omega_{0}}-\frac{2}{3}\approx 0.046$ to $\frac{3-\omega_{0}}{2}\approx 0.10$. In terms of latency, the cost of $\frac{n^{3}}{PM^{3/2}}+\log P$ exceeds the latency cost of CAPS by a factor ranging from $\log P$ to $P^{(3-\omega_{0})/2}\approx P^{0.10}$, depending on the relative problem size.

5 Strassen-2.5D

Performance Results

We have implemented CAPS using MPI on a Cray XT4, and compared it to various previous classical and Strassen-based algorithms. The benchmarking data is shown in Figure 3.

The nodes of the Cray XT4 have 8GB of memory and a quad-core AMD “Bupdapest” processor running at 2.3GHz. We treat the entire node as a single processor, and when we use the classical algorithm we call the optimized threaded BLAS in Cray’s LibSci to provide parallelism between the four cores in a node. The peak flop rate is 9.2 GFLOPS per core, or 36.8 GFLOPS per node. The machine consists of 9,572 nodes. All the data in Figure 3 is for multiplying two square matrices with $n=94080$.

2 Performance

Note that the vertical scale of Figure 3 is “effective GFLOPS”, which is a useful measure for comparing classical and fast matrix multiplication algorithms. It is calculated as

For classical algorithms, which perform $2n^{3}$ floating point operations, this gives the actual GFLOPS. For fast matrix multiplication algorithms it gives the performance relative to classical algorithms, but does not accurately represent the number of floating point operations performed.

Our algorithm outperforms all previous algorithms, and attains performance as high as 49.1 effective GFLOPS/node, which is $33\%$ above the theoretical maximum for all classical algorithms. Compared with the best classical implementation, our speedup ranges from $51\%$ for small values of $P$ up to $94\%$ when using most of the machine. Compared with the best previous parallel Strassen algorithms, our speedup ranges from $24\%$ up to $184\%$. Unlike previous Strassen algorithms, we are able to attain substantial speedups over the entire range of processors.

3 Strong scaling

Figure 3 is a strong scaling plot: the problem size is fixed and each algorithm is run with $P$ ranging from the minimum that provides enough memory up to the largest allowed value of $P$ smaller than the size of the machine. Perfect strong scaling corresponds to a horizontal line in the plot. As the communication analysis predicts, CAPS exhibits better strong scaling than any of the other algorithms (with the exception of ScaLAPACK, which obtains very good strong scaling by having poor performance for small values of $P$).

4 Details of the implementations

This implementation is the CAPS algorithm, with a few modifications from the presentation in Section 3. First, when computing locally it switches to classical matrix multiplication below some size $n_{0}$. Second, it is generalized to run on $P=c7^{k}$ processors for $c\in\{1,2,3,6\}$ rather than just $7^{k}$ processors. As a result, the base-case classical matrix multiplication is done on $c$ processors rather than 1. Finally, implementation uses the Winograd variant of Strassen; see Appendix A for more details. Every point in the plot is tuned to use the best interleaving pattern of BFS and DFS steps, and the best total number of Strassen steps. For points in the figure, the optimal total number of Strassen steps is always 5 or 6.

4.2 ScaLAPACK

We use ScaLAPACK as optimized by Cray in LibSci. This is an implementation of the SUMMA algorithm, and can run on an arbitrary number of processors. It should give the best performance if $P$ is a perfect square so the processors can be placed on a square 2D grid. All the runs shown in Figure 3 are with $P$ a perfect square.

4.3 2.5D classical

This is the code of . It places the $P$ processors in a grid of size $\sqrt{P/c}\times\sqrt{P/c}\times c$, and requires that $\sqrt{P/c}$ and $c$ are integers with $1\leq c\leq P^{1/3}$, and $c$ divides $\sqrt{P/c}$. Additionally, it gets the best performance if $c$ is as large as possible, given the constraint that $c$ copies of the input and output matrices fit in memory. In the case that $c=1$ this code is an optimized implementation of SUMMA. The values of $P$ and $c$ for the runs in Figure 3 are chosen to get the best performance.

4.4 Strassen-2D

Following the algorithm of , this implementation uses the DFS code from the implementation of CAPS at the top level, and then uses the optimized SUMMA code from the 2.5D implementation with $c=1$. Since the classical code requires that $P$ is a perfect square, this requirement applies here as well. The number of Strassen steps taken is tuned to give the best performance for each $P$ value, and the optimal number varies from 0 to 2.

4.5 2D-Strassen

Following the algorithm of , the 2D-Strassen implementation is analagous to the Strassen-2D implementation, but with the classical algorithm run before taking local Strassen steps. Similarly the same code is used for local Strassen steps here and in our implementation of CAPS. This code also requires that $P$ is a perfect square. The number of Strassen steps is tuned for each $P$ value, and the optimal number varies from 0 to 3.

4.6 2.5D-Strassen

This implementation uses the 2.5D implementation to reduce the problem to one processor, then takes several Strassen steps. The processor requirements are the same as for the 2.5D implementation. The number of Strassen steps is tuned for each number of processors, and the optimal number varies from 0 to 3. We also tested the Strassen-2.5D algorithm, but its performance was always lower than 2.5D-Strassen in our experiments.

Conclusions/Future Work

CAPS has the same stability properties as sequential versions of Strassen. For a complete discussion of the stability of fast matrix multiplication algorithms, see . We highlight a few main points here. The tightest error bounds for classical matrix multiplication are component-wise: $|C-\hat{C}|\leq n\epsilon|A|\cdot|B|,$ where $\hat{C}$ is the computed result and $\epsilon$ is the machine precision. Strassen and other fast algorithms do not satisfy component-wise bounds but do satisfy the slightly weaker norm-wise bounds: $\|C-\hat{C}\|\leq f(n)\epsilon\|A\|\|B\|,$ where $\|A\|=\max_{i,j}A_{ij}$ and $f$ is polynomial in $n$ . Accuracy can be improved with the use of diagonal scaling matrices: $D_{1}CD_{3}=D_{1}AD_{2}\cdot D_{2}^{-1}BD_{3}$. It is possible to choose $D_{1},D_{2},D_{3}$ so that the error bounds satisfy either $|C_{ij}-\hat{C}_{ij}|\leq f(n)\epsilon\|A(i,:)\|\|B(:,j)\|$ or $\|C-\hat{C}\|\leq f(n)\epsilon\||A|\cdot|B|\|$. By scaling, the error bounds on Strassen become comparable to those of many other dense linear algebra algorithms, such as LU and QR decomposition . Thus using Strassen for the matrix multiplications in a larger computation will often not harm the stability at all.

2 Hardware scaling

Although Strassen performs asymptotically less computation and communication than classical matrix multiplication, it is more demanding on the hardware. That is, if one wants to run matrix multiplication near the peak CPU speed, Strassen is more demanding of the memory size and communication bandwidth. This is because the ratio of computational cost to bandwidth cost is lower for Strassen than for classical. From the lower bounds in Section 2.4, the asymptotic ratio of computational cost to bandwidth cost is $M^{\omega_{0}/2-1}$ for Strassen-based algorithms, versus $M^{1/2}$ for classical algorithms. This means that it is harder to run Strassen near peak than it is to run classical matrix multiplication near peak. In terms of the machine parameters $\beta$ and $\gamma$ introduced in Section 2.1, the condition to be able to be computation-bound is $\gamma M^{1/2}\geq c\beta$ for classical matrix multiplication and $\gamma M^{\omega_{0}/2-1}\geq c^{\prime}\beta$ for Strassen. Here $c$ and $c^{\prime}$ are constants that depend on the constants in the communication and computational costs of classical and Strassen-based matrix multiplication.

The above inequalities may guide hardware design as long as classical and Strassen matrix multiplication are considered important computations. They apply both to the distributed case, where $M$ is the local memory size and $\beta$ is the inverse network bandwidth, and to the sequential/shared-memory case where $M$ is the cache size and $\beta$ is the inverse memory-cache bandwidth.

3 Optimizing on-node performance

Note that although our implementation performs above the classical peak performance, it performs well below the corresponding Strassen-Winograd peak, defined by the time it takes to perform $c_{s}n^{\omega_{0}}/P$ flops at the peak speed of each processor. To some extent, this is because Strassen is more demanding on the hardware, as noted above. However we have not yet analyzed whether the amount our performance is below Strassen peak can be entirely accounted for based on machine parameters. It is also possible that a high performance shared-memory Strassen implementation might provide substantial speedups for our implementation, as well as for 2D-Strassen and 2.5D-Strassen.

4 Testing on various architectures

We have implemented and benchmarked CAPS on only one architecture, a Cray XT4. It remains to check that it outperforms other matrix multiplication algorithms on a variety of architectures. On some architectures it may be more important to consider the topology of the network and redesign the algorithm to minimize contention, which we have not done.

5 Improvements to the algorithm

To be practically useful, it is important to generalize the number of processors on which CAPS can run. To attain the communication lower bounds, CAPS as presented in Section 3 must run on $P$ a power of seven processors. Of course, CAPS can then be run on any number of processors by simply ignoring no more than $\frac{6}{7}$ of them and incurring a constant factor overhead. Thus we can run on arbitrary $P$ and attain the communication and computation lower bounds up to a constant factor. However the computation time is still dominant in most cases, and it is preferable to attain the computation lower bound exactly. It is an open question whether any algorithm can run on arbitrary $P$, attain the computation lower bound exactly, and attain the communication lower bound up to a constant factor.

Moreover, although the communication costs of this algorithm match the lower bound up to a constant factor in bandwidth, and up to a $\log P$ factor in latency, it is an open question to determine the optimal constant in the lower bound and perhaps provide a new algorithm that matches it exactly. Note that in the analysis of CAPS in Section 3, the constants could be slightly improved.

6 Parallelizing other algorithms

We can apply our parallelization approach to recursive classical matrix multiplication to obtain a communication-optimal algorithm. This algorithm has the same asymptotic communication costs as the 2.5D algorithm . We observed comparable performance to the 2.5D algorithm on our experimental platform. As with CAPS, this algorithm has not been optimized for contention, whereas the 2.5D algorithm is very well optimized for contention on torus networks.

6.2 Other fast matrix multiplication algorithms

Our approach of executing a recursive algorithm in parallel by traversing the recursion tree in DFS (sequential) or BFS (parallel) manners is not limited to Strassen’s algorithm. All fast square matrix multiplication algorithms are built out of ways to multiply $n_{0}\times n_{0}$ matrices using $q<n_{0}^{3}$ multiplications. Like with Strassen and Strassen-Winograd, they compute $q$ linear combinations of entries of each of $A$ and $B$, multiply these pairwise, then compute the entries of $C$ as linear combinations of these.By , all fast matrix multiplication algorithms can be expressed in this bilinear form. CAPS can be easily generalized to any such multiplication, with the following modifications:

The number of processors $P$ is a power of $q$.

The data layout must be such that all $n_{0}^{2}$ blocks of $A$, $B$, and $C$ are distributed equally among the $P$ processors with the same layout.

The BFS and DFS determine whether the $q$ multiplications are performed in parallel or sequentially.

The communication costs are then exactly as above, but with $\omega_{0}=\log_{n_{0}}q$.

It is unclear whether any of the faster matrix multiplication algorithms are useful in practice. One reason is that the fastest algorithms are not explicit. Rather, there are non-constructive proofs that the algorithms exist. To implement these algorithms, they would have to be found, which appears to be a difficult problem. The generalization of CAPS described in this section does apply to all of them, so we have proved the existence of a communication-avoiding non-explicit parallel algorithm corresponding to every fast matrix multiplication algorithm. We conjecture that the algorithms are all communication-optimal, but that is not yet proved since the lower bound proofs in may not apply to all fast matrix multiplication algorithms. In cases where the lower bounds do apply, they match the performance of the generalization of CAPS, and so they are communication-optimal.

References

Appendix A Strassen-Winograd Algorithm

The Strassen-Winograd algorithm is usually preferred to Strassen’s algorithm in practice since it requires fewer additions. We use it for our implementation of CAPS. Divide the input matrices $A,B$ and output matrix $C$ into 4 submatrices:

Then form 7 linear combinations of the submatrices of each of $A$ and $B$, call these $T_{i}$ and $S_{i}$, respectively; multiply them pairwise; then form the submatrices of $C$ as linear combinations of these products:

This is one step of Strassen-Winograd. The algorithm is recursive since it can be used for each of the 7 smaller matrix multiplications. In practice, one often uses only a few steps of Strassen-Winograd, although to attain $O(n^{\omega_{0}})$ computational cost, it is necessary to recursively apply it all the way down to matrices of size $O(1)\times O(1)$. The precise computational cost of Strassen-Winograd is

Here $c_{s}$ is a constant depending on the cutoff point at which one switches to the classical algorithm. For a cutoff size of $n_{0}$, the constant is $c_{s}=(2n_{0}+4)/n_{0}^{\omega_{0}-2}$ which is minimized at $n_{0}=8$ yielding a computational cost of approximately $3.73n^{\omega_{0}}-5n^{2}$. If using Strassen-Winograd with cutoff $n_{0}$ this should be substituted into the computational cost expressions of Section 3.

Appendix B Communication-cost ratios

In this section we derive the ratio $R$ of classical to Strassen-based bandwidth cost lower bounds that appear in the beginning of Section 4. Note that both classical and Strassen-based lower bounds are attained by optimal algorithms. Similar derivations apply to the other ratios quoted in that section. Because the bandwidth cost lower bounds are different in the memory-dependent and the memory-independent cases, and the threshold between these is different for the classical and Strassen-based bounds, it is necessary to consider three cases.

$M=\Omega(n^{2}/P)$ and $M=O(n^{2}/P^{2/\omega_{0}})$. The first condition is necessary for there to be enough memory to hold the input and output matrices; the second condition puts both classical and Strassen-based algorithms in the memory-dependent case. Then the ratio of the bandwidth costs is:

Using the two bounds that define this case, we obtain $R=O(P^{(3-\omega_{0})/2})$ and $R=\Omega(P^{3/\omega_{0}-1})$.

$M=\Omega(n^{2}/P^{2/\omega_{0}})$ and $M=O(n^{2}/P^{2/3})$. This means that in the classical case the memory-dependent lower bound dominates, but in the Strassen-based case the memory-independent lower bound dominates. Then the ratio is:

Using the two bounds that define this case, we obtain $R=O(P^{3/\omega_{0}-1})$ and $R=\Omega(P^{2/\omega_{0}-2/3})$.

$M=O(P^{2/3})$. This means that both the classical and Strassen-based lower bounds are dominated by the memory-independent cases. Then the ratio is:

Overall, depending on the ratio of the problem size to the available memory, the factor by which the classical bandwidth costs exceed the Strassen-based bandwidth costs is between $\Theta(P^{2/\omega_{0}-2/3})$ and $\Theta(P^{(3-\omega_{0})/2})$.