Graph Expansion and Communication Costs of Fast Matrix Multiplication

Introduction

The communication of an algorithm (e.g., transferring data between the CPU and memory devices, or between parallel processors, a.k.a. I/O-complexity) often costs significantly more time than its arithmetic. It is therefore of interest (1) to obtain lower bounds for the communication needed, and (2) to design and implement algorithms attaining these lower bounds. Communication also requires much more energy than arithmetic, and saving energy may be even more important than saving time.

Communication time varies by orders of magnitude, from $O\left(10^{-9}\right)$ second for an L1 cache reference, to $O\left(10^{-2}\right)$ second for disk access. The variation can be even more dramatic when communication occurs over networks or the internet. While Moore’s Law predicts an exponential increase of hardware density in general, the annual improvement rate of time-per-arithmetic-operation has, over the years, consistently exceeded that of time-per-word read/write [Graham et al. (2004), Fuller and Millett (2011)]. The fraction of running time spent on communication is thus expected to increase further.

We model communication costs of sequential and parallel architecture as follows. In the sequential case, with two levels of memory hierarchy (fast and slow), communication means reading data items (words) from slow memory (of unbounded size), to fast memory (of size $M$) and writing data from fast memory to slow memorySee [Ballard et al. (2010)] for definition of a model with memory hierarchy, and a reduction from the two-level model. All bounds in this paper thus apply to the model with memory hierarchy as well.. Words that are stored contiguously in slow memory can be read or written in a bundle which we will call a message. We assume that a message of $n$ words can be communicated between fast and slow memory in time $\alpha+\beta n$ where $\alpha$ is the latency (seconds per message) and $\beta$ is the inverse bandwidth (seconds per word). We define the bandwidth cost of an algorithm to be the total number of words communicated and the latency cost of an algorithm to be the total number of messages communicated. We assume that the input matrices initially reside in slow memory, and are too large to fit in the smaller fast memory. Our goal then is to minimize both bandwidth and latency costs.The sequential communication model used here is sometimes called the two-level I/O model or disk access machine (DAM) model (see [Aggarwal and Vitter (1988), Bender et al. (2007), Chowdhury and Ramachandran (2006)]). Our bandwidth cost model follows that of [Hong and Kung (1981)] and [Irony et al. (2004)] in that it assumes the block-transfer size is one word of data ($B=1$ in the common notation). However, our model allows message sizes to vary from one word up to the maximum number of words that can fit in fast memory.

In the parallel case, we assume $p$ processors, each with memory of size $M$ (or with larger memory size, as long as we never use more than $M$ in each processor). We are interested in the communication among the processors. As in the sequential case, we assume that a message of $n$ consecutively stored words can be communicated in time $\alpha+\beta n$. This cost includes the time required to “pack” non-contiguous words into a single message, if necessary. We assume that the input is initially evenly distributed among all processors, so $M\cdot p$ is at least as large as the input. Again, the bandwidth cost and latency cost are the word and message counts respectively. However, we count the number of words and messages communicated along the critical path as defined in [Yang and Miller (1988)] (i.e., two words that are communicated simultaneously are counted only once), as this metric is closely related to the total running time of the algorithm. As before, our goal is to minimize the number of words and messages communicated.

We assume that (1) the cost per flop is the same on each processor and the communication costs ($\alpha$ and $\beta$) are the same between each pair of processors (this assumption is for ease of presentation and can be dropped, using [Ballard et al. (2011)]; see Section 6.3), (2) all communication is “blocking”: a processor can send/receive a single message at a time, and cannot communicate and compute a flop simultaneously (the latter assumption can be dropped, affecting the running time by a factor of two at most), and (3) there is no communication resource contention among processors. For example, if processor 0 sends a message of size $n$ to processor 1 at time 0, and processor 2 sends a message of size $n$ to processor 3 also at time 0, the cost along the critical path is $\alpha+\beta n$. However, if both processor 0 and processor 1 try to send a message to processor 2 at the same time, the cost along the critical path will be the sum of the costs of each message.

2 The Computation Graph and Implementations of an Algorithm

The computation performed by an algorithm on a given input can be modeled (see Section 3) as a computation directed acyclic graph (CDAG) : We have a vertex for each input / intermediate / output argument, and edges according to direct dependencies (e.g., for the binary arithmetic operation $x:=y+z$ we have a directed edge from $v_{y}$ to $v_{x}$ and from $v_{z}$ to $v_{x}$, where the vertices $v_{x},v_{y},v_{z}$ stand for the arguments $x,y,z$, respectively).

An implementation of an algorithm determines, in the parallel model, which arithmetic operations are performed by which of the $p$ processors. This corresponds to partitioning the corresponding CDAG into $p$ parts. Edges crossing between the various parts correspond to arguments that are in the possession of one processor, but are needed by another processor, therefore relate to communication. In the sequential model, an implementation determines the order of the arithmetic operations, in a way that respects the partial ordering of the CDAG (see Section 3 relating this to communication cost).

Implementations of an algorithm may vary greatly in their communication costs. The I/O-complexity of an algorithm is the minimum bandwidth cost of the algorithm, over all possible implementations. The I/O-complexity of a problem is defined to be the minimum I/O-complexity of all algorithms for this problem. A lower bound of the I/O-complexity of an algorithm is therefore a results of the form: any implementation of algorithm $Alg$ requires at least $X$ communication. An upper bound is of the form: there is an implementation for algorithm $Alg$ that requires at most $X$ communication. We detail below some of the I/O-complexity lower and upper bounds of specific algorithms, or a class of algorithms. I/O-complexity lower bounds for a problem are claims of the form: any algorithm for a problem $P$ requires at least $X$ communication. These are much harder to find (but see for example [Demmel, Grigori, Hoemmen, and Langou, 2008]).

The lower bounds in this paper are for all implementations for a family of algorithms: “Strassen-like” fast matrix multiplication. Generally speaking, a “Strassen-like” algorithm utilizes an algorithm for multiplying two constant-size matrices in order to recursively multiply matrices of arbitrary size; see Section 5 for precise definition and technical assumptions.

3 Previous Work

Consider the classical $\Theta(n^{3})$ algorithm for matrix multiplicationBy which we mean any algorithm that computes using the $n^{3}$ multiplications, whether this is done recursively, iteratively, block-wise or any other way.. While naïve implementations are communication inefficient, communication-minimizing sequential and parallel variants of this algorithm were constructed, and proved optimal, by matching lower bounds [Cannon (1969), Hong and Kung (1981), Frigo et al. (1999), Irony et al. (2004)].

In [Ballard et al. (2010), Ballard et al. (2011c)] we generalize the results of [Hong and Kung (1981), Irony et al. (2004)] regarding matrix multiplication, to obtain new I/O-complexity lower bounds for a much wider variety of algorithms. Most of our bounds are shown to be tight. This includes all “classical” algorithms for $LU$ factorization, Cholesky factorization, $LDL^{T}$ factorization, and many for the $QR$ factorization, and eigenvalues and singular values algorithms. Thus we essentially cover all direct methods of linear algebra. The results hold for dense matrix algorithms (most of them have $O(n^{3})$ complexity), as well as sparse matrix algorithms (whose running time depends on the number of non-zero elements, and their locations). They apply to sequential and parallel algorithms, to compositions of linear algebra operations (like computing the powers of a matrix), and to certain graph-theoretic problemsSee [Michael et al. (2002)] for bounds on graph-related problems, and our [Ballard et al. (2011c)] for a detailed list of previously known and recently designed sequential and parallel algorithms that attain the above mentioned lower bounds..

The optimal algorithms for square matrix multiplication are well known. Optimal algorithms for dense LU, Cholesky, QR, eigenvalue problems and the SVD are more recent. These include [\citeNPGustavson97; \citeNPToledo97; \citeNPElmrothGustavson98; \citeNPFrigoLeisersonProkopRamachandran99; \citeNPAhmedPingali00; \citeNPFrensWise03; [Demmel, Grigori, Hoemmen, and Langou, 2008]; [Demmel, Grigori, and Xiang, 2008]; \citeNPBallardDemmelHoltzSchwartz09a; \citeNPDavidDemmelGrigoriPeyronnet10; \citeNPDemmelGrigoriXiang10; \citeNPBallardDemmelDumitriu10], and are not part of standard libraries like LAPACK [Anderson et al. (1992)] and ScaLAPACK [Blackford et al. (1997)]. See [Ballard et al. (2011c)] for more details.

In [Ballard et al. (2010), Ballard et al. (2011c)] we use the approach of [Irony et al. (2004)], based on the Loomis-Whitney geometric theorem [Loomis and Whitney (1949), Burago and Zalgaller (1988)], by embedding segments of the computation process into a three-dimensional cube. This approach, however, is not suitable when distributivity is used, as is the case in Strassen [Strassen (1969)] and other fast matrix multiplication algorithms (e.g., [Coppersmith and Winograd (1990), Cohn et al. (2005)]).

While the I/O-complexity of classic matrix multiplication and algorithms with similar structure is quite well understood, this is not the case for algorithms of more complex structure. The problem of minimizing communication in parallel classical matrix multiplication was addressed [Cannon (1969)] almost simultaneously with the publication of Strassen’s fast matrix multiplication [Strassen (1969)]. Moreover, an I/O-complexity lower bound for the classical matrix multiplication algorithm has been known for three decades [Hong and Kung (1981)]. Nevertheless, the I/O-complexity of Strassen’s fast matrix multiplication and similar algorithms has not been resolved.

In this paper we obtain first communication cost lower bounds for Strassen’s and other fast matrix multiplication algorithms, in the sequential and parallel models. These bounds are attainable both for sequential and for parallel algorithms and so optimal.

4 Communication Costs of Fast Matrix Multiplication

The I/O-complexity $IO(n)$ of Strassen’s algorithm (see Algorithm 1, Appendix A), applied to $n$-by-$n$ matrices on a machine with fast memory of size $M$, can be bounded above as follows (for actual uses of Strassen’s algorithm, see [Douglas et al. (1994), Huss-Lederman et al. (1996), Desprez and Suter (2004)]): Run the recursion until the matrices are sufficiently small. Then, read the two input sub-matrices into the fast memory, perform the matrix multiplication inside the fast memory, and write the result into the slow memoryHere we assume that the recursion tree is traversed in the usual depth-first order.. We thus have $IO(n)\leq 7\cdot IO\left(\frac{n}{2}\right)+O(n^{2})$ and $IO\left(\frac{\sqrt{M}}{3}\right)=O(M)$. Thus

4.2 Lower bound

In this paper, we obtain a tight lower bound:

(Main Theorem) The I/O-complexity $IO(n)$ of Strassen’s algorithm on a machine with fast memory of size $M$, assuming that no intermediate values are computed twiceWe assume no recomputation throughout the paper., is

It holds for any implementation and any known variant of Strassen’s algorithmThis lower bound for the sequential case seems to contradict the upper bound from FOCS’99 [Frigo et al. (1999), Blelloch et al. (2008)]), due to a miscalculation (see [Leiserson (2008)] for details). ,To obtain the lower bounds for latency costs we divide the bandwidth costs by the maximal message length, $M$. This holds for all the lower bounds here, both in the sequential and parallel models.. This includes Winograd’s $O(n^{\lg 7})$ variant that uses 15 additions instead of 18, which is the most used fast matrix multiplication algorithm in practice [Douglas et al. (1994), Huss-Lederman et al. (1996), Desprez and Suter (2004)].

For parallel algorithms, using a reduction from the sequential to the parallel model (see e.g., [Irony et al. (2004)] or our [Ballard et al. (2011c)]) this yields:

Let $IO(n)$ be the I/O-complexity of Strassen’s algorithm, run on a machine with $p$ processors, each with a local memory of size $M$. Assume that no intermediate values are computed twice. Then

We can extend these bounds to a wider class of all “Strassen-like” fast matrix multiplication algorithms. Note that this class does not include all fast matrix multiplication algorithms (see Section 5.1 for definition of “Strassen-like” algorithms, and in particular the technical assumption in Section 5.1.1). Let $Alg$ be any “Strassen-like” matrix multiplication algorithm that runs in time $O(n^{\omega_{0}})$ for some $2<\omega_{0}<3$. Then, using the same arguments that lead to (1), the I/O-complexity of $Alg$ can be shown to be $IO(n)=O\left(\left(\frac{n}{\sqrt{M}}\right)^{\omega_{0}}\cdot M\right)$. We obtain a matching lower bound:

The I/O-complexity $IO(n)$ of a recursive “Strassen-like” fast matrix multiplication algorithm with $O(n^{\omega_{0}})$ arithmetic operations, on a machine with fast memory of size $M$ is

Note that or the cubic recursive algorithm for matrix multiplication, $\omega_{0}=\lg 8=3$, and the above formula is $IO(n)=\Omega\left(\left(\frac{n}{\sqrt{M}}\right)^{3}\cdot M\right)=\Omega\left(\frac{n^{3}}{\sqrt{M}}\right)$ and identifies with the lower bounds of [Hong and Kung (1981)] and [Irony et al. (2004)]. While the lower bounds for $\omega_{0}=3$ and for $\omega_{0}<3$ have the same form, the proofs are completely different, and it is not clear whether our approach can be used to prove their lower bounds and vice versa.

Let $IO(n)$ be the I/O-complexity of a “Strassen-like” algorithm (with arithmetic performed as in Theorem 1.3), run on a machine with $p$ processors, each with a local memory of size $M$. Assume that no intermediate values are computed twice. Then

5 The Expansion Approach

The proof of the main theorem is based on estimating the edge expansion of the computation directed acyclic graph (CDAG) of an algorithm. The I/O-complexity is shown to be closely related to the edge expansion properties of this graph. As the graph has a recursive structure, the expansion can be analyzed directly (combinatorially, similarly to what is done in [Mihail (1989), Alon et al. (2008), Koucky et al. (2010)]) or by spectral analysis (in the spirit of what was done for the Zig-Zag expanders [Reingold et al. (2002)]). There is, however, a new technical challenge. The replacement product and the Zig-Zag product act similarly on all vertices. This is not what happens in our case: multiplication and addition vertices behave differently.

The expansion approach is similar to the one taken by Hong and Kung [Hong and Kung (1981)]. They use the red-blue pebble game to obtain tight lower bounds on the I/O-complexity of many algorithms, including classical $\Theta(n^{3})$ matrix multiplication, matrix-vector multiplication, and FFT. The proof is obtained by showing that the size of any subset of the vertices of the CDAG is bounded by a function of the size of its dominator set (recall that a dominator set $D$ for $S$ is a set of vertices such that every path from an input vertex to a vertex in $S$ contains some vertex in $D$).

On the one hand, their dominator set technique has the advantage of allowing recomputation of any intermediate value. We were not able to allow recomputation using our edge expansion approach. On the other hand, the dominator set requires large input or output. Such an assumption is not needed by the edge expansion approach, as the bounds are guaranteed by edge expansion of many (internal) parts of the CDAG. In that regard, one can view the approach of [Irony et al. (2004)] (also in [Demmel, Grigori, Hoemmen, and Langou, 2008; \citeNPBallardDemmelHoltzSchwartz10a; \citeNPBallardDemmelHoltzSchwartz11a]) as an edge expansion assertion on the CDAGs of the corresponding classical algorithms.

The study of expansion properties of a CDAG was also suggested as one of the main motivations of Lev and Valiant [Lev and Valiant (1983)] in their work on superconcentrators. They point out many papers proving that classes of algorithms computing DFT, matrix inversion and other problems all have to have CDAGs with good expansion properties, thus providing lower bounds on the number of the arithmetic operations required.

Other papers study connections between bounded space computation, and combinatorial expansion-related properties of the corresponding CDAG (see e.g., [Savage (1994), Bilardi and Preparata (1999), Bilardi et al. (2000)] and references therein).

6 Paper organization

Section 2 contains preliminaries on the notions of graph expansion. In Section 3 we state and prove the connection between I/O-complexity and the expansion properties of the computation graph. In Section 4 we analyze the expansion of the CDAG of Strassen’s algorithm. We discuss the generalization of the bounds to other algorithms in Section 5, and present conclusions and open problems in Section 6.

Preliminaries

The edge expansion $h(G)$ of a $d$-regular undirected graph $G=(V,E)$ is:

0.2 When G𝐺G is not regular

Note that CDAGs are typically not regular. If a graph $G=(V,E)$ is not regular but has a bounded maximal degree $d$, then we can add ($<d$) loops to vertices of degree $<d$, obtaining a regular graph $G^{\prime}$. We use the convention that a loop adds 1 to the degree of a vertex. Note that for any $S\subseteq V$, we have $|E_{G}(S,V\setminus S)|=|E_{G^{\prime}}(S,V\setminus S)|$, as none of the added loops contributes to the edge expansion of $G^{\prime}$.

0.3 Expansion of small sets

For many graphs, small sets expand better than larger sets. Let $h_{s}(G)$ denote the edge expansion for sets of size at most $s$ in $G$:

In many cases, $h_{s}(G)$ does not depend on $|V(G)|$, although it may decrease when $s$ increases. One way of bounding $h_{s}(G)$ is by decomposing $G$ into small subgraphs of large edge expansion.

Let $G=(V,E)$ be a $d$-regular graph that can be decomposed into edge-disjoint (but not necessarily vertex-disjoint) copies of a $d^{\prime}$-regular graph $G^{\prime}=(V^{\prime},E^{\prime})$. Then the edge expansion of $G$ for sets of size at most $|V^{\prime}|/2$ is $h(G^{\prime})\cdot\frac{d^{\prime}}{d}$, namely

For proving this claim, recall the definition of graph decomposition:

We say that the set of graphs $\{G^{\prime}_{i}=(V_{i},E_{i})\}_{i\in[l]}$ is an edge-disjoint decomposition of $G=(V,E)$ if $V=\bigcup_{i}V_{i}$ and $E=\biguplus_{i}E_{i}$.

(of Claim 1) Let $U\subseteq V$ be of size $U\leq|V^{\prime}|/2$. Let $\{G^{\prime}_{i}=(V_{i},E_{i})\}_{i\in[l]}$ be an edge-disjoint decomposition of $G$, where every $G_{i}$ is isomorphic to $G^{\prime}$. Then

Therefore $\frac{|E_{G}(U,V\setminus U)|}{d\cdot|U|}\geq h(G^{\prime})\cdot\frac{d^{\prime}}{d}~{}.$

I/O-Complexity and Edge Expansion

In this section we recall the notion of computation graph of an algorithm, then show how a partition argument connects the expansion properties of the computation graph and the I/O-complexity of the algorithm. A similar partition argument already appeared in [Irony et al. (2004)], and then in our [Ballard et al. (2011c)]. In both cases it is used to relate I/O-complexity to the Loomis-Whitney geometric bound [Loomis and Whitney (1949)], which can be viewed, in this context, as an expansion guarantee for the corresponding graphs.

For a given algorithm, we consider the computation (directed) graph $G=(V,E)$, where there is a vertex for each arithmetic operation (AO) performed, and for every input element. $G$ contains a directed edge $(u,v)$, if the output operand of the AO corresponding to $u$ (or the input element corresponding to $u$), is an input operand to the AO corresponding to $v$. The in-degree of any vertex of $G$ is, therefore, at most 2 (as the arithmetic operations are binary). The out-degree is, in general, unboundedAs the lower bounds are derived for the bounded out-degree case, we will show how to convert the corresponding CDAG to obtain constant out-degree, without affecting the I/O-complexity too much., i.e., it may be a function of $|V|$. We next show how an expansion analysis of this graph can be used to obtain the I/O-complexity lower bound for the corresponding algorithm.

2 The partition argument

Let $M$ be the size of the fast memory. Let $O$ be any total ordering of the vertices that respects the partial ordering of the CDAG $G$, i.e., all the edges are going up in the total order. This total ordering can be thought of as the actual order in which the computations are performed. Let $P$ be any partition of $V$ into segments $S_{1},S_{2},...$, so that a segment $S_{i}\in P$ is a subset of the vertices that are contiguous in the total ordering $O$.

Let $R_{S}$ and $W_{S}$ be the set of read and write operands, respectively (see Figure 1). Namely, $R_{S}$ is the set of vertices outside $S$ that have an edge going into $S$, and $W_{S}$ is the set of vertices in $S$ that have an edge going outside of $S$. Then the total I/O-complexity due to reads of AOs in $S$ is at least $|R_{S}|-M$, as at most $M$ of the needed $|R_{S}|$ operands are already in fast memory when the execution of the segment’s AOs starts. Similarly, $S$ causes at least $|W_{S}|-M$ actual write operations, as at most $M$ of the operands needed by other segments are left in the fast memory when the execution of the segment’s AOs ends. The total I/O-complexity is therefore bounded below byOne can think of this as a game: the first player orders the vertices. The second player partitions them into contiguous segments. The objective of the first player (e.g., a good programmer) is to order the vertices so that any consecutive partitioning by the second player leads to a small communication count.

3 Edge expansion and I/O-complexity

Consider a segment $S$ and its read and write operands $R_{S}$ and $W_{S}$ (see Figure 1). If the graph $G$ containing $S$ has $h(G)$ edge expansionThe direction of the edges does not matter much for the expansion-bandwidth argument: treating all edges as undirected changes the I/O-complexity estimate by a factor of 2 at most. For simplicity, we will treat $G$ as undirected., maximum degree $d$ and at least $2|S|$ vertices, then (using the definition of $h(G)$), we have

$|R_{S}|+|W_{S}|\geq\frac{1}{2}\cdot h(G)\cdot|S|$ .

We have $|E(S,V\setminus S)|\geq h(G)\cdot d\cdot|S|$. Either (at least) half of the edges $E(S,V\setminus S)$ touch $R_{S}$ or half of them touch $W_{S}$. As every vertex is of degree $d$, we have $|R_{S}|+|W_{S}|\geq\max\{|R_{S}|,|W_{S}|\}\geq\frac{1}{d}\cdot\frac{1}{2}\cdot|E(S,V\setminus S)|\geq h(G)\cdot|S|/2$.

Combining this with (6) and choosing to partition $V$ into $|V|/s$ segments of equal size $s$, we obtain: $IO\geq\max_{s}\frac{|V|}{s}\cdot\left(\frac{h(G)\cdot s}{2}-2M\right)=\Omega\left(|V|\cdot h(G)\right)$. In many cases $h(G)$ is too small to attain the desired I/O-complexity lower bound. Typically, $h(G)$ is a decreasing function in $|V(G)|$, namely the edge expansion deteriorates with the increase of the input size and with the running time of the corresponding algorithm. This is the case with matrix multiplication algorithms: the cubic, as well as the Strassen and “Strassen-like” algorithms. In such cases, it is better to consider the expansion of $G$ on small sets only: $IO\geq\max_{s}\frac{|V|}{s}\cdot\left(\frac{h_{s}(G)\cdot s}{2}-2M\right)$. ChoosingThe existence of a value $s$ that satisfies the condition is not always guaranteed. In the next section we confirm this for Strassen, for sufficiently large $|V(G)|$ (in particular, $|V(G)|$ has to be larger than $M$). Indeed this is the interesting case, as otherwise all computations can be performed inside the fast memory, with no communication, except for reading the input once. the minimal $s$ so that

In some cases, the computation graph $G$ does not fit this analysis: it may not be regular, it may have vertices of unbounded degree, or its edge expansion may be hard to analyze. In such cases, we may consider some subgraph $G^{\prime}$ of $G$ instead to obtain a lower bound on the I/O-complexity :

Let $G=(V,E)$ be a computation graph of an algorithm $Alg$. Let $G^{\prime}=(V^{\prime},E^{\prime})$ be a subgraph of $G$, i.e., $V^{\prime}\subseteq V$ and $E^{\prime}\subseteq E$. If $G^{\prime}$ is $d$-regular and $\alpha=\frac{|V^{\prime}|}{|V|}$, then the I/O-complexity of $Alg$ is

where $s$ is chosen so that ${h_{s}(G^{\prime})\cdot\alpha s\over 2}\geq 3M~{}$.

The correctness of this claim follows from Equations (7) and (8), and from the fact that at least an $\alpha/2$ fraction of the segments have at least $\frac{\alpha}{2}\cdot s$ of their vertices in $G^{\prime}$ (otherwise $V^{\prime}<\frac{\alpha}{2}\cdot V/s\cdot s+(1-\frac{\alpha}{2})\cdot V/s\cdot\frac{\alpha}{2}s<\alpha V$). We therefore have:

Let $Alg$ be an algorithm with $AO(N)$ arithmetic operations ($N$ being the total input size, $N=\Theta(n^{2})$ for matrix multiplication) and computation graph $G(N)=(V,E)$. Let $G^{\prime}(N)=(V^{\prime},E^{\prime})$ be a regular constant degree subgraph of $G$, with $\frac{|V^{\prime}|}{|V|}=\Theta(1)$. Then the I/O-complexity of $Alg$In Strassen’s algorithm, $N=2n^{2}$ is the number of input matrices elements and $T(N)=\Theta\left(n^{\omega_{0}}\right)=\Theta\left(N^{\omega_{0}/2}\right)$. $G^{\prime}$ is the graph $Dec_{k}C$ for $k=\lg M$, see Section 4 for the definition of $Dec_{k}C$. on a machine with fast memory of size $M$ is

As $AO(N)=\Theta(|V^{\prime}|)$ and $h_{s}(G^{\prime}(N))$ for $s=AO(M)$ is $\Theta(h(G^{\prime}(M)))$ (recall Claim 1) we obtain, equivalently,

Expansion Properties of Strassen’s Algorithm

Recall Strassen’s algorithm for matrix multiplication (see Algorithm 1 in Appendix A) and consider its computation graph (see Figure 2). Let $H_{i}$ be computation graph of Strassen’s algorithm for recursion of depth $i$, hence $H_{\lg n}$ corresponds to the computation for input matrices of size $n\times n$. $H_{\lg n}$ has the following structure:

Encode $A$: generate weighted sums of elements of $A$ (this corresponds to the left factors of lines 5-11 of the algorithm).

Similarly encode $B$ (this corresponds to the right factors of lines 5-11 of the algorithm).

Then multiply the encodings of $A$ and $B$ element-wise (this corresponds to line 2 of the algorithm).

Finally, decode $C$, by taking weighted sums of the products (this corresponds to lines 12-15 of the algorithm).

$Dec_{1}C$ is presented, for simplicity, with vertices of in-degree larger than two (but constant). A vertex of degree larger than two, in fact, represents a full binary (not necessarily balanced) tree. Note that replacing these high in-degree vertices with trees changes the edge expansion of the graph by a constant factor at most (as this graph is of constant size, and connected). Moreover, there is no change in the number of input and output vertices. Therefore the arguments in the following proof of Lemma 4.2 still hold.

Assume w.l.o.g. that $n$ is an integer power of $2$. Denote by $Enc_{\lg n}A$ the part of $H_{\lg n}$ that corresponds to the encoding of matrix $A$. Similarly, $Enc_{\lg n}B$, and $Dec_{\lg n}C$ correspond to the parts of $H_{\lg n}$ that compute the encoding of $B$ and the decoding of $C$, respectively.

We next construct the computation graph $H_{i+1}$ by constructing $Dec_{i+1}C$ (from $Dec_{i}C$ and $Dec_{1}C$) and similarly constructing $Enc_{i+1}A$ and $Enc_{i+1}B$, then composing the three parts together.

Identify the $4\cdot 7^{i}$ output vertices of the copies of $Dec_{1}C$ with the $4\cdot 7^{i}$ input vertices of the copies of $Dec_{i}C$:

Recall that each $Dec_{1}C$ has four output vertices.

The first output vertex of the $7^{i}$ $Dec_{1}C$ graphs are identified with the $7^{i}$ input vertices of the first copy of $Dec_{i}C$.

The second output vertex of the $7^{i}$ $Dec_{1}C$ graphs are identified with the $7^{i}$ input vertices of the second copy of $Dec_{i}C$. And so on.

We make sure that the $j$th input vertex of a copy of $Dec_{i}C$ is identified with an output vertex of the $j$th copy of $Dec_{1}C$.

We similarly obtain $Enc_{i+1}A$ from $Enc_{i}A$ and $Enc_{1}A$,

and $Enc_{i+1}B$ from $Enc_{i}B$ and $Enc_{1}B$.

For every $i$, $H_{i}$ is obtained by connecting edges from the $j$th output vertices of $Enc_{i}A$ and $Enc_{i}B$ to the $j$th input vertex of $Dec_{i}C$.

This completes the construction. Let us note some properties of these graphs.

The graph $Dec_{1}C$ has no vertices which are both input and output. As all out-degrees are at most 4 and all in degree are at most 2 (Recall Comment 4.1) we have:

All vertices of $Dec_{\lg n}C$ are of degree at most $6$.

However, $Enc_{1}A$ and $Enc_{1}B$ have vertices which are both input and output (e.g., $A_{11}$), therefore $Enc_{\lg n}A$ and $Enc_{\lg n}B$ have vertices of out-degree $\Theta(\lg n)$. All in-degrees are at most $2$, as an arithmetic operation has at most two inputs.

As $H_{\lg n}$ contains vertices of large degrees, it is easier to consider $Dec_{\lg n}C$: it contains only vertices of constant bounded degree, yet at least one third of the vertices of $H_{\lg n}$ are in it.

(Main lemma) The edge expansion of $Dec_{k}C$ is

The proof follows below, but first note that it suffices to deduce the expansion of $Dec_{\lg n}C$ on small sets. Assume w.l.o.g. that $n$ is an integer power of $\sqrt{M}$.We may assume this, as we are dealing with a lower bound here, so it suffices to prove the assertion for an infinite number of $n$’s. Alternatively, in the following decomposition argument, we leave out a few of the top or bottom levels of vertices of $Dec_{\lg n}C$, so that $n$ is an integer power of $\sqrt{M}$ and so that at most $|S|/2$ vertices of $S$ are cut off. Then $Dec_{\lg n}C$ can be split into edge-disjoint copies of $Dec_{\frac{1}{2}\lg M}C$. Using Claim 1, we thus have:

$s\cdot h_{s}(Dec_{\lg n}C)\geq 3M$ for $s=9\cdot M^{\lg 7/2}$.

As $Dec_{\lg n}C$ contains $\alpha=\frac{1}{3}$ of the vertices of $H_{\lg n}$, Lemma 3.2 now yields Main Theorem 1.1. Note that $Dec_{\lg n}C$ has no input vertices, so no restriction on input replication is needed.

1.2 Combinatorial Estimation of the Expansion

Let $G_{k}=(V,E)$ be $Dec_{k}C$, and let $S\subseteq V,|S|\leq|V|/2$. We next show that $|E(S,V\setminus S)|\geq c\cdot d\cdot|S|\cdot\left(\frac{4}{7}\right)^{k}$, where $c$ is some universal constant, and $d$ is the constant degree of $Dec_{k}C$ (after adding loops to make it regular).

The proof works as follows. Recall that $G_{k}$ is a layered graph (with layers corresponding to recursion steps), so all edges (excluding loops) connect between consecutive levels of vertices. We argue (in Claim 8) that each level of $G_{k}$ contains about the same fraction of $S$ vertices, or else we have many edges leaving $S$. We also observe (in Fact 9) that such homogeneity (of a fraction of $S$ vertices) does not hold between distinct parts of the lowest level, or, again, we have many edges leaving $S$. We then show that the homogeneity between levels, combined with the heterogeneity of the lowest level, guarantees that there are many edges leaving $S$.

Let $l_{i}$ be the $i$th level of vertices of $G_{k}$, so $4^{k}=|l_{1}|<|l_{2}|<\cdots<|l_{i}|=4^{k-i+1}7^{i-1}<\cdots<|l_{k+1}|=7^{k}$. Let $S_{i}\equiv S\cap l_{i}$. Let $\sigma=\frac{|S|}{|V|}$ be the fractional size of $S$ and $\sigma_{i}=\frac{|S_{i}|}{|l_{i}|}$ be the fractional size of $S$ at level $i$. Due to averaging, we observe the following:

There exist $i$ and $i^{\prime}$ such that $\sigma_{i}\leq\sigma\leq\sigma_{i^{\prime}}$.

so $\frac{3}{7}\leq\frac{|l_{k+1}|}{|V|}\leq\frac{3}{7}\cdot\frac{1}{1-\left(\frac{4}{7}\right)^{k+2}}$, and $\frac{3}{7}\cdot\left(\frac{4}{7}\right)^{k}\leq\frac{|l_{1}|}{|V|}\leq\frac{3}{7}\cdot\left(\frac{4}{7}\right)^{k}\cdot\frac{1}{1-\left(\frac{4}{7}\right)^{k+2}}.$

There exists $c^{\prime}=c^{\prime}(G_{1})$ so that $|E(S,V\setminus S)\cap E(l_{i},l_{i+1})|\geq c^{\prime}\cdot d\cdot|\delta_{i}|\cdot|l_{i}|$.

Let $G^{\prime}$ be a $G_{1}$ component connecting $l_{i}$ with $l_{i+1}$ (so it has four vertices in $l_{i}$ and seven in $l_{i+1}$). $G^{\prime}$ has no edges in $E(S,V\setminus S)$ if all or none of its vertices are in $S$. Otherwise, as $G^{\prime}$ is connected, it contributes at least one edge to $E(S,V\setminus S)$. The number of such $G_{1}$ components with all their vertices in $S$ is at most $\min\{\sigma_{i},\sigma_{i+1}\}\cdot\frac{|l_{i}|}{4}$. Therefore, there are at least $|\sigma_{i}-\sigma_{i+1}|\cdot\frac{|l_{i}|}{4}$ $G_{1}$ components with at least one vertex in $S$ and one vertex that is not.

If there exists $i$ so that $\frac{|\sigma-\sigma_{i}|}{\sigma}\geq\frac{1}{10}$, then

where $c>0$ is some constant depending on $G_{1}$ only.

Assume that there exists $j$ so that $\frac{|\sigma-\sigma_{j}|}{\sigma}\geq\frac{1}{10}$. Let $\delta_{i}\equiv\sigma_{i+1}-\sigma_{i}$. By Claim 7, we have

By the initial assumption, there exists $j$ so that $\frac{|\sigma-\sigma_{j}|}{\sigma}\geq\frac{1}{10}$, therefore $\max_{i}\sigma_{i}-\min_{i}\sigma_{i}\geq\frac{\sigma}{10}$, then

for any $c\leq\frac{c^{\prime}}{10}\cdot\frac{3}{7}$.

Let $T_{k}$ be a tree corresponding to the recursive construction of $G_{k}$ in the following way (see Figure 3): $T_{k}$ is a tree of height $k+1$, where each internal node has four children. The root $r$ of $T_{k}$ corresponds to $l_{k+1}$ (the largest level of $G_{k}$). The four children of $r$ correspond to the largest levels of the four graphs that one can obtain by removing the level of vertices $l_{k+1}$ from $G_{k}$. And so on. For every node $u$ of $T_{k}$, denote by $V_{u}$ the set of vertices in $G_{k}$ corresponding to $u$. We thus have $|V_{r}|=7^{k}$ where $r$ is the root of $T_{k}$, $|V_{u}|=7^{k-1}$ for each node $u$ that is a child of $r$; and in general we have $4^{i}$ tree nodes $u$ corresponding to a set of size $|V_{u}|=7^{k-i+1}$. Each leaf $l$ corresponds to a set of size $1$.

For a tree node $u$, let us define $\rho_{u}=\frac{|S\cap V_{u}|}{|V_{u}|}$ to be the fraction of $S$ nodes in $V_{u}$, and $\delta_{u}=|\rho_{u}-\rho_{p(u)}|$, where $p(u)$ is the parent of $u$ (for the root $r$ we let $p(r)=r$). We let $t_{i}$ be the $i$th level of $T_{k}$, counting from the bottom, so $t_{k+1}$ is the root and $t_{1}$ are the leaves.

As $V_{r}=l_{k+1}$ we have $\rho_{r}=\sigma_{k+1}$. For a tree leaf $u\in t_{1}$, we have $|V_{u}|=1$. Therefore $\rho_{u}\in\{0,1\}$. The number of vertices $u$ in $t_{1}$ with $\rho_{u}=1$ is $\sigma_{1}\cdot|l_{1}|$.

Let $u_{0}$ be an internal tree node, and let $u_{1},u_{2},u_{3},u_{4}$ be its four children. Then

where $c^{\prime\prime}=c^{\prime\prime}(G_{1})$.

The proof follows that of Claim 7. Let $G^{\prime}$ be a $G_{1}$ component connecting $V_{u_{0}}$ with $\bigcup_{i\in}V_{u_{i}}$ (so it has seven vertices in $V_{u_{0}}$ and one in each of $V_{u_{1}}$,$V_{u_{2}}$,$V_{u_{3}}$,$V_{u_{4}}$). $G^{\prime}$ has no edges in $E(S,V\setminus S)$ if all or none of its vertices are in $S$. Otherwise, as $G^{\prime}$ is connected, it contributes at least one edge to $E(S,V\setminus S)$. The number of $G_{1}$ components with all their vertices in $S$ is at most $\min\{\rho_{u_{0}},\rho_{u_{1}},\rho_{u_{2}},\rho_{u_{3}},\rho_{u_{4}}\}\cdot\frac{|V_{u_{1}}|}{4}$. Therefore, there are at least $\max_{i\in}\{|\rho_{u_{0}}-\rho_{u_{i}}|\}\cdot\frac{|V_{u_{1}}|}{4}\geq\frac{1}{16}\cdot\sum_{i\in}|\rho_{u_{i}}-\rho_{u_{0}}|\cdot|V_{u_{i}}|$ $G_{1}$ components with at least one vertex in $S$ and one vertex that is not.

Other Algorithms

We now discuss the applicability of our approach to other algorithms, starting with other fast matrix multiplication algorithms.

A “Strassen-like” algorithm has a recursive structure that utilizes a base case: multiplying two $n_{0}$-by-$n_{0}$ matrices using $m(n_{0})$ multiplications. Given two matrices of size $n$-by-$n$, it splits them into $n_{0}^{2}$ blocks (each of size $\frac{n}{n_{0}}$-by-$\frac{n}{n_{0}}$), and works blockwise, according to the base case algorithm. Additions (and subtractions) in the base case are interpreted as additions (and subtractions) of blocks. These are performed element-wise. Multiplications in the base case are interpreted as multiplications of blocks. These are performed by recursively calling the algorithm. The arithmetic count of the algorithm is then $T(n)=m(n_{0})\cdot T\left(\frac{n}{n_{0}}\right)+O(n^{2})$, so $T(n)=\Theta(n^{\omega_{0}})$ where $\omega_{0}=\log_{n_{0}}m(n_{0})$.

This is the structure of all the fast matrix multiplication algorithms that were obtained since Strassen’s [Pan (1980), Bini (1980), Schönhage (1981), Romani (1982), Coppersmith and Winograd (1982), Strassen (1987), Coppersmith and Winograd (1987)], (see [Bűrgisser et al. (1997)] for discussion of these algorithms), as well as [Cohn et al. (2005)], where the base case utilizes a novel group-theoretic approach. In fact, any fast matrix multiplication algorithm can be converted into this form [Raz (2003)], and can even be made numerically stable while preserving this form [Demmel, Dumitriu, Holtz, and Kleinberg, 2007].

For our technique to work, we further demand that the $Dec_{1}C$ part of the computation graph is a connected graph, in order to be “Strassen-like” (this was assumed in the proof of Claim 7). Thus the “Strassen-like” class includes Winograd’s variant of Strassen’s algorithm [Winograd (1971)], which uses 15 additions rather than 18, but not the cubic algorithm, where $Dec_{1}C$ is composed of four disconnected graphs (corresponding to the four outputs). We conjecture that $Dec_{1}C$ is indeed connected for all existing fast matrix-multiplication algorithms. We note that the demand of connectivity of $Dec_{1}C$ may be waved in some cases (see [Ballard et al. (2011e)]).

1.2 The communication costs of “Strassen-like” algorithms

To prove Theorem 1.3, which generalizes the I/O-complexity lower bound of Strassen’s algorithm (Theorem 1.1) to all “Strassen-like” algorithms, we note the following: The entire proof of Theorem 1.1, and in particular, the computations in the proof of Lemma 4.2, hold for any “Strassen-like” algorithm, where we plug in $n_{0}^{2},m(n_{0})$, and $\frac{n_{0}}{m(n_{0})}$ instead of $4,7$, and $\frac{4}{7}$. For bounding the asymptotic I/O-complexity , we do not care about the number of internal vertices of $Dec_{1}C$; we need only to know that $Dec_{1}C$ is connected (this critical technical assumption is used in the proof of Claim 7), and to know the sizes $n_{0}$ and $m(n_{0})$. The only nontrivial adjustment is to show the equivalent of Fact 4: that the graph $Dec_{\log n}C$ is of bounded degree.

The $Dec_{\log n}C$ graph of any “Strassen-like” algorithm is of degree bounded by a constant.

If the set of input vertices of $Dec_{1}C$, and the set of its output vertices are disjoint, then $Dec_{\log n}C$ is of constant bounded degree (its maximal degree is at most twice the largest degree of $Dec_{1}C$).

Assume (towards contradiction) that the base graph $Dec_{1}C$ has an input vertex which is also an output vertex. An output vertex represents the inner product of two $n_{0}$-long vectors, i.e., the corresponding row-vector of $A$ and column vector of $B$. The corresponding bilinear polynomial is irreducible. This is a contradiction, since an input vertex represents the multiplication of a (weighted) sum of elements of $A$ with a (weighted) sum of elements of $B$.

2 Uniform, Non-stationary Fast Matrix Multiplication Algorithms

Another class of matrix multiplication algorithms, the uniform, non-stationary algorithms, allows mixing of schemes of the previous (“Strassen-like”) class. In each recursive level, a different scheme may be used. The CDAG has a repeating structure inside one level, but the structure may differ between two distinct levels. This class includes, for example, algorithms that optimize for input sizes (for sizes that are not an integer power of a constant integer). The class also includes algorithms that cut the recursion off at some point, and then switch to the classical algorithm. For these and other implementation issues, see [Douglas et al. (1994), Huss-Lederman et al. (1996)] (sequential model) and [Desprez and Suter (2004)] (parallel model). The I/O-complexity lower bound generalizes to this class, and will appear in a separate note [Ballard et al. (2011e)].

3 Non-uniform, Non-stationary Fast Matrix Multiplication Algorithms

A third class, the non-uniform, non-stationary algorithms, allows recursive calls to have different structure, even when they refer to multiplication of matrices in the same recursive level. It is not clear how to analyze the expansion of the CDAG of an algorithm in the third class, although we are not aware of any algorithms in this class. Such an analysis, applied to the base case of [Cohn et al. (2005)], may improve the I/O-complexity lower bound for fast matrix multiplication by a (large) constant.

4 Multiplying Rectangular Matrices

Multiplication of rectangular matrices have seen a series of increasingly fast algorithms culminating in Coppersmith’s algorithm [Coppersmith (1997)]. It is possible to extend our approach and obtain the first lower bounds on the communication costs for these algorithms, and show that in some cases they are attainable, and therefore optimal [Ballard et al. (2011d)].

5 Other Algorithms

Fast matrix multiplication algorithms are basic building blocks in many fast algorithms in linear algebra, such as algorithms for LU, QR, and solving the Sylvester equation [Demmel, Dumitriu, and Holtz, 2007]. Therefore, I/O-complexity lower bounds for these algorithms can be derived from our lower bounds for fast matrix multiplication algorithms [Ballard et al. (2011a)]. For example, a lower bound on LU (or QR, etc.) follows when the fast matrix multiplication algorithm is called by the LU algorithm on sufficiently large subblocks of the matrix. This is the case in the algorithms of [Demmel, Dumitriu, and Holtz, 2007], and we can then deduce matching lower and upper bounds [Ballard et al. (2011a)].

Conclusions and Open Problems

We obtained a tight lower bound for the I/O-complexity of Strassen’s and “Strassen-like” fast matrix multiplication algorithms. These bounds are optimal for the sequential model with two memory levels and with memory hierarchy. The lower bounds extend to the parallel model and other models. Recently these bounds were attained (up to an $O(\log p)$ factor) by new parallel implementations, for Strassen’s algorithm and for “Strassen-like” algorithms [Ballard et al. (2011)].

Some (parallel) algorithms require very little, up to a constant factor extra memory beyond what is necessary to keep the input and output. These are sometimes called linear space algorithms. One class of such algorithms are the “2D” algorithms for classical matrix multiplication, that use two-dimensional grid of processors. Here we allow $M=\Theta\left(\frac{n^{2}}{p}\right)$ local memory use (recall that $p$ is the number of processors, and $n$ the dimension of the matrices), thus no replication of the input matrices is allowed [Cannon (1969)].

If the underlying grid of $p$ processors is a three-dimensional mesh, and the available memory per processor is larger by a factor of $p^{\frac{1}{3}}$ than the minimum necessary to store the input and output matrices, then a “3D” algorithm can be used (see [Dekel et al. (1981), Aggarwal et al. (1990), Agarwal et al. (1995), McColl and Tiskin (1999)]). These “3D” algorithms can reduce the communication cost by a factor of $p^{1/6}$, down to $\Theta\left(\frac{n^{2}}{p^{\frac{2}{3}}}\right)$ [Aggarwal et al. (1990)], attaining the lower bounds [Irony et al. (2004), Ballard et al. (2011c)] that take into account any amount of replication.

Recently, Demmel and Solomonik [Solomonik and Demmel (2011)] showed how to combine these two extremes into one algorithm (named “2.5D”) and obtained a communication efficient implementation for classical matrix multiplication, for local memory size $M=\Theta\left(c\cdot\frac{n^{2}}{p}\right)$ for any $1\leq c\leq p^{\frac{1}{3}}$. See Table 1.

Using Corollaries 1.2 and 1.4, and plugging in $M=\Theta\left(\frac{n^{2}}{p}\right)$, $M=\Theta\left(\frac{n^{2}}{p^{\frac{2}{3}}}\right)$ , and $M=\Theta\left(c\cdot\frac{n^{2}}{p}\right)$ we obtain corresponding lower bounds for “Strassen-like” algorithm with various restriction local memory sizes. These were recently attained by a new parallel implementation for Strassen and “Strassen-like” algorithms [Ballard et al. (2011)] (see Table 1). Interestingly, the numerators here do not depend on $\omega_{0}$. Thus, an improvement of $\omega_{0}$ (the exponent of the arithmetic cost of the algorithm) affects only the power of $p$ in the denominator.

2 Recursive Implementations

In some cases, the simplest recursive implementation of an algorithm turns out to be communication-optimal (e.g., in the cases of matrix multiplication [Frigo et al. (1999)] and Cholesky decomposition [Ahmed and Pingali (2000), Ballard et al. (2010)], but not in the case of LU decomposition [Toledo (1997)], which is bandwidth optimal but not latency optimal). This leads to the question: when is the communication-optimality of an algorithm determined by the expansion properties of the corresponding computation graphs? In this work we showed that such is the case for “Strassen-like” fast matrix multiplication algorithms.

3 Other Hardware

It is of great interest to construct new models general enough to capture the rich and evolving design space of current and predicted future computers. Such models can be homogeneous, consisting of many layers, where the components of each layer are the same (e.g., a supercomputer with many identical multi-core chips on a board, many identical boards in a rack, many identical racks, and many identical levels of associated memory hierarchy); or heterogeneous, with components with different properties residing on the same level (e.g., CPUs alongside GPUs, where the latter can do some computations very quickly, but are much slower to communicate with).

Some experience has been acquired with such systems (see the MAGMA project [(7)], and also [Volkov and Demmel (2008)] for using GPU assisted linear algebra computation ). A first step in analyzing such systems has been recently introduced by Ballard, Demmel, and Gearhart [Ballard et al. (2011)], where they modeled heterogenous shared memory architectures, such as mixed GPU/CPU architecture, and obtained tight lower and upper bounds for $O(n^{3})$ matrix multiplication.

Note that we can similarly generalize Corollaries 1.2 and 1.4 to other models, such as the heterogenous model and shared memory model. The reduction is achieved by observing the communication of a single processor.

However, there is currently no systematic theoretic way of obtaining upper and lower bounds for arbitrary hardware models. Expanding such results to other architectures and algorithmic techniques is a challenging goal. For example, recursive algorithms tend to be cache oblivious and communication optimal for the sequential hierarchy model. Finding an equivalent technique that would work for an arbitrary architecture is a fundamental open problem.

We thank Eran Rom, Edgar Solomonik, and Chris Umans for helpful discussions.

References

Appendix A Strassen’s Fast Matrix Multiplication Algorithm

Strassen’s original algorithm follows [Strassen (1969)]. See [Winograd (1971)] for Winograd’s variant, which reduces the number of additions.