Low-Bandwidth Matrix Multiplication: Faster Algorithms and More General Forms of Sparsity
Chetan Gupta Janne H. Korhonen Jan Studený Jukka Suomela Hossein Vahidi

Casimir energy on the sphere and 6D CFT trace anomaly
    R Aros^1, F Bugini^2, DE Díaz^3 and C Núñez-Barra^4
    
=======================================================




    In prior work, Gupta et al. (SPAA 2022) presented a distributed algorithm for multiplying sparse n × n matrices, using n computers. They assumed that the input matrices are uniformly sparse—there are at most d non-zeros in each row and column—and the task is to compute a uniformly sparse part of the product matrix. Initially each computer knows one row of each input matrix, and eventually each computer needs to know one row of the product matrix. In each communication round each computer can send and receive one O(log n)-bit message. Their algorithm solves this task in O(d^1.907) rounds, while the trivial bound is O(d^2).

    We improve on the prior work in two dimensions: First, we show that we can solve the same task faster, in only O(d^1.832) rounds. Second, we explore what happens when matrices are not uniformly sparse. We consider the following alternative notions of sparsity: row-sparse matrices (at most d non-zeros per row), column-sparse matrices, matrices with bounded degeneracy (we can recursively delete a row or column with at most d non-zeros), average-sparse matrices (at most dn non-zeros in total), and general matrices.

    We show that we can still compute X = AB in O(d^1.832) rounds even if one of the three matrices (A, B, or X) is average-sparse instead of uniformly sparse. We present algorithms that handle a much broader range of sparsity in O(d^2 + log n) rounds, and present conditional hardness results that put limits on further improvements and generalizations.


Casimir energy on the sphere and 6D CFT trace anomaly
    R Aros^1, F Bugini^2, DE Díaz^3 and C Núñez-Barra^4
    
=======================================================





§ INTRODUCTION


In this work we study sparse matrix multiplication in a distributed setting, for the low-bandwidth model (which is closely related to the node-capacitated clique model). Our work improves over the prior work by <cit.> in two ways:

    
  * we present a faster algorithm, and
    
  * we study more general forms of sparsity.




 §.§ Setting and prior work


The task is to compute the matrix product X = AB for n × n matrices using a network of n computers. Initially each computer holds its own part of A and B. Computation proceeds in rounds, and in one round each computer can send one O(log n)-bit message to another computer and receive one message from another computer; we will assume that the elements of A, B, and X fit in one message. Eventually each computer has to report its own part of the product matrix X. How many communication rounds are needed to solve the task?

The trivial solution for dense matrices takes O(n^2) rounds: everyone sends all information to computer number 1, which solves the task locally and then distributes the solution to other computers. However, we can do better: for matrix multiplication over semirings there is an algorithm that runs in O(n^) rounds, and for matrix multiplication over fields there is an algorithm that runs in O(n^2-2/ω) rounds, where ω is the exponent of centralized matrix multiplication <cit.>. By plugging in the latest value of ω < 2.371552 <cit.>, we obtain O(n^) rounds.

The key question is how much better we can do when we multiply sparse matrices. Let us first look at the case of uniformly sparse matrices, with at most d nonzero elements in each row and column. We also assume that we are only interested in a uniformly sparse part of the product matrix X. Now it is natural to assume that each computer initially holds one row of A and one row of B, and it needs to know one row of X. We assume the supported setting: the sparsity structures of A, B, and X are known in advance, while the values of the nonzero elements are revealed at run time. There are two algorithms from prior work that are applicable in this setting: for moderately large values of d we can use the algorithm by <cit.>, which runs in O(dn^1/3) rounds, while for small values of d the fastest algorithm is due to <cit.>, and the round complexity is O(d^) for semirings and O(d^) for fields; see <ref> for an overview.





 §.§ Contribution 1: faster algorithm


Our first contribution is improvements in the running time: we design a faster algorithm that solves the case of semirings in O(d^) rounds and the case of fields in O(d^) rounds.

While there are no non-trivial unconditional lower bounds, we point out that O(d^) for semirings or O(d^) for fields would imply major breakthroughs for dense matrix multiplication (by simply plugging in d = n). In the following figure, we illustrate the progress we make towards these milestones:


    
    < g r a p h i c s >




The algorithm from prior work <cit.> is based on the idea of processing triangles. Here a triangle is a triple {i,j,k} such that A_ij and B_jk are nonzeros and X_ik is one of the elements of interest, and we say that we process the triangle if we add the product A_ijB_jk to the sum X_ik. Once all such triangles are processed, we have computed all elements of interest in the product matrix.

There are two strategies that <cit.> uses for processing triangles:

    
  * If there are many triangles in total, we can find a dense cluster of them, we can interpret such a cluster as a tiny instance of dense matrix multiplication, and we can apply dense matrix multiplication algorithms to batch-process such clusters.
    
  * If there are few triangles in total, we can afford to process the triangles one by one.


Our improvements are in the second phase. Trivially, if we have a uniformly sparse instance, there are at most d^2n triangles in total, and we can process them in a naive manner in O(d^2) rounds. The second phase of <cit.> shows that if the number of triangles is d^2- n, one can process them in O(d^2-/2) rounds. Here the key obstacle is the expression /2 in the exponent: even if we eliminate many triangles in the first phase, the second phase is still relatively expensive. We design a new algorithm that is able to process d^2- n triangles in O(d^2-) rounds. By combining this with the general strategy of <cit.> and optimizing the trade-off between using time in the first phase vs. using time in the second phase, we obtain the new round complexities O(d^) and O(d^)—the difference between semirings and fields is in the dense matrix multiplication routine that we use in the first phase.



 §.§ Contribution 2: beyond uniform sparsity


So far we have followed the lead of <cit.> and discussed uniformly sparse matrices (i.e., each row and column has at most d nonzeros). This is a rather restrictive notion of sparsity—can we extend prior work to more general notions of sparsity?

Perhaps the most natural candidate to consider would be average sparsity: we would simply assume that the total number of nonzeros in A and B and the total number of elements of interest in X is at most dn. Unfortunately, it turns out that this notion of sparsity is way too general, and we cannot expect fast algorithms for this case: in <ref> we show that if we could solve average-sparse matrix multiplication for n× n matrices using n computers in time that is independent of n or only mildly depends on n (say, polylogarithmic in n), it would imply major breakthroughs in our understanding of dense matrix multiplication.

In summary, in order to have any real hope of pushing the complexity of sparse matrix multiplication down to O(d^) or even something more modest like O(d^2 + log n), we need to explore some intermediate notions of sparsity. In this work we make use of the following notions of sparsity:

    
  * (d) = uniformly sparse: at most d nonzeros per row and column.
    
  * (d) = row-sparse: at most d nonzeros per row.
    
  * (d) = column-sparse: at most d nonzeros per column.
    
  * (d) = bounded degeneracy: we can recursively eliminate the matrix so that at each step we delete a row or column with at most d nonzeros.
    
  * (d) = average-sparse: at most dn nonzeros in total.
    
  *  = general matrices.

We will omit (d) when it is clear from the context, and we will use e.g. expressions such as × = to refer to the task of computing AB = X such that A ∈(d), B ∈(d), and X̂∈(d); here we write X̂ for the matrix that indicates which elements of the product X we are interested in.

Family  may at first look rather unusual in the context of linear algebra. However, if we interpret a matrix A as a bipartite graph G (with an edge {i,j} whenever A_ij is nonzero), then A ∈(d) corresponds to the familiar graph-theoretic notion that G is d-degenerate. This is a widely used notion of sparsity in the graph-theoretic setting, and also closely connected with other notions of sparsity such as bounded arboricity.

We also point out that any matrix A ∈(d) can be written as a sum A = X + Y such that X ∈ RS(d) and Y ∈ CS(d); to see this, eliminate A by deleting sparse rows or sparse columns, and put each sparse row in X and put each sparse column in Y. In particular, matrix multiplication × decomposes into operations of the form ×, ×, ×, and ×. In summary, we have

    ⊆{
    }⊆⊆⊆.


As we will see, all of our algorithmic results are symmetric w.r.t. the three matrices. To help present such results, we will introduce the following shorthand notation. Let , , and  be families of matrices. We write [::] to refer to the following six operation: computing × =, × =, × =, × =, × =, and × =.

Now equipped with this notation, we can rephrase the main result from the prior work <cit.>: they show that [ :  : ] can be solved in O(d^) rounds for semirings and in O(d^) rounds for fields. On the other hand, in <ref> we argue why it is unlikely to extend this all the way to [ :  : ]. The key question is what happens between the two extremes,  and ? In particular, could we solve [::] efficiently?

Our answer to this question consists of four parts:

    
  * We generalize prior work and show that we can tolerate one average-sparse matrix: [ :  : ] can be computed in O(d^) rounds for semirings and in O(d^) rounds for fields. We emphasize that here the running time is still independent of n.
    
  * As soon as two of the matrices are  or more general, the task becomes at least as hard as computing the sum of values distributed across n computers, or broadcasting a single value to n computers. Hence, it is natural to expect an additional O(log n) term in the running time.
    
  * We show that many combinations all the way to [ :  : ] and [ :  : ] can be solved in O(d^2 + log n) rounds, for both semirings and fields.
    
  * Finally, we provide evidence that this is essentially as general as what we can hope for.

The results are summarized in <ref>.





 §.§ Related work and applications


Our model of computing is the low-bandwidth model; this is closely related to the model known as node-capacitated clique or node-congested clique in the literature <cit.>—in the low-bandwidth model each computer can send and receive one message, while the node-capacitated clique model is usually defined so that each computer can send and receive O(log n) messages per round. Both the low-bandwidth model and the node-capacitated clique model can be interpreted as variants of the congested clique model <cit.>. Indeed, any algorithm that runs in T(n) rounds in the congested clique model can be simulated in nT(n) rounds in the low-bandwidth model, and for many problems (such as the dense matrix multiplication) such a simulation also results in the fastest known algorithms for the low-bandwidth model. Both the low-bandwidth model and the congested clique model can be interpreted as special cases of the classic bulk synchronous parallel model <cit.>.

Matrix multiplication in this context has been primarily studied in the congested clique model <cit.>. However, as observed by <cit.>, the congested clique model is poorly-suited for the study of sparse matrix multiplication: one can only conclude that for sufficiently sparse matrices, the problem is solvable in O(1) rounds, and one cannot explore more fine-grained differences between different algorithms.

Similar to <cit.>, we work in the supported version of the low-bandwidth model. In general, the supported model <cit.> refers to a setting in which the structure of the input is known in advance (we can do arbitrary preprocessing based on the structure of the input), while the specific instance is revealed at run time. In the case of graph problems, we can see an unweighted graph G in advance, while the adversary reveals a weighted version G' of G (a special case being a 0/1-weighted graph, which can be interpreted as a subgraph of G). For matrix multiplication, we know the sparsity structure in advance (i.e., which elements are potentially nonzero), while the adversary reveals the concrete values at run time.

One of the main applications that we have in mind here is triangle detection, which has been extensively studied in the distributed setting <cit.>. If we can multiply matrices, we can also easily detect triangles in a graph. Moreover, [::] corresponds to triangle detection in a bounded-degree graph, while e.g. [::] corresponds to triangle detection in a sparse graph.



 §.§ Open questions for future work


In this work we initiated the study of different notions of sparsity in the context of low-bandwidth matrix multiplication. This work also gives rise to a number of questions for future work.

First, our running times—O(d^) rounds for semirings and O(d^) for fields—are still far from the conditional lower bounds Ω(d^) and Ω(d^), respectively. Pushing the running time further down is a major challenge for future work. If we follow the two-phase approach of <cit.>, our present work essentially makes the second phase optimal; the main hope for future improvements now lies in the first phase of the algorithm.

Our algorithms are designed for the supported model. Eliminating the knowledge of the support (i.e., knowledge of the sparsity structure) is a major challenge for future work.

Our hardness results are conditional. Turning these into unconditional lower bounds is yet another question for future work.

Finally, there are some gaps in <ref>. For example, we do not have conditional hardness result for [::], and also <ref> do not cover all permutations of the matrix families.



§ PRELIMINARIES


We work in the low-bandwidth model. There are n computers. Initially each computer holds its own part of A and B, and eventually each computer has to report its own part of X. When we study sparse matrices with at most dn nonzeros, we assume that each computer holds at most d elements, while when we study dense matrices, we assume that each computer holds at most n elements.

For our algorithms, it does not matter how the input and output is distributed among the computers—with an additional O(d) time we can permute the input and output as appropriate. Our conditional lower bounds hold for any fixed distribution of input and output.



 §.§ Supported model and indicator matrices


When we study matrix multiplication in the supported model, we assume that we know in advance indicator matrices Â, B̂, and X̂ that encode the structure of our instance: Â_ij = 0 implies A_ij = 0, B̂_jk = 0 implies B_jk = 0, and X̂_ik = 0 indicates that we do not need to compute the value of X_ik. When we make assumptions on sparsity, our assumptions refer to the sparsity of Â, B̂, and X̂. For example, when we study × =, we assume that Â∈(d), B̂∈(d), and X̂∈(d). It then also follows that A ∈(d) and B ∈(d).



 §.§ Tripartite graph and triangles


It will be convenient to assume that our matrices are indexed with indices from three disjoint sets I, J, and K, each of size n. For example, elements of matrix A are indexed with A_ij where i ∈ I and j ∈ J. Following <cit.>, we write  for the set of all triangles, which are triples {i,j,k} with i ∈ I, j ∈ J, and k ∈ K such that Â_ij 0, B̂_jk 0, and X̂_ik 0.

For a collection of triangles , we write G() for the tripartite graph G() = (V,E), where the set of nodes is V = I ∪ J ∪ K and there is an edge {u,v}∈ E if there is a triangle T ∈ with {u,v}⊆ T. For a set of nodes U ⊆ V, we write [U] = { T ∈ : T ⊆ U } for the set of triangles induced by U.

We say that we process a triangle {i,j,k} if we have added the product A_ijB_jk to the sum X_ik. The key observation is that processing all triangles is exactly equivalent to computing all values of interest in the product matrix X = AB.



 §.§ Clusters and clusterings


We say that U ⊆ V is a cluster if U = I' ∪ J' ∪ K' for I' ⊆ I, J' ⊆ J, K' ⊆ K, |I'| = d, |J'| = d, and |K'| = d. We say that a collection of triangles  is clustered if there are disjoint clusters U_1, …, U_k such that = [U_1] ∪…∪[U_k].

A clustered collection of triangles can be processed efficiently by applying an algorithm for dense matrix multiplication in parallel to each cluster <cit.>. By applying the strategy of <cit.> and the latest value of ω from <cit.>, we obtain:

    A clustered instance of matrix multiplication can be solved in O(d^4/3) rounds over semirings, and in O(d^1.156671) rounds over fields.





§ HANDLING FEW TRIANGLES FAST


We will start by a general result that will form the foundation for all of our upper-bound results. <cit.> observed that if we have a balanced instance in which each node touches at most κ triangles, we can solve matrix multiplication in O(κ) rounds. However, the key challenge is what to do in unbalanced instances: we have at most κ n triangles, but some nodes can touch many more than κ triangles. The following lemma shows how to handle also such instances efficiently:


    Let  be a set of triangles with || ≤κ n, and assume that for each pair of nodes u,v ∈ V there are at most m triangles T ∈ with {u,v}∈ T. Assume that each computer initially holds at most d elements of A and at most d elements of B, and it will need to hold at most d elements of X. Then all triangles in  can be processed (over fields or semirings) in O(κ + d + log m) rounds.


We emphasize that this strictly improves over <cit.> in two different ways: our result is applicable in a more general setting, and we avoid the factor-2 loss in the exponent of the running time.



 §.§ High-level plan


Our high-level plan for proving <ref> goes as follows:

    
  * We start with an arbitrary collection of κ n triangles  over the original set of 3n nodes V = I ∪ J ∪ K.
    
  * We construct a new virtual collection of κ n triangles ' over a new set of O(n) virtual nodes V'.
    
  * The new virtual collection of triangles is balanced in the sense that each virtual node touches only κ triangles.
    
  * The virtual nodes are assigned to real computers so that each real computer is responsible for only O(1) virtual nodes.
    
  * We show how we can route original data to the virtual nodes, so that each virtual node only needs to do work for each triangle that it touches, and we can also aggregate the results back to the original nodes that need it.

The most challenging part is the final step, routing.



 §.§ Virtual instance


For each node v ∈ V, let (v) be the set of triangles in T ∈ with v ∈ T. Let t(v) = |(v)| and ℓ(v) = ⌈ t(v)/κ⌉. Our set of virtual nodes is

    V'    = I' ∪ J' ∪ K', 
    
        I'    = { (i,x) | i ∈ I, x ∈{1,2,…,ℓ(i)}}, 
    
        J'    = { (k,x) | j ∈ J, x ∈{1,2,…,ℓ(j)}}, 
    
        K'    = { (j,x) | k ∈ K, x ∈{1,2,…,ℓ(k)}}.

The key observation is that the set of virtual nodes is not too large:

    |V'| = ∑_v ℓ(v) ≤∑_v (t(v)/κ + 1) ≤ 3||/κ + n = 4n.

Now we can construct a new balanced collection of triangles ' over V' such that each virtual node v' ∈ V' touches at most κ triangles: for the first κ triangles in (v) we replace v with the virtual node (v,1), for the next κ triangles we replace v with (v,2), etc.

We assign the virtual nodes to real computers so that each computer is responsible for at most 4 virtual nodes. Therefore, in what follows, we can assume that the virtual nodes are computational entities, and we can simulate their work in the real computer network with constant overhead.



 §.§ Routing


Now let us see how we will process the triangles. Consider a triangle {i,j,k}∈, and let its unique virtual copy be {i',j',k'}∈'. Conceptually, we proceed as follows:

    
  * The real computer p(i,j) that initially holds A_ij transmits it to virtual computer i'.
    
  * The real computer p(j,k) that initially holds B_jk transmits it to virtual computer i'.
    
  * Virtual computer i' computes the product A_ij B_jk and transmits it to the real computer p(i,k) that is responsible for storing X_ik.

If we did this in a naive manner, it would be prohibitively expensive; a more careful routing scheme is needed; see <ref> for an illustration.



First consider step (1): the transmission of A_ij from p = p(i,j) to i'. By assumption, i,j is part of at most m triangles, and a single real computer initially holds at most d distinct values of matrix A. On the other hand, a single virtual computer needs to receive up to κ messages, so the total number of messages to transmit is bounded by κ n.

The key idea is that we will do intermediate-point routing. Form an array of all triples (i,j,i') such that A_ij is needed by i', and order the array lexicographically (i,j,i'). For each (i,j) the first triple (i,j,i') is called the anchor element for (i,j). The array has at most κ n triples, and we can assign the triples to our real computers so that each computer has at most κ triples. For each (i,j) we write q(i,j) for the computer that holds the anchor element for (i,j).

Now we first route A_i,j from computer p(i,j) to computer q(i,j). This can be completed in O(d+κ) rounds, as each computer has d outgoing messages and κ incoming messages. We can implement this, for example, by considering the bipartite graph in which on one side we have sender, on one side we have recipients: the sender-side has maximum degree d and the recipient-side has maximum degree κ, and hence we can find a proper edge coloring with O(d+κ) colors and use the color classes to schedule the messages. This way in each round each computer sends and receives at most one message.

Next, we will spread the messages so that each computer that holds a triple (i,j,i') knows the value of A_ij. Thanks to the array being sorted, this is now easy. Let us focus on some pair (i,j). Recall that q(i,j) is the first computer that holds any triple of the form (i,j,·), and it already knows the value of A_ij. Let r(i,j) be the last computer that holds any triple of the form (i,j,·). Now our task is to simply spread A_ij from q(i,j) to all computers in the range q(i,j)+1, …, r(i,j). If r(i,j) = q(i,j), there is nothing to be done. On the other hand, if r(i,j)  q(i,j), we have a convenient situation: computer q(i,j) only needs to spread the value of A_ij, and computers q(i,j)+1, …, r(i,j) only need to receive the value of A_ij. Hence, we can first let q(i,j) inform q(i,j)+1, and then q(i,j)+1, …, r(i,j) form a broadcast tree that distributes A_ij to each of the computers. Note that all these broadcast trees involve disjoint sets of computers, and we can implement broadcast in parallel. We have at most m triples of form (i,j,·), and hence (using a very sloppy estimate) at most m computers involved in each broadcast operation. A broadcast tree of depth O(log m) suffices.

At this point each real computer holds up to κ triples of (i,j,i') and knows the corresponding value A_ij. Then we route A_ij to virtual computer i'; as everyone needs to send and receive up to κ messages, this can be implemented in O(κ) rounds.

Step (2) is essentially identical to step (1); instead of values A_ij and triples (i,j,i'), we transmit values B_jk and use triples (j,k,i').

For step (3) we do the converse of step (1). We now form triples (i,k,i') that indicate that virtual computer i' holds a product A_ij B_jk that needs to be accumulated to X_ik. We sort the array as above, and then first use O(κ) rounds to route the products from virtual computers to real computers. Now each real computer that holds multiple triples of the form (i,k,·) can locally aggregate all these products into a single sum. Then we define q(i,j) and r(i,j) as above, and construct a convergecast tree with which computers q(i,j)+1, …, r(i,j) compute the grand total of the products A_ij B_jk they have received, and finally q(i,j)+1 relays this information to the anchor q(i,j), which can compute the sum X_ij. Finally, we route X_ij to the computer that needs to report it, using additional O(κ + d) rounds.

This completes the proof of <ref>. In what follows, we will present applications of this result.



§ ALGORITHM FOR [US : US : AS]


In our terminology, <cit.> proved the following statement:


    In the supported low-bandwidth model, sparse matrix multiplication of the form [ :  : ] can be computed in O(d^) rounds over semirings and in O(d^) rounds over fields.


In this section, we both improve the running time and widen the scope of applicability:


    In the supported low-bandwidth model, sparse matrix multiplication of the form [ :  : ] can be computed in O(d^) rounds over semirings and in O(d^) rounds over fields.




 §.§ Preliminaries


We start with the following technical lemma:


    Let  be defined by an instance of [ :  : ]. For each node x ∈ V there are at most d^2 triangles T ∈ with x ∈ T.


    Consider first the case of × =. If x ∈ I, then x is incident to at most d nodes j ∈ J and at most d nodes k ∈ K, and hence in total there can be at most d^2 triangles of the form {x,j,k}. If x ∈ J, then x is incident to at most d nodes i ∈ I, and each of them is incident to at most d nodes k ∈ K, and hence in total there can be at most d^2 triangles of the form {i,x,k}. The case of x ∈ K is similar to x ∈ J. For the cases of × = and × =, permute the roles of I, J, and K as appropriate.


We have two simple corollaries that will be useful shortly:


    Let  be defined by an instance of [ :  : ]. For each pair of nodes u,v ∈ V there are at most d^2 triangles T ∈ with {u,v}∈ T.



    The total number of triangles in any [ :  : ] instance is bounded by d^2n.




 §.§ Part 1: handling dense parts


In the first part we follow the basic structure of <cit.>, with two changes:

    
  * we generalize it to [::],
    
  * we choose different parameter values so that we can fully benefit from the new algorithm for part 2.

We will then diverge from <cit.> in part 2, where we will apply <ref>.

Similar to <cit.>, we first show that if there are many triangles, we can find a dense cluster:


    Let  be defined by an instance of [ :  : ], and let ⊆. If || ≥ d^2- n for some > 0, then there exists a cluster U ⊆ V with |[U]| ≥ d^3-4/24.


    This is a generalization of <cit.>, which holds for [ :  : ]. The original proof makes use of the following facts about G():
    
        
  * there are at most dn edges between J and K,
        
  * each i ∈ I is incident to at most d nodes of J,
        
  * each i ∈ I is incident to at most d nodes of K,
        
  * each j ∈ J is incident to at most d nodes of I,
        
  * each k ∈ K is incident to at most d nodes of I.
    
    All of these also hold for × =. For the cases of × = and × =, permute the roles of I, J, and K as appropriate.


Then by plugging in <ref> instead of <cit.> in the proof of <cit.>, we obtain the following corollary:


    Let δ > 0 and > 0, and assume that d is sufficiently large.
    Let  be defined by an instance of [ :  : ], and let ⊆. Assume || ≥ d^2-n. Then we can partition  into disjoint sets  and ' such that  is clustered and || ≥ d^2-5-4δn.


    We follow the strategy in the proof of <cit.>: Apply <ref> to find a cluster U. Put the triangles T ∈[U] that are fully contained in U to , and put the triangles T ∈ that only partially touch U to '. Repeat until there are sufficiently few triangles left; put all remaining triangles to '.

    The original proof makes use of the fact that each node x ∈ V is contained in at most d^2 triangles. This is trivial for [ :  : ]. For [ :  : ] we apply <ref>, and then the original analysis holds verbatim.


Finally, by plugging in <ref> instead of <cit.> in the proof of <cit.>, we obtain the following corollary:


    Let δ > 0 and 0 ≤γ <, and assume that d is sufficiently large.
    Let  be defined by an instance of [ :  : ], and let ⊆.
    Assume || ≤ d^2-γn. Then we can partition  into disjoint sets = _1 ∪…∪_L ∪' such that each _i is clustered, L ≤ 144 d^5 - γ + 4δ and |'| ≤ d^2- n.


    We follow the proof of <cit.>: Apply <ref> repeatedly to construct clusterings _1, _2, …. Stop once the residual set of triangles ' is sufficiently small. The original analysis holds verbatim.


Now we are ready to prove the following result. Note that unlike the analogous result from <cit.>, we have got here a convenient formulation in which the exponent in the time complexity matches the exponent in the size of the residual part:

    Let  be defined by an instance of [ :  : ]. We can partition  into _1 ∪_2 such that:
    
        
  * For matrix multiplication over semirings, _1 can be solved in O(d^) rounds and |_2| ≤ d^ n.
        
  * For matrix multiplication over fields, _1 can be solved in O(d^) rounds and |_2| ≤ d^ n.
    


    For the case of semirings, we apply <ref> with the parameters δ, γ, and  shown in <ref>. In the first step, we start with γ = 0; note that by <ref> we satisfy the assumption that our initial set of triangles has size at most d^2-γn = d^2n. We pick a very small δ and choose a suitable . Then we use <ref> to process each of _i. The total processing time will be O(d^α), where α = 5-γ+4δ+4/3. Our choice of  is optimized so that α≤ and hence we do not exceed our time budget. After processing all of _i, we are left with the residual set of triangles of size |'| ≤ d^β n, where β = 2 -. We then repeat the same process, starting with γ = 2-β, and again optimizing . We continue until β≤.

    For the case of fields, the idea is the same, but we make use of the parameter values in <ref>.








 §.§ Part 2: handling few triangles


Let us recap: We started with an arbitrary instance, with possibly up to d^2 n triangles. Then we have used O(d^α) rounds to process “dense” parts, and we are left with only d^α n triangles; here α = for semirings and α = for fields.

Now we can apply <ref> to handle all these remaining d^α n triangles in O(d^α) rounds. When we apply the lemma, we set κ = d^α, and by <ref> we have m = d^2. The overall running time is O(d^α + d + log d^2) = O(d^α). This completes the proof of <ref>.




§ MORE GENERAL ALGORITHMS


We will now step beyond [::] and consider more general settings. We will show that we can handle much more general notions of sparsity in O(d^2 +log n) rounds. For each case we prove a bound on the total number of triangles, and then apply <ref>.



 §.§ Algorithm for [US : AS : GM]



    Let  be defined by an instance of [::]. Then || ≤ d^2n.


    We give a proof for × =; the remaining cases can be proved similarly.
    We can bound the number of triangles of the form {i,j,k} with i ∈ I, j ∈ J, and k ∈ K as follows: there are at most dn edges of the form {j,k}, since B ∈. For each such edge {j,k} there are at most d edges of the form {i,j}, since A ∈. It follows that the number of triangles is at most d^2 n.



    In the supported low-bandwidth model, sparse matrix multiplication of the form [::] can be computed in O(d^2 + log n) rounds over semirings and fields.


    Follows from <ref> and <ref> by setting κ = d^2 and m = n.




 §.§ Algorithm for [BD : AS : AS]


When we consider [::], it will be convenient to decompose  into  and .


    Let  be defined by an instance of [::]. Then || ≤ d^2n.


    Let us consider the case × =; the other permutations are similar.
    We can bound the number of triangles of the form {i,j,k} with i ∈ I, j ∈ J, and k ∈ K as follows: there are at most dn edges of the form {i,k}, since X ∈. For each such edge {i,k} there are at most d edges of the form {i,j}, since A ∈. It follows that the number of triangles is at most d^2 n.



    Let  be defined by an instance of [::]. Then || ≤ d^2n.


    Transpose the matrices in <ref>.



    Let  be defined by an instance of [::]. Then || ≤ 2d^2n.


    Let us consider the case × =; the other permutations are similar. As we discussed in <ref>, we can decompose X = AB into X = A_1B + A_2B, where A_1 ∈ and A_2 ∈. The claim follows by applying <ref> to X_1 = A_1B and <ref> to X_2 = A_2B.



    In the supported low-bandwidth model, sparse matrix multiplication of the form [::] can be computed in O(d^2 + log n) rounds over semirings and fields.


    Follows from <ref> and <ref> by setting κ = 2d^2 and m = n.




§ CONDITIONAL HARDNESS


In this section we connect the hardness of sparse matrix multiplication with simpler primitives, such as broadcasting.



 §.§ Broadcasting and aggregation

For settings between [::] and [::] our upper bound from <ref> has an additive O(log n) term, which comes from broadcast and convergecast operations. In the following lemma, we show that some broadcasting and/or aggregation is needed in all of these cases. We note that a simple low-bandwidth adaptation of the proof in <cit.> shows an Ω(log n)-round lower bound for computing a sum over values held by all computers, and hence this result also implies an unconditional lower bound.



    In the supported low-bandwidth model, sparse matrix multiplication of the form [::] is at least as hard as computing a sum of n values distributed among n computers, or broadcasting a single value to n computers, even for d = 1.


    Consider first × =. A special case of this is computing AB = X such that all nonzeros in A are in row 1, all nonzeros in B are in column 1, and we are only interested in computing element (1,1) of X. Furthermore, let all nonzeros of B be 1; we have the following task (where “?” indicates elements of the result matrix that we are not interested in):
    
    [ a_1 a_2   … a_n;   0   0   …   0;   4;   0   0   …   0;     ]×[ 1 0 … 0; 1 0 … 0; 4; 1 0 … 0;   ]
            =
            [ x ? … ?; ? ? … ?; 4; ? ? … ?;   ].

    Now this is equivalent to computing a sum x = ∑_j a_j, in a setting in which initially each computer holds one a_j.

    Then consider × =. A special case of this is computing AB = X such that all nonzeros in A are in column 1, the only nonzero of B is at (1,1), and we are only interested in computing elements in the first column of X. Furthermore, let all nonzeros of A be 1; we have the following task:
    
    [ 1 0 … 0; 1 0 … 0; 4; 1 0 … 0;   ]×[ b 0 … 0; 0 0 … 0; 4; 0 0 … 0;   ]
            =
            [ x_1   ?   …   ?; x_2   ?   …   ?;   4; x_4   ?   …   ?;     ].

    Now we need to output x_1 = … = x_n = b, and each computer has to hold one of these values. Hence, this is equivalent to broadcasting the value b to all computers.

    The case of × = is also a broadcast task, by a similar argument.




 §.§ Dense matrix multiplication


Let us now argue why [::] is unlikely to be solvable in O(d^2 + log n) rounds, or even in f(d) + O(n^) rounds for sufficiently small values of .


    Computing [::] for d = 1 in T(n) rounds implies an algorithm that solves dense matrix multiplication in T'(n) = n T(n^2) rounds.


    Assume  solves [::] in T(n) rounds for d = 1. Let m = √(n). We use  as a black box to design an algorithm ' that multiplies dense m × m matrices using m computers in T'(m) rounds. Algorithm ' works as follows: given an input X' = A'B' with dimensions m × m, construct an instance X = AB of dimensions n × n by padding with zeros:
    
    [ A'  0;  0  0;    ]×[ B'  0;  0  0;    ]
            =
            [ X'  ?;  ?  ?;    ].

    This is now an average-sparse instance with m^2 = n nonzeros. Then apply algorithm , so that each of the m computers simulates m = √(n) computers; this way we have n virtual computers available for running . The running time of  is T(n) = T(m^2), and hence the simulation completes in time T'(m) = m T(m^2).


In particular, if we could achieve T(n) = o(n^1/6) for sparse matrix multiplication over semirings, then we could do dense matrix multiplication over semirings in T'(n) = o(n^4/3) rounds, and it would imply major breakthroughs in our understanding of dense matrix multiplication.

Similarly, if we could achieve T(n) = O(n^0.01) over fields, then we could do dense matrix multiplication over fields in T'(n) = O(n^1.02) rounds, which is far from the current upper bound of O(n^).

This is the main justification why our upper bound in <ref> is applicable up to [::], but not all the way to [::].



 §.§ Routing


Next we consider problem settings that are at least as hard as routing many elements to a single computer. We will assume here that the assignment of the input and output values to computers only depends on the structure of the input, and not on the numerical values. We will write A^v, B^v, and X^v to denote the set of elements of matrices A, B, and X that are held by computer v.

We first consider [::]. We prove a hardness result for × =; the case of × = is symmetric, but × = is left for future work.


    In the supported low-bandwidth model, to solve × =, at least one node needs to output Ω(√(n)) values originally held by other nodes.



    For the sake of simplicity, we give a proof when matrix A contains at most 2n nonzero elements, i.e. d=2. Let all the elements of A except a_i,i, a_i,(i  mod n)+1 for all i ∈ [n] are equals to 0; that is, our task is to compute
    
    [ a_1,1 a_1,2     0     0     …     0     0;     0 a_2,2 a_2,3     0     …     0     0;     7; a_n,1     0     0     0     …     0 a_n,n;       ]×
            B
            =
            X,

    where B is a general (dense) matrix, and we are interested in all values of X.
    
    Each computer in the network contains two (nonzero) elements of A, n elements of B and n elements of X. We divide the analysis into the following two cases; consider a computer v:
    
        
  * If X^v contains at least √(n) elements from one of the columns of X: let us assume that X^v contains √(n) elements from column j of X i.e.,  
    x_i_1,j, x_i_2,j, …, x_i_√(n),j∈ X^v.
 We know that 
    x_i,j = a_i,i b_i,j + a_i,(i  mod n)+1 b_(i  mod n)+1,j.
 Let us set b_i,j =1, and a_i,(i  mod n)+1 = 0 for all i,j ∈ [n]. This implies x_i_l,j = a_i_l,i_l for all i_l ∈ [√(n)]. Since computer v initially holds only two of these a_i_l,i_l, we need to route at least Ω(√(n)) values from other computers to computer v.

        
  * If X^v contains less than √(n) elements from every column of X: let 
    I ={(i,j) | x_i,j∈ X^v}.
 If  we set  a_i,i =1 and a_i,((i-1)  mod n)+1 = 0 for all i ∈ [n], we have X= B. Thus, v must output all the elements of b_i,j such that (i,j) ∈ I. Therefore, v must receive all these elements.        
        Now it might happen that these n elements of B are already stored at v. However, we will prove that computer v also requires at least √(n) more elements stored at other computers in order to produce correct output for different values of nonzero elements of A and B. Note that if X^v contains less than √(n) nonzero elements from each column of X, then the number of columns of X from which X^v contains an element is more than √(n). Let C^v be the set of columns of X such that X^v contains some element from these columns. From above, we can say that in each C_i ∈ C^v there exist two elements x_j,i,x_j+1,i such that x_j,i∈ X^v and x_j+1,i∉ X^v (because X^v contains less than √(n) elements from each column). Let 
    J= {j | x_j,i∈ X^v  and  x_j+1,i∉ X^v }.
 Now we take a different assignment where we keep matrix B the same, and in matrix A we set a_j,j = 0 and a_j,(j  mod n +1) = 1 for all j ∈ J. This will result in x_j,i = b_j+1,i for all i and j such that x_j,i∈ X^v and x_j+1,i∉ X^v. We know that there exist more than √(n) such x_j,i. From the above paragraph, we know that v must contain all b_j,i such that x_j,i∈ X^v. Therefore, we can say that v does not contain b_j+1,i such that x_j,i∈ X^v and x_j+1,i∉ X^v. Thus, v needs to receive Ω(√(n)) values from other computers.


Next we consider [::]. We prove a hardness result for × =; the case of × = and × = is left for future work. We consider the special case of × =; the case of × = is at least as hard:


    In the supported low-bandwidth model, to solve × =, at least one node needs to output Ω(√(n)) values originally held by other nodes.



We prove the lower bound for the case when d=1, i.e. each node in the network contains one element of A, one element of B and n elements of X. Let A and B be matrices such that all the elements of these matrices are zero except a_i,1 and b_1,i for all i ∈ [n]:

    [ a_1,1     0     …     0; a_2,1     0     …     0;     4; a_n,1     0     …     0;       ]×[ b_i,1 b_1,2     … b_1,n;     0     0     …     0;     4;     0     0     …     0;       ]
        =
        X.

Similar to the above, we divide the analysis into two cases.

    
  * If X^v contains at least √(n) elements from one of the columns of X: let us assume that X^v contains √(n) elements from column j of X i.e.,  
    x_i_1,j, x_i_2,j, … , x_i_√(n),j∈ X^v.
 Let us set b_1,i = 1 for all i ∈ n. This implies that x_i_k, j = a_i_k,1 for all k ∈√(n). We know that v contains only one element of A, and therefore we need to receive Ω(√(n)) elements from other computers.
    
  * If X^v contains less than √(n) elements from each column of X: in such cases we know that there exist at least √(n) elements x_i_1,j_1, x_i_2,j_2, … x_i_√(n),j_√(n) such that j_i≠ j_k for i ≠ k and i,k ∈ [√(n)]. Now consider an assignment where we set all a_i,1 = 1. This implies x_i_k,j_k = b_1,j_k for all k ∈√(n). Again we can conclude that we need to receive Ω(√(n)) elements from other computers.

 
plainurl



