The I/O complexity of hybrid algorithms for square matrix multiplication

TL;DR

This paper establishes tight lower bounds on the I/O complexity of hybrid matrix multiplication algorithms combining Strassen-like and standard methods, providing new insights into their efficiency limits in hierarchical memory systems.

Contribution

It introduces the first I/O lower bounds for non-uniform, non-stationary hybrid matrix multiplication algorithms, extending the understanding of their computational complexity.

Findings

Derived tight lower bounds for hybrid algorithms in hierarchical memory.

Extended bounds to non-uniform, non-stationary algorithms with recursive variations.

Discussed implications for parallel processing models.

Abstract

Asymptotically tight lower bounds are derived for the I/O complexity of a general class of hybrid algorithms computing the product of $n \times n$ square matrices combining ``\emph{Strassen-like}'' fast matrix multiplication approach with computational complexity $\Theta{n^{\log_2 7}}$, and ``\emph{standard}'' matrix multiplication algorithms with computational complexity $\Omega\left(n^3\right)$. We present a novel and tight $\Omega\left(\left(\frac{n}{\max\{\sqrt{M},n_0\}}\right)^{\log_2 7}\left(\max\{1,\frac{n_0}{M}\}\right)^3M\right)$ lower bound for the I/O complexity a class of ``\emph{uniform, non-stationary}'' hybrid algorithms when executed in a two-level storage hierarchy with $M$ words of fast memory, where $n_0$ denotes the threshold size of sub-problems which are computed using standard algorithms with algebraic complexity $\Omega\left(n^3\right)$. The lower bound is actually derived for the more general class of ``\emph{non-uniform, non-stationary}'' hybrid algorithms which allow recursive calls to have a different structure, even when they refer to the multiplication of matrices of the same size and in the same recursive level, although the quantitative expressions become more involved. Our results are the first I/O lower bounds for these classes of hybrid algorithms. All presented lower bounds apply even if the recomputation of partial results is allowed and are asymptotically tight. The proof technique combines the analysis of the Grigoriev's flow of the matrix multiplication function, combinatorial properties of the encoding functions used by fast Strassen-like algorithms, and an application of the Loomis-Whitney geometric theorem for the analysis of standard matrix multiplication algorithms. Extensions of the lower bounds for a parallel model with $P$ processors are also discussed.

Equations

$$|D|\ge \max \{|\mathcal{Y}{}_i^{\prime }\cap A_i|,|\mathcal{Y}{}_i^{\prime }\cap B_i|\}.$$

$$I{O}_{\mathcal{A}}\left(n,M,B\right)\ge \max \{2n^2,c|\mathcal{T}|M^{-1/2},{\nu }_{2}M\}{B}^{-1}$$

$$I{O}_{\mathcal{A}}\left(n,M,B_m,P\right)\ge \max \{c|\mathcal{T}|M^{-1/2},{\nu }_{2}M\}{\left(P{B}_m\right)}^{-1}$$

$$|{\mathcal{T}}_i^{(j)}|\le \sqrt{|\mathcal{Y}{}_{i,\mathbf{A}}^{\prime }||\mathcal{Y}{}_{i,\mathbf{B}}^{\prime }||{\mathcal{Z}}_i^{\prime }|}.$$

$$|D{}^{\prime }|=\sum _{i=1}^{{\nu }_1}|D{}_i^{\prime }|=\sum _{i=1}^{{\nu }_1}\max \{0,|{\mathcal{Z}}_i^{\prime }|-|{\mathcal{T}}_i^{(j)}|/n_i\}\ge |\mathcal{Z}|-|{\mathcal{T}}^{(j)}|/2\sqrt{M},$$

$$|Y|=\sum _{i=1}^{{\nu }_1}|\mathcal{Y}{}_i^{\prime }|\ge \sum _{i=1}^{{\nu }_1}\frac{|{\mathcal{T}}_i^{(j)}|}{\sqrt{|{\mathcal{Z}}_i^{\prime }|}}.$$

$$|D{}_{Y}|\geq\min\bigl{\{}2M,\frac{\sum_{i=1}^{\nu_{1}}|\mathcal{T}^{(j)}_{i}|}{\sum_{i=1}^{\nu_{1}}\sqrt{|\mathcal{Z}_{i}^{\prime}|}}\bigr{\}}=\min\bigl{\{}2M,\frac{|\mathcal{T}^{(j)}|}{\sqrt{|\mathcal{Z}^{\prime}|}}\bigr{\}}.$$

$$|D{}^{\prime\prime}|\geq\min\bigl{\{}2M,|\mathcal{T}^{(j)}|/\sqrt{|\mathcal{Z}^{\prime}|}\bigr{\}}.$$

$$|D{}|\geq|D{}^{\prime}|+|D{}^{\prime\prime}|\geq|\mathcal{Z}^{\prime}|-|\mathcal{T}^{(j)}|/2\sqrt{M}+\min\bigl{\{}2M,|\mathcal{T}^{(j)}|/\sqrt{|\mathcal{Z}^{\prime}|}\bigr{\}}.$$

$$|D|$$

$$IO_{\mathcal{A}}\left(n,M,B\right)\geq\max\{2n^{2},\left(\frac{n}{\max\{n_{0}{},2\sqrt{M}\}}\right)^{\log_{2}7}\left(\max\bigl{\{}1,\frac{n_{0}{}}{2\sqrt{M}}\bigr{\}}\right)^{3}M\}B^{-1}$$

$$IO_{\mathcal{A}}\left(n,M,B_{m}{},P\right)\geq\left(\frac{n}{\max\{n_{0}{},2\sqrt{M}\}}\right)^{\log_{2}7}\left(\max\Bigl{\{}1,\frac{n_{0}{}}{2\sqrt{M}}\Bigr{\}}\right)^{3}\frac{M}{P{}B_{m}{}}.$$

$${\nu }_1=7^i=7^{{\log }_{2}n/n_0}={\left(\frac{n}{n_0}\right)}^{{\log }_{2}7}$$

$$|\mathcal{T}|={\nu }_1{n}_0{}^3={\left(\frac{n}{n_0}\right)}^{{\log }_{2}7}{n}_0{}^3.$$

$${\nu }_2=7^i=7^{{\log }_{2}n/2\sqrt{M}}={\left(\frac{n}{2\sqrt{M}}\right)}^{{\log }_{2}7}$$

$$\frac{a_1}{\sqrt{b_1}}\ge \frac{a_2}{\sqrt{b_2}}\ge ....\ge \frac{a_i}{\sqrt{b_i}}$$

$$\frac{a_1}{\sqrt{b_1}}+\frac{1}{2}\sum _{j=2}^i\frac{a_j}{\sqrt{b_j}}\ge \frac{{\sum }_{j=1}^i{a}_i}{\sqrt{{\sum }_{j=1}^j{b}_j}}$$

$$\frac{a_1}{\sqrt{b_1}}\ge \frac{a_2}{\sqrt{b_2}}\ge ....\ge \frac{a_i+1}{\sqrt{b_i+1}}$$

$$\frac{a_1}{\sqrt{b_1}}+\frac{1}{2}\sum _{j=2}^{i+1}\frac{a_j}{\sqrt{b_j}}\ge \frac{{\sum }_{j=1}^i{a}_j}{\sqrt{{\sum }_{j=1}^i{b}_j}}+\frac{1}{2}\frac{a_{i+1}}{\sqrt{b_{i+1}}}.$$

$$\frac{a^{\prime }}{\sqrt{b^{\prime }}}+\frac{1}{2}\frac{a_{i+1}}{\sqrt{b_{i+1}}}\ge \frac{a^{\prime }+a_{i+1}}{\sqrt{b^{\prime }+b_{i+1}}}=\frac{a^{\prime }}{\sqrt{b^{\prime }+b_{i+1}}}+\frac{a_{i+1}}{\sqrt{b^{\prime }+b_{i+1}}}$$

$$\frac{a^{\prime }}{\sqrt{b^{\prime }}}+\frac{1}{2}\frac{a_{i+1}}{\sqrt{b_{i+1}}}$$

$$=\frac{a}{\sqrt{b}}\left(\frac{x}{\sqrt{y}}+\frac{1-x}{2\sqrt{1-y}}\right)$$

$$\frac{x}{\sqrt{y}}+\frac{1-x}{\sqrt{1-y}}=\frac{2x\sqrt{1-y}+(1-x)\sqrt{y}}{2\sqrt{y}\sqrt{1-y}}\ge 1.$$

$$2x\sqrt{1-y}+\left(1-x\right)\sqrt{y}\ge 2\sqrt{y}\sqrt{1-y},$$

$$2\sqrt{1-y}\left(x-\sqrt{y}\right)\ge \sqrt{y}\left(x-1\right),$$

$$|K^{(1)}|+\frac{1}{2}\sum _{j=2}^7|K^{(j)}|$$

$$\ge \frac{{\sum }_{j=1}^7{a}_j}{\sqrt{{\sum }_{j=1}^7{b}_j}},$$

$${\left(4\sqrt{M}\left(\sqrt{{\delta }^{(1)}}+\frac{1}{2}\sum _{j=2}^7\sqrt{{\delta }^{(j)}}\right)\right)}^2=16M\left({\delta }^{(1)}+\sqrt{{\delta }^{(1)}}\sum _{j=2}^7\sqrt{{\delta }^{(j)}}+{\left(\frac{1}{2}\sum _{j=2}^7\sqrt{{\delta }^{(j)}}\right)}^2\right).$$

$${\left(4\sqrt{M}\left(\sqrt{{\delta }^{(1)}}+\frac{1}{2}\sum _{j=2}^7\sqrt{{\delta }^{(j)}}\right)\right)}^2\ge 16M\sum _{j=1}^7{\delta }^{(j)}\ge {\left(4\sqrt{M\left(|Z|-2|Q|\right)}\right)}^2.$$

$${\left(4\sqrt{M\left(|Z|-2\left(|D|-|D^{\prime }|\right)\right)}\right)}^2$$

$$=16M\left(|Z|-2|D|\right)+32M|D^{\prime }|.$$

$${\left(\varphi +|D^{\prime }|\right)}^2={\left(4\sqrt{M\left(|Z|-2\left(|D|-|D^{\prime }|\right)\right)}\right)}^2>|D^{\prime }{|}^2.$$

Full text

The I/O complexity of hybrid algorithms

for square matrix multiplication

Lorenzo De Stefani lorenzo@cs.brown.edu Department of Computer Science, Brown University

Asymptotically tight lower bounds are derived for the I/O complexity of a general class of hybrid algorithms computing the product of $n\times n$ square matrices combining “Strassen-like” fast matrix multiplication approach with computational complexity $\Theta\left(n^{\log_{2}7}\right)$, and “standard” matrix multiplication algorithms with computational complexity $\Omega\left(n^{3}\right)$. We present a novel and tight $\Omega\left(\left(\frac{n}{\max\{\sqrt{M},n_{0}{}\}}\right)^{\log_{2}7}\left(\max\{1,\frac{n_{0}{}}{M}\}\right)^{3}M\right)$ lower bound for the I/O complexity a class of “uniform, non-stationary” hybrid algorithms when executed in a two-level storage hierarchy with $M$ words of fast memory, where $n_{0}{}$ denotes the threshold size of sub-problems which are computed using standard algorithms with algebraic complexity $\Omega\left(n^{3}\right)$.

The lower bound is actually derived for the more general class of “non-uniform, non-stationary” hybrid algorithms which allow recursive calls to have a different structure, even when they refer to the multiplication of matrices of the same size and in the same recursive level, although the quantitative expressions become more involved. Our results are the first I/O lower bounds for these classes of hybrid algorithms. All presented lower bounds apply even if the recomputation of partial results is allowed and are asymptotically tight.

The proof technique combines the analysis of the Grigoriev’s flow of the matrix multiplication function, combinatorial properties of the encoding functions used by fast Strassen-like algorithms, and an application of the Loomis-Whitney geometric theorem for the analysis of standard matrix multiplication algorithms. Extensions of the lower bounds for a parallel model with $P{}$ processors are also discussed.

1 Introduction

Data movement plays a critical role in the performance of computing systems, in terms of both time and energy. This technological trend [28] appears destined to continue, as physical limitations on minimum device size and on maximum message speed lead to inherent costs when moving data, whether across the levels of a hierarchical memory system or between processing elements of a parallel system [10]. While the communication requirements of algorithms have been widely investigated in literature, obtaining significant and tight lower bounds based on such requirements remains an important and challenging task.

In this paper, we focus on the I/O complexity of a general class of hybrid algorithms for the computing the product of square matrices which combine. Such algorithms combine fast algorithms with base case $2\times 2$ similar to Strassen’s matrix multiplication algorithm [35] with algebraic (or computational) complexity $\mathcal{O}\left(n^{\log_{2}7}\right)$ with standard (or classic) matrix multiplication algorithms with algebraic complexity $\Omega\left(n^{3}\right)$. Further, these algorithms allow recursive calls to have a different structure, even when they refer to the multiplication of matrices in the same recursive level and of the same input size. These algorithms are referred in literature as “non-uniform, non-stationary”. This class includes, for example, algorithms that optimize for input sizes [15, 16, 20]. Matrix multiplication is a pervasive primitive utilized in many applications.

While of actual practical importance, to the best of our knowledge, no characterization of the I/O complexity of such algorithms has presented before this work. This is likely due to the the fact that the irregular nature of hybrid algorithms and, hence, the irregular structure of the corresponding Computational Directed Acyclic Graphs (CDAGs), complicates the analysis of the combinatorial properties of the CDAG which is the foundation of many of I/O lower bound technique presented in literature (e.g., [8, 18, 29]).

The technique used in this work overcomes such challenges and yields asymptotically tight I/O lower bounds which hold even if recomputation of intermediate values is allowed.

Previous and Related Work:

Strassen [35] showed that two $n\times n$ matrices can be multiplied with $O(n^{\omega})$ operations, where $\omega=\log_{2}7\approx 2.8074$, hence with asymptotically fewer than the $n^{3}$ arithmetic operations required by the straightforward implementation of the definition of matrix multiplication. This result has motivated a number of efforts which have lead to increasingly faster algorithms, at least asymptotically, with the current record being at $\omega<2.3728639$ [24].

I/O complexity has been introduced in the seminal work by Hong and Kung [18]; it is essentially the number of data transfers between the two levels of a memory hierarchy with a fast memory of $M$ words and a slow memory with an unbounded number of words. Hong and Kung presented techniques to develop lower bounds to the I/O complexity of computations modeled by computational directed acyclic graphs (CDAGs). The resulting lower bounds apply to all the schedules of the given CDAG, including those with recomputation, that is, where some vertices of the CDAG are evaluated multiple times. Among other results, they established a $\Omega\left(n^{3}/\sqrt{M}\right)$ lower bound to the I/O complexity of standard, definition-based matrix multiplication algorithms, which matched a known upper bound [13]. The techniques of [18] have also been extended to obtain tight communication bounds for the definition-based matrix multiplication in some parallel settings [4, 21, 33] and for the special case of “sparse matrix multiplication” [27]. Ballard et al. generalized the results on matrix multiplication of Hong and Kung [18] in [7, 6] by using the approach proposed in [21] based on the Loomis-Whitney geometric theorem [25, 36].

In an important contribution, Ballard et al. [8], obtained an $\Omega((n/\sqrt{M})^{\log_{2}7}M)$ I/O lower bound for Strassen’s algorithm, using the “edge expansion approach”. The authors extend their technique to a class of “Strassen-like” fast multiplication algorithms and to fast recursive multiplication algorithms for rectangular matrices [5]. This result was later generalized to increasingly broader classes of “Strassen-like” algorithms by Scott et. al [31] using the “path routing” technique, and De Stefani [14] using a combination the concept of Grigoriev’s flow of a function and the “dichotomy width” technique [9]. While the previously mentioned results hold only under the restrictive assumption that no intermediate result may be more than once (i.e., the no-recomputation assumption), in [11] Bilardi and De Stefani introduced the first asymptotically tight I/O lower bound which holds if recomputation is allowed. Their technique was later extended to the analysis of Strassen-like algorithms with base case $2\times 2$ [26], and to the analysis of Toom-Cook integer multiplication algorithms [12].

A parallel, “communication avoiding” implementation of Strassen’s algorithm whose performance matches the known lower bound [8, 31], was proposed by Ballard et al. [3]. A communication efficient algorithm for the special case of sparse matrices based on Strassen’s algorithm was presented in [22].

In [32], Scott derived a lower bound for the I/O complexity of a class of uniform, non-stationary algorithms combining Strassen-like algorithm with recursive standard algorithms. This result holds only under the restrictive no-recomputation assumption.

To the best of our knowledge, ours is the first work presenting asymptotically tight I/O lower bounds for non-uniform, non-stationary hybrid algorithms for matrix multiplication that hold when recomputation is allowed.

Our results:

We present the first I/O lower bound for a class $\mathfrak{H}{}$ of non-uniform, non-stationary hybrid matrix multiplication algorithms when executed in a two-level storage hierarchy with $M$ words of fast memory. Algorithms in $\mathfrak{H}{}$ combine fast Strassen-like algorithms with base case $2\times 2$ with algebraic complexity $\Theta\left(n^{\log_{2}7}\right)$, and standard algorithms based on the definition with algebraic complexity $\Omega\left(n^{3}\right)$. These algorithms allow recursive calls to have a different structure, even when they refer to the multiplication of matrices in the same recursive level and of the same input size. The result in Theorem 7 relates the I/O complexity of algorithms in $\mathfrak{H}{}$ to the number and the input size of an opportunely selected set of the sub-problems generated by the algorithms themselves.

We also present, in Theorem 9, a novel $\Omega\left(\left(\frac{n}{\max\{\sqrt{M},n_{0}{}\}}\right)^{\log_{2}7}\left(\max\{1,\frac{n_{0}{}}{M}\}\right)^{3}M\right)$ lower bound for the I/O complexity of algorithms in a subclass $\mathfrak{UH}\left(n_{0}{}\right){}$ of $\mathfrak{H}{}$ composed by uniform non-stationary hybrid algorithms where $n_{0}{}$ denotes the threshold size of sub-problems which are computed using standard algorithms with algebraic complexity $\Omega\left(n^{3}\right)$.

The previous result by Scott [32] covers only a sub-class of $\mathfrak{UH}\left(n_{0}{}\right){}$ composed by uniform, non-stationary algorithms combining Strassen-like algorithms with the recursive standard algorithm, and holds only assuming that no intermediate value is recomputed. Instead, all our bounds allow for recomputation of intermediate values and are asymptotically tight. As the matching upper bounds do not recompute any intermediate value, we conclude that using recomputation may reduce the I/O complexity of the considered classes of hybrid algorithms by at most a constant factor.

Our proof technique is of independent interest since it exploits to a significant extent the “divide and conquer” nature exhibited by many algorithms. Our approach combines elements from the “G-flow” I/O lower bound technique originally introduced by Bilardi and De Stefani, with an application of the Loomis-Whitney geometric theorem, which has been used by Irony et al. to study the I/O complexity of standard matrix multiplication algorithms [21], to recover an information which relates to the concept of Minimum set introduced in Hong and Kung’s method. We follow the dominator set approach pioneered by Hong and Kung in [18]. However, we focus the dominator analysis only on a select set of target vertices, which, depending on the algorithm structure, correspond either to the outputs of the sub-CDAGs that correspond to sub-problems of a suitable size (i.e., chosen as a function of the fast memory capacity $M$) computed using a fast Strassen-like algorithm, or to the the vertices corresponding to the elementary products evaluated by a standard (definition) matrix multiplication algorithm.

We derive our main results for the hierarchical memory model (or external memory model). Our results generalize to parallel models with P processors. For these parallel models, we derive lower bounds for the “bandwith cost”, that is for the number of messages (and, hence, the number of memory) that must be either sent or received by at least one processor during the CDAG evaluation.

Paper organization:

In Section 2 we outline the notation and the computational models used in the rest of the presentation. In Section 3 we rigorously define the class of hybrid matrix multiplication algorithms $\mathfrak{H}{}$ being considered. In Section 4 we discuss the construction and important properties of the CDAGs corresponding to algorithms in $\mathfrak{H}{}$. In Section 5 we introduce the concept of Maximal Sup-Problem (MSP) and describe their properties which lead to the I/O lower bounds for algorithms in $\mathfrak{H}{}$ in Section 6.

2 Preliminaries

We consider algorithms that compute the product of two square matrices $\mathbf{A}\times\mathbf{B}=\mathbf{C}$ with entries from a ring $\mathcal{R}{}$. We use $A$ to denote the set variables each corresponding to an entry of matrix $\mathbf{A}$. We refer to the number of entries of a matrix $\mathbf{A}$ as its “size” and we denote it as $|A|$. We denote the entry on the $i$-th row of the $j$-th column of matrix $\mathbf{A}$ as $\mathbf{A}[i][j]$.

In this work we focus on algorithms whose execution can be modeled as a computational directed acyclic graph (CDAG) $G=(V,E)$. Each vertex $v\in V$ represents either an input value or the result of a unit-time operation (i.e., an intermediate result or one of the output values) which is stored using a single memory word. For example, each of the input (resp., output) vertices of $G$ corresponds to one of the $2n^{2}$ entries of the factor matrices $\mathbf{A}$ and $\mathbf{B}$ (resp., to the $n^{2}$ entries of the product matrix $\mathbf{C}$). The directed edges in $E$ represent data dependencies. That is, a pair of vertices $u,v\in V$ are connected by an edge $(u,v)$ directed from $u$ to $v$ if and only if the value corresponding to $u$ is an operand of the unit time operation which computes the value corresponding to $v$. A directed path connecting vertices $u,v\in V$ is an ordered sequence of vertices starting with $u$ and ending with $v$, such that there is an edge in $E$ directed from each vertex in the sequence to its successor.

We say that $G^{\prime}=(V^{\prime},E^{\prime})$ is a sub-CDAG of $G=(V,E)$ if $V^{\prime}\subseteq V$ and $E^{\prime}\subseteq E\cap(V^{\prime}\times V^{\prime})$. Note that, according to this definition, every CDAG is a sub-CDAG of itself. We say that two sub-CDAGs $G^{\prime}=(V^{\prime},E^{\prime})$ and $G^{\prime\prime}=(V^{\prime\prime},E^{\prime\prime})$ of $G$ are vertex disjoint if $V^{\prime}\cap V^{\prime\prime}=\emptyset$. Analogously, two directed paths in $G$ are vertex disjoint if they do not share any vertex.

When analyzing the properties of CDAGs we make use of the concept of dominator set originally introduced in [18]. We use the following – slightly different – definition:

Definition 1** (Dominator set).**

Given a CDAG $G=(V,E)$, let $I\subset V$ denote the set of its input vertices. A set $D\subseteq V$ is a dominator set for $V^{\prime}\subseteq V$ with respect to $I^{\prime}\subseteq I$ if every path from a vertex in $I^{\prime}$ to a vertex in $V^{\prime}$ contains at least a vertex of $D$. When $I^{\prime}=I$, $D$ is simply referred as “a dominator set for $V^{\prime}$”.

I/O model and machine models:

We assume that sequential computations are executed on a system with a two-level memory hierarchy, consisting of a fast memory or cache of size $M$ (measured in memory words) and a slow memory of unlimited size. An operation can be executed only if all its operands are in cache. We assume that each entry of the input and intermediate results matrices (including entries of the output matrix) is maintained in a single memory word (the results trivially generalize if multiple memory words are used).

Data can be moved from the slow memory to the cache by read operations and, in the other direction, by write operations. These operations are also called I/O operations. We assume the input data to be stored in slow memory at the beginning of the computation. The evaluation of a CDAG in this model can be analyzed by means of the “red-blue pebble game” [18]. The number of I/O operations executed when evaluating a CDAG depends on the “computational schedule”, that is, it depends on the order in which vertices are evaluated and on which values are kept in/discarded from cache.

The I/O complexity $IO_{G}(M)$ of a CDAG $G$ is defined as the minimum number of I/O operations over all possible computational schedules. We further consider a generalization of this model known as the “External Memory Model” by Aggarwal and Vitter [2], where $B\geq 1$ values can be moved between cache and consecutive slow memory locations with a single I/O operation. For $B=1$, this model clearly reduces to the red-blue pebble game.

Given an algorithm $\mathcal{A}$, we only consider “parsimonious execution schedules”, that is schedules such that: (i) each time an intermediate result (excluding the output entries of $\mathbf{C}$) is computed, such value is then used to computed to compute at least one of the values of which it is an operand before being removed from the memory (either the cache or slow memory); and (ii) any time an intermediate result is read from slow to cache memory, such value is then used to computed to compute at least one of the values of which it is an operand before being removed from the memory or moved back to slow memory using a write I/O operation. Clearly, any non-parsimonious schedule $\mathcal{C}$ can be reduced to a parsimonious schedule $\mathcal{C}^{\prime}$ by removing all the steps which violate the definition of parsimonious computation. $\mathcal{C}^{\prime}$ has therefore less computational and I/O operations than $\mathcal{C}$. Hence, restricting the analysis to parsimonious computations leads to no loss of generality.

We also consider a parallel model where $P$ processors, each with a local memory of size $2n^{2}/P\leq M<n^{2}$, are connected by a network. We do not, however, make any assumption on the initial distribution of the input data nor regarding the balance of the computational load among the $P{}$ processors. Processors can exchange point-to-point messages, with every message containing up to $B_{m}{}$ memory words. For this parallel model, we derive lower bounds for the number of messages that must be either sent or received by at least one processor during the CDAG evaluation. The notion of “parsimonious execution schedules” straightforwardly extends to this parallel model.

3 Hybrid matrix multiplication algorithms

In this work, we consider a family of hybrid matrix multiplication algorithms obtained by hybridizing the two following classes of algorithms:

Standard matrix multiplication algorithms: This class includes all the square matrix multiplication algorithms which, given the input factor matrices $\mathbf{A},\mathbf{B}\in\mathcal{R}{}^{n\times n}$, satisfy the following properties:

• •

The $n^{3}$ elementary products $\mathbf{A}[i][j]\mathbf{B}[j][i]$, for $i,j=0,\ldots,n-1$, are directly computed;

• •

Each of the $\mathbf{C}[i][j]$ is computed by summing the values of the $n$ elementary products $\mathbf{A}[i][z]\mathbf{B}[z][j]$, for $z=0,\ldots,n-1$, through a summation tree by additions and subtractions only;

• •

The evaluations of the $\mathbf{C}[i][j]$’s are independent of each other. That is, internal vertex sets of the summation trees of all the $\mathbf{C}[i][j]$’s are disjoint from each other.

This class, also referred in literature as classic, naive or conventional algorithms, correspond to that studied for the by Hong and Kung [18] (for the sequential setting) and by Irony et al. [21] (for the parallel setting). Algorithms in this class have computational complexity $\Omega\left({n^{3}}\right)$. This class includes, among others, the sequential iterative definition algorithm, the sequential recursive divide and conquer algorithm based on block partitioning, and parallel algorithms such as Cannon’s “2D” algorithm [13], the “2.5D” algorithm by Solomonik and Demmel [34], and “3D” algorithms [1, 23].

Fast Strassen-like matrix multiplication algorithms with base case $2\times 2$: This class includes algorithms following a structure similar to that of Strassen’s [35] (see Appendix B) and Winograd’s variation [37] (which reduces the leading coefficient of the arithmetic complexity reduced from 7 to 6). Algorithms in this class generally follow three steps:

Encoding: Generate the inputs, of size $n/2\times n/2$ of seven sub-problems, as linear sums of the input matrices; 2. 2.

Recursive multiplications: Compute (recursively) the seven generated matrix multiplication sub-products; 3. 3.

Decoding: Computing the entries of the product matrix $\mathbf{C}$ via linear combinations of the output of the seven sub-problems.

Algorithms in this class have algebraic complexity $\mathcal{O}\left(n^{\log_{2}7}\right)$, which is optimal for algorithms with base case $2\times 2$ [37].

Remarkably, the only properties of relevance for the analysis of the I/O complexity of algorithms in these classes are those used in the characterization of the classes themselves.

In this work we consider a general class of non-uniform, non-stationary hybrid square matrix multiplication algorithms, which allow mixing of schemes from the fast Strassen-like class with algorithms from the standard class. Given an algorithm $\mathcal{A}$ let $P$ denote the problem corresponding to the computation of the product of the input matrices $\mathcal{A}$ and $\mathcal{B}$. Consider an “instruction function” $f_{\mathcal{A}}(P)$, which, given as input $P$ returns either (a) indication regarding the algorithm from the standard class which is to be used to compute $P$, or (b) indication regarding the fast Strassen-like algorithm to be used to recursively generate seven sub-problems $P_{1},P_{2},\ldots,P_{7}$ and the instruction functions $f_{\mathcal{A}}(P_{i})$ for each of the seven sub-problems. We refer to the class of non-uniform, non-stationary algorithms which can be characterized by means of such instruction functions as $\mathfrak{H}{}$. Algorithms in $\mathfrak{H}{}$ allow recursive calls to have a different structure, even when they refer to the multiplication of matrices in the same recursive level. E.g., some of the sub-problems with the same size may be computed using algorithms form the standard class while others may be computed using recursive algorithms from the fast class. This class includes, for example, algorithms that optimize for input sizes, (for sizes that are not an integer power of a constant integer).

We also consider a sub-class $\mathfrak{UH}\left(n_{0}{}\right){}$ of $\mathfrak{H}{}$ constituted by uniform, non-stationary hybrid algorithms which allow mixing of schemes from the fast Strassen-like class for the initial $\ell$ recursion levels, and then cut the recursion off once the size the generated sub-problems is smaller or equal to a set threshold $n_{0}{}\times n_{0}{}$, and switch to using algorithm form the standard class. Algorithms in this class are uniform, i.e., sub-problems of the same size are all either recursively computed using a scheme form the fast class, or are all computed using algorithms from the standard class.

This corresponds to actual practical scenarios, as the use of Strassen-like algorithms is mostly beneficial for large input size. As the size of the input of the recursively generated sub-problems decreases, the asymptotic benefit of fast algorithms is lost due to the increasing relative impact of the constant multiplicative factor, and algorithms in the standard class exhibit lower actual algebraic complexity. For a discussion on such hybrid algorithms and their implementation issues we refer the reader to [16, 20] (sequential model) and [15] (parallel model).

4 The CDAG of algorithms in $\mathfrak{H}{}$

Let $G^{\mathcal{A}}=(V_{\mathcal{A}},E_{\mathcal{A}})$ denote the CDAG that corresponds to an algorithm $\mathcal{A}\in\mathfrak{H}{}$ used to multiply input matrices $\mathbf{A},\mathbf{B}\in\mathcal{R}{}^{n\times n}$. The challenge in the characterization of $G^{\mathcal{A}}$ comes from the fact that rather than considering to a single algorithm, we want to characterize the CDAG corresponding to the class $\mathfrak{H}{}$. Further, the class $\mathfrak{H}{}$ is composed by a rich variety of vastly different and irregular algorithm. Despite such variety, we show a general template for the construction of $G^{\mathcal{A}}$ and we identify some of it properties which crucially hold regardless of the implementation details of $\mathcal{A}$ and, hence, of $G^{\mathcal{A}}$.

Construction:

$G^{\mathcal{A}}$ can be obtained by using a recursive construction that mirrors the recursive structure of the algorithm itself. Let $P$ denote the entire matrix multiplication problem computed by $\mathcal{A}$. Consider the case for which, according to the instruction function $f_{\mathcal{A}}(P)$, $P$ is to be computed using an algorithm from the standard class. As we do not fix a specific algorithm, we do not correspondingly have a fixed CDAG. The only feature of interest for the analysis is that, in this case, the CDAG $G^{\mathcal{A}}$ corresponds to the execution of an algorithm from the standard class for input matrices of size $n\times n$.

Consider instead the case for which, according to $f_{\mathcal{A}}(P)$, $P$ is to be computed using an algorithm from the fast class. In the base case for $n=2$ the problem $P$ is computed without generating any further sub-problems. As an example, we present in Figure 1(a) the base case for Strassen’s original algorithm [35]. If $n>2$, then $f_{\mathcal{A}}(P)$ specifies the divide and conquer scheme to be used to generate the seven sub-problems $P_{1},P_{2},\ldots,P_{7}$, and the instruction function for each of them. The sub-CDAGs of $G^{\mathcal{A}}$ corresponding to each of the seven sub-problems $P_{i}$, denoted as $G^{\mathcal{A}}_{P_{i}}$ are constructed according to $f_{\mathcal{A}}(P_{i})$, following recursively the steps discussed previously. $G^{\mathcal{A}}$ can then be constructed by composing the seven sub-CDAGs $G^{\mathcal{A}}_{P_{i}}$. $n^{2}$ disjoint copies of an encoder sub-CDAG $Enc_{A}$ (resp., $Enc_{B}$) are used to connect the input vertices of $G^{2n\times 2n}$, which correspond to the values of the input matrix $\mathbf{A}$ (resp., $\mathbf{B}$) to the appropriate input vertices of the seven sub-CDAGs $G^{\mathcal{A}}_{P_{i}}$; the output vertices of the sub-CDAGs $G^{\mathcal{A}}_{P_{i}}$ (which correspond to the outputs of the seven sub-products) are connected to the appropriate output vertices of the entire $G^{\mathcal{A}}$ CDAG using $n^{2}$ copies of the decoder sub-CDAG $Dec$. We present an example of such recursive construction in Figure 1(b).

Properties of $G^{\mathcal{A}}$:

While the actual internal structure $G^{\mathcal{A}}$, and, in particular, the structure of encoder and decoder sub-CDAGs depends on the specific Strassen-like algorithm being used by $\mathcal{A}$, all versions share some properties of great importance. Let $G(X,Y,E)$ denote an encoder CDAG for a fast multiplication algorithm $2\times 2$ base case, with $X$ (resp., $Y$) denoting the set of input (resp., output) vertices, and $E$ denoting the set of edges directed from $X$ to $Y$.

Lemma 1** (Lemma 3.3 [26]).**

Let $G=(X,Y,E)$ denote an encoder graph for a fast matrix multiplication algorithm with base case $2\times 2$. There are no two vertices in Y with identical neighbors sets.

While the correctness of this Lemma can be simply verified by inspection in the case of Strassen’s algorithm [35], Lemma 1 generalizes the statement to all encoders corresponding to fast matrix multiplication algorithms with base case $2\times 2$. From Lemma 1 we have:

Lemma 2**.**

Let $\mathcal{A}\in\mathfrak{H}{}$ and let $P_{1}$ and $P_{2}$ be any two sub-problems generated by $\mathcal{A}$ with input size greater than $n_{0}{}\times n_{0}{}$, such that $P_{2}$ is not recursively generated while computing $P_{1}$ and vice versa. Then, the sub-CDAGs of $G^{\mathcal{A}}$ corresponding, respectively, to $P_{1}$ and to $P_{2}$ are vertex disjoint.

The following lemma, originally introduced for Strassen’s algorithm in [11] and then generalized for Strassen-like algorithms with base case $2\times 2$ in in [26], captures a connectivity property of encoder sub-CDAGs.

Lemma 3** (Lemma 3.1 [26]).**

Given an encoder CDAG for any Strassen-like algorithm with base case $2\times 2$, for any subset $Y$ of its output vertices, there exists a subset $X$ of its input vertices, with $\min\{|Y|,1+\lceil\left(|Y|-1\right)/2\rceil\}\leq|X|\leq|Y|$, such that there exist $|X|$ vertex-disjoint paths connecting the vertices in $X$ to vertices in $Y$.

The proofs of Lemma 1 and Lemma 3 are based on an argument originally presented by Hopcroft and Kerr [19]. We refer the reader to [26] for the proofs.

5 Maximal sub-problems and their properties

For an algorithm $\mathcal{A}\in\mathfrak{H}{}$, let $P^{\prime}$ denote a sub-problem generated by $\mathcal{A}$. In our presentation we consider the entire matrix multiplication problem an improper sup-problem generated by $\mathcal{A}$. We refer as the ancestor sub-problems of $P^{\prime}$ as the sequence of sub-problems $P^{\prime}_{0},P^{\prime}_{2},\ldots,P^{\prime}_{i}$ generated by $\mathcal{A}$ such that $P^{\prime}_{j+1}$ was recursively generated to compute $P^{\prime}_{j}$ for $j=0,1,\ldots,i-1$, and such that $P$ was recursively generated to compute $P^{\prime}_{i}$. Clearly, if $P^{\prime}$ is the entire problem, then $P^{\prime}$ has no ancestors.

Towards studying the I/O complexity of algorithms in $\mathfrak{H}{}$ we focus on the analysis of a particular set of sub-problems.

Definition 2** (Maximal Sub-Problems (MSP)).**

Let $\mathcal{A}\in\mathfrak{H}{}$ be an algorithm used to multiply matrices $\mathbf{A},\mathbf{B}\in\mathcal{R}{}^{n\times n}$. If $n\leq 2\sqrt{M}$ we say that $\mathcal{A}$ does not generate any Maximal Sub-Problem (MSP).

Let $P_{i}$ be a sub-problem generated by $\mathcal{A}$ with input size $n_{i}\times n_{i}$, with $n_{i}\geq 2M$, and such that all its ancestors sub-problems are computed, according to $\mathcal{A}$ using algorithms from the fast class. We say that:

• •

$P_{i}$* is a Type 1 MSP of $\mathcal{A}$ if, according to $\mathcal{A}$, is computed using an algorithm from the standard class. If the entire problem is to be solved using an algorithm for the standard class, we say that the entire problem is the unique (improper) Type 1 MSP generated by $\mathcal{A}$.*

• •

$P_{i}$* is a Type 2 MSP of $\mathcal{A}$ if, according to $\mathcal{A}$, is computed by generating 7 sub-problems according to the recursive scheme corresponding to an algorithm from the fast (Strassen-like) class, and if the generated sub-problems have input size strictly smaller than $2\sqrt{M}\times 2\sqrt{M}$. If the entire problem uses a recursive algorithm from the fast class to generate 7 sub-problems with input size smaller than $2\sqrt{M}\times 2\sqrt{M}$, we say that the entire problem is the unique, improper, Type 2 MSP generated by $\mathcal{A}$.*

In the following we denote as $\nu_{1}$ (resp., $\nu_{2}$) the number of Type 1 (resp., Type 2) MSPs generated by $\mathcal{A}$. Let $P_{i}$ denote the $i$-th MSP generated by $\mathcal{A}$ and let $G^{\mathcal{A}}_{P_{i}}$ denote the corresponding sub-CDAG of $G^{\mathcal{A}}$. We denote as $\mathbf{A}_{i}$ and $\mathbf{B}_{i}$ (resp., $\mathbf{C}_{i}$) the input factor matrices (resp., the output product matrix) of $P_{i}$.

Properties of MSPs and their corresponding sub-CDAGs:

By Definition 2, we have that for each pair of distinct MMSPs $P_{1}$ and $P_{2}$, $P_{2}$ is not recursively generated by $\mathcal{A}$ in order to compute $P_{1}$ or vice versa. Hence, by Lemma 2, the sub-CDAGs of $G^{\mathcal{A}}$ that correspond each to one of the MSPs generated by $\mathcal{A}$ are vertex, disjoint.

In order to obtain our I/O lower bound for algorithms in $\mathfrak{H}{}$, we characterize properties regarding the minimum dominator size of an arbitrary subset of $\mathcal{Y}{}$ and $\mathcal{Z}{}$.

Lemma 4** (Proof in Appendix A.1).**

Let $G^{\mathcal{A}}$ be the CDAG corresponding to an algorithm $\mathcal{A}\in\mathfrak{H}{}$ which admits $n_{1}$ Type 1 MSPs. For each Type 1 MSP $P_{i}$ let $\mathcal{Y}{}_{i}$ denote the set of input vertices of the associated sub-CDAG $G^{\mathcal{A}}_{P_{i}}$ which correspond each to an entry of the input matrices $\mathbf{A}_{i}$ and $\mathbf{B}_{i}$. Further, we define $\mathcal{Y}{}=\cup_{i=1}^{\nu_{1}}\mathcal{Y}{}_{i}$.

Let $Y\subseteq\mathcal{Y}{}$ in $G^{\mathcal{A}}$ such that $|Y\cap\mathcal{Y}{}_{i}|=a_{i}/\sqrt{b_{i}}$, with $a_{i},b_{i}\in\mathbb{N}$, $a_{i}\geq b_{i}$ for $i=1,2,\ldots,\nu_{1}$, and such that $b_{i}=0$ if and only if $a_{i}=0$.111Here we use as convention that $0/0=0$. Any dominator set $D{}$ of $Y$ satisfies $|D{}|\geq\min\{2M,\sum_{i=1}^{\nu_{1}}a_{i}/\sqrt{\sum_{i=1}^{\nu_{1}}b_{i}}\}$.

Lemma 5** (Proof in Appendix A.2).**

Let $G^{\mathcal{A}}$ be the CDAG corresponding to an algorithm $\mathcal{A}\in\mathfrak{H}{}$ which admits $n_{2}$ Type 2 MSPs. Further let $\mathcal{Z}$ denote the set of vertices corresponding to the entries of the output matrices of the $n_{2}$ Type 2 MSPs. Given any subset $Z\subseteq\mathcal{Z}{}$ in $G^{\mathcal{A}}$ with $|Z|\leq 4M$, any dominator set $D$ of $Z$ satisfies $|D|\geq|Z|/2$.

For each Type 1 MSP $P_{i}$ generated by $\mathcal{A}$, with input size222In general, different Type 1 MSP may have different input sizes $n_{i}\times n_{i}$, we denote as $T_{i}$ the set of variables whose value correspond to the $n_{i}^{3}$ elementary products $\mathbf{A}_{i}[j][k]\mathbf{B}_{i}[k][j]$ for $j,k=0,1,\ldots,n_{i}-1$. Further, we denote as $\mathcal{T}_{i}$ the set of vertices corresponding to the variables in $T_{i}$, and we define $\mathcal{T}=\cup_{i=1}^{\nu_{1}}\mathcal{T}_{i}$.

Lemma 6** (Proof in Appendix A.3).**

For any Type 1 MSPs generated by $\mathcal{A}$ consider $T^{\prime}_{i}\subseteq T_{i}$. Let $\mathcal{Y}{}_{i}^{(\mathbf{A})}\subseteq\mathcal{Y}{}_{i}$ (resp., $\mathcal{Y}{}_{i}^{(\mathbf{B})}\subseteq\mathcal{Y}{}_{i}$) denote a subset of the vertices corresponding to entries of $\mathbf{A}_{i}$ (resp., $\mathbf{B}_{i}$) which are multiplied in at least one of the elementary products in $T^{\prime}_{i}$. Then any dominator $D$ of the vertices corresponding to $T^{\prime}_{i}$ with respect to the the vertices in $\mathcal{Y}{}_{i}$ is such that

[TABLE]

6 I/O lower bounds for algorithm in $\mathfrak{H}{}$ and $\mathfrak{UH}\left(n_{0}{}\right){}$

Theorem 7**.**

Let $\mathcal{A}\in\mathfrak{H}{}$ be an algorithm to multiply two square matrices $\mathbf{A},\mathbf{B}\in\mathcal{R}^{n\times n}$. If run on a sequential machine with cache of size $M$ and such that up to $B$ memory words stored in consecutive memory locations can be moved from cache to slow memory and vice versa using a single memory operation, $\mathcal{A}$’s I/O complexity satisfies:

[TABLE]

for $c=0.38988157484$, where $|\mathcal{T}|$ denotes the total number of internal elementary products computed by the Type 1 MSPs generated by $\mathcal{A}$ and $\nu_{2}$ denotes the total number of Type 2 MSPs generated by $\mathcal{A}$.

If run on $P{}$ processors each equipped with a local memory of size $M<n^{2}$ and where for each I/O operation it is possible to move up to $B_{m}{}$ words, $\mathcal{A}$’s I/O complexity satisfies:

[TABLE]

Proof.

We prove the result in (1) (resp., (2)) for the case $B=1$ (resp., $B_{m}{}=1$). The result then trivially generalizes for a generic $B$ (resp., $B_{m}{}$). We first prove the result for the sequential case in in (1). The bound for the parallel case in (2) will be obtained as a simple generalization.

The fact that $IO_{\mathcal{A}}(n,M,1)\geq 2n^{2}$ follows trivially from the fact that as in our model the input matrices $\mathbf{A}$ and $\mathbf{B}$ are initially stored in slow memory, it will necessary to move the entire input to the cache at least once using at least $2n^{2}$ I/O operations. If $\mathcal{A}{}$ does not generate any MSPs the statement in (1) is trivially verified. In the following, we assume $\nu_{1}+\nu_{2}\geq 1$.

Let $G^{\mathcal{A}}$ denote the CDAG associated with algorithm $\mathcal{A}$ according to the construction in Section 4. By definition, and from Lemma 2, the $\nu_{1}+\nu_{2}$ sub-CDAGs of $G^{\mathcal{A}}$ corresponding each to one of the MSPs generated by $\mathcal{A}$ are vertex-disjoint. Hence, the $\mathcal{T}_{i}$’s are a partition of $\mathcal{T}$ and $|\mathcal{T}|=\sum_{i=1}^{\nu_{1}}|\mathcal{T}_{i}|$.

By Definition 2, the MSP generated by $\mathcal{A}{}$ have input (resp., output) matrices of size greater or equal to $2\sqrt{M}\times 2\sqrt{M}$. Recall that we denote as $\mathcal{Z}{}$ the set of vertices which correspond to the outputs of the $\nu_{2}$ Type 2 MSPs, we have $|\mathcal{Z}{}|\geq 4M\nu_{l}$.

Let $\mathcal{C}$ be any computation schedule for the sequential execution of $\mathcal{A}$ using a cache of size $M$. We partition $\mathcal{C}$ into non-overlapping segments $\mathcal{C}_{1},\mathcal{C}_{2},\ldots$ such that during each $\mathcal{C}_{j}$ either (a) exactly $M^{3/2}$ distinct values 333For simplicity of presentation, we assume $M^{3/2}\in\mathbb{N}^{+}$. corresponding to vertices in $\mathcal{T}$, denoted as $\mathcal{T}^{(j)}$, are explicitly computed (i.e., not loaded from slow memory), or (b) $4M$ distinct values corresponding to vertices in $\mathcal{Z}{}$ (denoted as $\mathcal{Z}{}_{j}$) are evaluated for the first time. Clearly there are at least $\max\{|\mathcal{T}|/{M^{3/2}},\nu_{2}\}$ such segments.

Below we show that the number $g_{j}$ of I/O operations executed during each $\mathcal{C}_{j}$ satisfies $g_{j}\geq cM$ for case (a) and $g_{j}\geq M$ for case (b), from which the theorem follows.

Case (a): For each Type 1 MSP $P_{i}$ let $\mathcal{T}^{(j)}_{i}=\mathcal{T}^{(j)}\cap\mathcal{T}_{i}$. As the $\nu_{1}$ sub-CDAGs corresponding each to one of the Type 1 MSPs are vertex-disjoint, so are the sets $\mathcal{T}_{i}$. Hence, the $\mathcal{T}^{(j)}_{i}$’s constitute a partition of $\mathcal{T}^{(j)}$. Let $\mathbf{A}_{i}$ and $\mathbf{B}_{i}$ (resp., $\mathbf{C}_{i}$) denote the input matrices (resp., output matrix) of $P_{i}$ with $\mathbf{A}_{i},\mathbf{B}_{i},\mathbf{C}_{i}\in\mathcal{R}^{n_{i}\times n_{i}}$, and let $A_{i}$ and $B_{i}$ (resp., $C_{i}$) denote the set of the variables corresponding to the entries of $\mathbf{A}_{i}$ and $\mathbf{B}_{i}$ (resp., $\mathbf{C}_{i}$). Further, we denote as $T_{i}$ the set of values corresponding to the vertices in $\mathcal{T}_{i}$. For $r,s=0,1,\ldots,n_{i}-1$, we say that $\mathbf{C}_{i}[r][s]$ is “active during $\mathcal{C}_{j}$” if any of the elementary multiplications $\mathbf{A}_{i}[r][k]\mathbf{B}_{i}[k][s]$, for $k=0,1,\ldots,n_{i}-1$, correspond to any of the vertices in $\mathcal{T}^{(j)}_{i}$. Further we say that a $\mathbf{A}_{i}[r][s]$ (resp., $\mathbf{B}_{i}[r][s]$) is “accessed during $\mathcal{C}_{j}$” if any of the elementary multiplications $\mathbf{A}_{i}[r][s]\mathbf{B}_{i}[s][k]$ (resp., $\mathbf{A}_{i}[k][r]\mathbf{B}_{i}[r][s]$), for $k=0,1,\ldots,n_{i}-1$, correspond to any of the vertices in $\mathcal{T}^{(j)}_{i}$.

Our analysis makes use of the following property of standard matrix multiplication algorithms:

Lemma 8** (Loomis-Whitney inequality [21, Lemma 2.2]).**

Let $\mathcal{Y}{}_{i,\mathbf{A}}^{\prime}$ (resp., $\mathcal{Y}{}_{i,\mathbf{B}}^{\prime}$) denote the set of vertices corresponding to the entries of $\mathbf{A}_{i}$ (resp., $\mathbf{B}_{i}$) which are accessed during $\mathcal{C}_{j}$, and let $\mathcal{Z}{}_{i}^{\prime}$ denote the set of vertices corresponding to the entries of $\mathbf{C}_{i}$ which are active during $\mathcal{C}_{j}$. Then

[TABLE]

Lemma 8 is a reworked version of a property originally presented by Irony et al. [21, Lemma 2.2], which itself is a consequence of the Loomis-Whitney geometric theorem [25].

Let $\mathbf{C}_{i}[r][s]$ be active during $\mathcal{C}_{j}$. In order to compute $\mathbf{C}_{i}[r][s]$ entirely during $\mathcal{C}_{j}$ (i.e., without using partial accumulation of the summation $\sum_{k=0}^{n_{i}-1}\mathbf{A}_{i}[r][k]\mathbf{B}[k][s]$), it will be necessary to evaluate all the $n_{i}$ elementary products $\mathbf{A}_{i}[r][k]\mathbf{B}_{i}[k][s]$, for $k=0,1,\ldots,n_{i}-1$, during $\mathcal{C}_{j}$ itself. Thus, at most $\lfloor|\mathcal{T}^{(j)}_{i}|/{n_{i}}\rfloor$ entries of $\mathbf{C}_{i}[r][s]$ can be entirely computed during $\mathcal{C}_{j}$.

Let $\mathbf{C}_{i}[r][s]$ denote an entry of $\mathbf{C}_{i}$ which is active but not entirely computed during $\mathcal{C}_{j}$. There are two possible scenarios:

• •

$\mathbf{C}_{i}[r][s]$ is computed during $\mathcal{C}_{j}$: The computation thus requires for a partial accumulation of $\sum_{k=0}^{n_{i}-1}\mathbf{A}_{i}[r][k]\mathbf{B}[k][s]$ to have been previously computed and either held in the cache at the beginning of $\mathcal{C}_{j}$, or to moved to cache using a read I/O operation during $\mathcal{C}_{j}$;

• •

$\mathbf{C}_{i}[r][s]$ is not computed during $\mathcal{C}_{j}$: As $\mathcal{C}$ is a parsimonious computation, the partial accumulation of $\sum_{k=0}^{n_{i}-1}\mathbf{A}_{i}[r][k]\mathbf{B}[k][s]$ obtained from the elementary products computed during $\mathcal{C}_{j}$ must either remain in the cache at the end of $\mathcal{C}_{j}$, or be moved to slow memory using a write I/O operation during $\mathcal{C}_{j}$;

In both cases, any partial accumulation either held in memory at the beginning (resp., end) of $\mathcal{C}_{j}$ or read from slow memory to cache (resp., written from cache to slow memory) during $\mathcal{C}_{j}$ is, by definition, not shared between multiple entries in $\mathbf{C}_{i}$.

Let $G^{\mathcal{A}}_{P_{i}}$ denote the sub-CDAG of $G^{\mathcal{A}}$ corresponding to the Type 1 MSP $P_{i}$. In the following, we refer as $D^{\prime}_{i}$ to the set of vertices of $G^{\mathcal{A}}_{P_{i}}$ corresponding to the values of such partial accumulators. For each of the least $|D{}^{\prime}_{i}|=\max\{0,|\mathcal{Z}_{i}^{\prime}|-|\mathcal{T}^{(j)}_{i}|/n_{i}\}$ entries of $\mathbf{C}_{i}$ which are active but not entirely computed during $\mathcal{C}_{j}$, either one of the entries of the cache must be occupied at the beginning of $\mathcal{C}_{j}$, or one I/O operation is executed during $\mathcal{C}_{j}$. Let $D{}^{\prime}=\cup_{i=1}^{\nu_{1}}D{}^{\prime}_{i}$. As, by Lemma 2, the sub-CDAGs corresponding to the $\nu_{1}$ Type 1 MSPs are vertex disjoint, so are the the sets $D{}^{\prime}_{i}$. Let $\mathcal{Z}^{\prime}=\sum_{i=1}^{\nu_{l}}|\mathcal{Z}_{i}^{\prime}|$. We have:

[TABLE]

where the last passage follows from the fact that, by Definition 2, $n_{i}\geq 2M$.

From Lemma 8, the set of vertices $\mathcal{Y}{}_{i,\mathbf{A}}^{\prime}$ (resp., $\mathcal{Y}{}_{i,\mathbf{B}}^{\prime}$) which correspond to entries of $\mathbf{A}_{i}$ (resp., $\mathbf{B}_{i}$) which are accessed during $\mathcal{C}_{j}$ satisfies $|\mathcal{Y}{}_{i,\mathbf{A}}^{\prime}||\mathcal{Y}{}_{i,\mathbf{B}}^{\prime}|\geq|\mathcal{T}^{(j)}_{i}|^{2}/|\mathcal{Z}_{i}^{\prime}|$. Hence, at least $|\mathcal{Y}{}_{i,\mathbf{A}}^{\prime}|+|\mathcal{Y}{}_{i,\mathbf{B}}^{\prime}|\geq 2|\mathcal{T}^{(j)}_{i}|/\sqrt{|\mathcal{Z}_{i}^{\prime}|}$ entries from the input matrices of $P_{i}$ are accessed during $\mathcal{C}_{j}$. Let $\mathcal{Y}$ denote the set of vertices corresponding to the entries of the input matrices $\mathbf{A}_{i},\mathbf{B}_{i}$ of $P_{i}$. From Lemma 6 we have that there exists a set $\mathcal{Y}{}_{i}^{\prime}\subseteq\mathcal{Y}{}_{i}$, with$|\mathcal{Y}{}_{i}^{\prime}|\geq\max\{|\mathcal{Y}{}_{i,\mathbf{A}}^{\prime}|,|\mathcal{Y}{}_{i,\mathbf{B}}^{\prime}|\}\geq|\mathcal{T}^{(j)}_{i}|/\sqrt{|\mathcal{Z}_{i}^{\prime}|}$, such that the vertices in $\mathcal{Y}{}_{i}^{\prime}$ are connected by vertex disjoint pats to the vertices in $\mathcal{T}^{(j)}_{i}$. Let $Y=\cup_{i=1}^{\nu_{1}}\mathcal{Y}{}_{i}^{\prime}$. As, by Lemma 2, the sub-CDAGs corresponding to the $\nu_{l}$ Type 1 MSPs are vertex disjoint, so are the the sets $\mathcal{Y}{}_{i}^{\prime}$ for $i=1,2,\ldots,\nu_{1}$. Hence

[TABLE]

From Lemma 4 any dominator $D{}_{Y}$ of $Y$, must be such that

[TABLE]

Hence, we can conclude that any dominator $D{}^{\prime\prime}$ of $\mathcal{T}^{(j)}$ must be such that

[TABLE]

Consider the set $D{}$ of vertices of $G_{\mathcal{A}}$ corresponding to the at most $M$ values stored in the cache at the beginning of $\mathcal{C}_{j}$ and to the at most $g_{j}$ values loaded into the cache form the slow memory (resp., written into the slow memory from the cache) during $\mathcal{C}_{j}$ by means of a read (resp., write) I/O operation. Clearly, $|D|\leq M+g_{j}$.

In order for the $M^{3/2}$ values from $\mathcal{T}^{(j)}$ to be computed during $\mathcal{C}_{j}$ there must be no path connecting any vertex in $\mathcal{T}^{(j)}$, and, hence, $Y$, to any input vertex of $G^{\mathcal{A}}$ which does not have at least one vertex in $D{}$, that is $D{}$ has to admit a subset $D{}^{\prime\prime}\subseteq D{}$ such that $D{}^{\prime\prime}$ is a dominator set of $\mathcal{T}^{(j)}$. Note that, as the values corresponding to vertices in $\mathcal{T}^{(j)}$ are actually computed during $C_{J}$ (i.e., not loaded from memory using a read I/O operation), $D{}^{\prime\prime}$ does not include vertices in $\mathcal{T}^{(j)}$ itself. Further, as motivated in the previous discussion, $D{}$ must include all the vertices in the set $D{}^{\prime}$ corresponding to values of partial accumulators of the active output values of Type 1 MSPs during $\mathcal{C}_{j}$.

By construction, $D{}^{\prime}$ and $D{}^{\prime\prime}$ are vertex disjoint. Hence, from (4) and (5) we have:

[TABLE]

As, by construction, $|\mathcal{T}^{(j)}|=M^{3/2}$, we have:

[TABLE]

By studying its derivative after opportunely accounting for the minimum, we have that (6) is minimized for $|\mathcal{Z}^{\prime}|=2^{-2/3}M$. Hence we have: $|D{}|>2^{-2/3}M+2^{1/3}M^{3/2}-M/2=1.38988157484M$. Whence $|D|\leq M+g_{j}$, which implies $g_{j}\geq|D|-M>0.38988157484M$, as stated above.

Case (b): In order for the $4M$ values from $\mathcal{Z}{}_{j}$ to be computed during $\mathcal{C}_{j}$ there must be no path connecting any vertex in $\mathcal{Z}{}_{j}$ to any input vertex of $G_{\mathcal{A}}$ which does not have at least one vertex in $D_{j}$, that is $D_{j}$ has to be a dominator set of $\mathcal{Z}{}_{j}$. From Lemma 5, any dominator set $D$ of any subset $Z\subseteq\mathcal{Z}{}$ with $|Z|\leq 4M$ satisfies $|D|\geq|Z|/2$, whence $M+g_{i}\geq|D_{i}|\geq|\mathcal{Z}{}_{j}|/2=2M$, which implies $g_{j}\geq M$ as stated above. This concludes the proof for the sequential case in (1).

The proof for the bound for the parallel model in  (2), follows from the observation that at least one of the $P{}$ processors, denoted as $P{}^{*}$, must compute at least $|\mathcal{T}|/P{}$ values corresponding to vertices in $\mathcal{T}$ or $|\mathcal{Z}{}|/P{}$ values corresponding to vertices in $\mathcal{Z}{}$ (or both). The bound follows by applying the same argument discussed for the sequential case to the computation executed by $P^{*}$. ∎

Note that if $\mathcal{A}$ is such that the product is entirely computed using an algorithm from the standard class (resp., a fast matrix multiplication algorithm), the bounds of Theorem 7 corresponds asymptotically to the results of [18] for the sequential case and [21] for the parallel case (resp., the results in [11]).

The bound in (2) can be further generalized to a slightly different model in which each of the $P{}$ processors is equipped with a cache memory of size $M$ and a slow memory of unbounded size. In such case, the I/O complexity of the algorithms in $\mathfrak{H}{}$ corresponds to the total number of both messages received and the number of I/O operations used to move data from cache to slow memory and vice versa (i.e., read and write) executed by at least one of the $P{}$ processors.

I/O lower bound for uniform, non stationary algorithms in $\mathfrak{UH}\left(n_{0}{}\right){}$:

For the sub-class of uniform, non stationary algorithms $\mathfrak{UH}\left(n_{0}{}\right){}$, given the values of $n$, $M$ and $n_{0}{}$ is possible to compute a closed form expression for the values of $\nu_{1},\nu_{2}$ and $|\mathcal{T}|$.444The constants terms in Theorem 9 assume that $n,n_{0}{}$ and $M$ are powers of two. If that not the case the statement holds with minor adjustments. Then, by applying Theorem 7 we have:

Theorem 9**.**

Let $\mathcal{A}\in\mathfrak{UH}\left(n_{0}{}\right){}$ be an algorithm to multiply two square matrices $\mathbf{A},\mathbf{B}\in\mathcal{R}^{n\times n}$. If run on a sequential machine with cache of size $M$ and such that up to $B$ memory words stored in consecutive memory locations can be moved from cache to slow memory and vice versa using a single memory operation, $\mathcal{A}$’s I/O complexity satisfies:

[TABLE]

If run on $P{}$ processors each equipped with a local memory of size $M<n^{2}$ and where for each I/O operation it is possible to move up to $B_{m}{}$ memory words, $\mathcal{A}$’s I/O complexity satisfies:

[TABLE]

Proof.

The statement follows by bounding the values $\nu_{1}$,$\nu_{2}$ and $|\mathcal{T}|$ for $\mathcal{A}\in\mathfrak{UH}\left(n_{0}{}\right){}$, and by applying the general result from Theorem 7. In order to simplify the presentation, in the following we assume that the values $n,n_{0}{},M$ are powers of two. If that is not the case, the theorem holds with some minor adjustments to the constant multiplicative factor.

Let let $i$ be the smallest value in $\mathbb{N}$ such that $n/2^{i}=\max\{n_{0}{},2\sqrt{M}\}$. By definition of $\mathfrak{H}{}$, at each of the $i$ recursive levels $\mathcal{A}$ generates $7^{i}$ sub-problems of size $n/2^{i}$.

• •

If $n_{0}{}>2\sqrt{M}$, $\mathcal{A}$ generates

[TABLE]

Type 1 MSP each with input size $n_{0}{}\times n_{0}{}$. As, by Definition 2, the Type 1 MSP are input dijoint we have:

[TABLE]

• •

Otherwise, if $n_{0}{}\leq 2\sqrt{M}$, $\mathcal{A}$ generates

[TABLE]

Type 2 MSP each with input size $2\sqrt{M}\times 2\sqrt{M}$.

The statement then follows by applying the result in Theorem 7. ∎

Theorem 9 extends the result by Scott [32] by expanding the class of hybrid matrix multiplication algorithms being considered (e.g., it does not limit the class of standard matrix multiplication to the divide and conquer algorithm based on block-partitioning), and by removing the assumption that no intermediate value may be recomputed.

On the tightness of the bound:

An opportune composition of the cache-optimal version of the Strassen’s algorithm [35] (as discussed in  [3]) with the standard cache-optimal divide and conquer algorithm for square matrix multiplication based on block-partitioning [13] leads to a sequential hybrid algorithms in $\mathfrak{H}{}$ (resp., $\mathfrak{UH}\left(n_{0}{}\right){}$) whose I/O cost asymptotically matches the I/O complexity lower bounds in Theorem 7 (1) (resp., Theorem 9 (7)).

Parallel algorithms in $\mathfrak{H}{}$ (resp., $\mathfrak{UH}\left(n_{0}{}\right){}$) asymptotically matching the I/O lower bounds in for the parallel case in Theorem 7 (2) (resp., Theorem 9 (8)) can be obtained by composing the communication avoiding version of Strassen’s algorithm by Ballard et al. [3] with the communication avoiding “2.5” standard algorithm by Solomonik and Demmel [34].

Hence, the lower bounds in Theorem 7 and Theorem 9 are asymptotically tight and the mentioned algorithms form $\mathfrak{H}{}$ and $\mathfrak{UH}\left(n_{0}{}\right){}$ whose I/O cost asymptotically match the lower bounds are indeed I/O optimal.

Further, as the mentioned I/O optimal algorithms from $\mathfrak{H}{}$ and $\mathfrak{UH}\left(n_{0}{}\right){}$ do not recompute any intermediate value 555We do not consider replication of the input used by the mentioned parallel algorithms as recomputation, but rather as a repeated access to the input values., we can conclude that using recomputation may lead to at most a constant factor reduction of the I/O cost of hybrid algorithms in $\mathfrak{H}{}$ and $\mathfrak{UH}\left(n_{0}{}\right){}$.

Generalization to fast matrix multiplication model with base other than $2\times 2$:

The general statement of Theorem 7 can be extended to by enriching $\mathfrak{H}{}$ to include any fast Strassen-like algorithm with base case other than $2\times 2$ provided that the associated encoder CDAG satisfies properties equivalent to those expressed by Lemma 2(i.e., the input disjointedness of the sub-problems generated at each recursive step) and Lemma 3 (i.e., the connectivity between input and output of the encoder CDAGs via vertex disjoint paths) for the $2\times 2$ base. If these properties hold, so does the general structure of Theorem 7, given an opportune adjustment of the definition of maximal sub-problem.

7 Conclusion

This work introduced the first characterization of the I/O complexity of hybrid matrix multiplication algorithms combining fast Strassen-like algorithms with standard algorithms. We established asymptotically tight lower bounds that hold even when recomputation is allowed. The generality of the technique used the analysis makes it promising for the analysis of other hybrid recursive algorithms, e.g., for hybrid algorithms for integer multiplication [12].

Our results contribute to the study of the effect of recomputation with respect to the I/O complexity of CDAG algorithms. While we are far from a characterization of those CDAGs for which recomputation is effective, this broad goal remains a fundamental challenge for any attempt toward a general theory of the communication requirements of computations.

Acknowledgments

The author would like to thank Gianfranco Bilardi at the University of Padova for the helpful suggestions and discussions.

Appendix A Proofs of technical lemmas

A.1 Proof of Lemma 4

The proof of Lemma 4 uses an analysis similar to that in [11](Lemma 6), albeit with several important variations. Before delving into the details of the proof of Lemma  4 we introduce the following technical lemma:

Lemma 10**.**

Let $a_{i},a_{2},\ldots a_{i}\in\mathbb{N}$ (resp., $b_{i},b_{2},\ldots b_{i}\in\mathbb{N}$) such that $b_{j}=0$ if and only if $a_{i}=0$ and such that:

[TABLE]

using the convention that $0/0=0$. Then we have:

[TABLE]

Proof.

The proof is by induction on the value of $i$. In the base case case $i=1$, hence $a=a_{1}$ and $b=b_{1}$. Thus the statement is trivially verified. We assume inductively that the statement holds for $i\geq 1$ and we show that it holds also for $i+1$. Let

[TABLE]

be can therefore apply the inductive hypothesis to the first $i$ elements and obtain:

[TABLE]

Let $a^{\prime}=\sum_{j=1}^{i}a_{j}$ and $b^{\prime}=\sum_{j=1}^{i}b_{j}$. If $a_{i+1}=0$ the statement is trivially verified. In the following we consider the case for which $a_{i+1}>0$ and, by assumption, $b_{i+1}>0$.

From (10), we have that the lemma is guaranteed to be verified if the following holds:

[TABLE]

Clearly (11) is verified for $b^{\prime}\geq b_{i+1}$. Since, by assumption, $b=b^{\prime}+b_{i+1}$, we can therefore conclude that the statement is verified for $b^{\prime}\geq 3b/4$.

In the following we consider the case for $b^{\prime}<3b/4$. Let $a^{\prime}=ax$ (resp., $b^{\prime}=by$), with $x\in(0,1]$. By assumption $x,y\neq 0$. Further, as we are considering the case such that $b^{\prime}<3b/4$, we have that $y\in(0,3/4)$. Since, by assumption, $a=a^{\prime}+a_{i+1}$ and $b=b^{\prime}+b_{i+1}$, we have:

[TABLE]

Therefore, the Lemma is clearly verified if the following holds:

[TABLE]

By multiplying both sides of the previous inequality by $2\sqrt{y}\sqrt{y-1}$ (which is clearly greater than [math]), we have that (12) holds if:

[TABLE]

which in turn holds if:

[TABLE]

As $y\in(0,3/4)$ we have $2\sqrt{1-y}>1>\sqrt{y}$ and $\left(x-\sqrt{y}\right)\left(x-1\right)$. Hence (14) holds ant the Lemma follows. ∎

Recall that given an algorithm $\mathcal{A}\in\mathfrak{H}{}$ used to multiply input squared matrices $\mathbf{A},\mathbf{B}\in\mathcal{R}^{n\times n}$, we denote as $G^{\mathcal{A}}$ the corresponding CDAG constructed according to the description in Section 4. Further, we denote as $\mathcal{X}{}$ the set of input vertices of $G^{\mathcal{A}}$. That is the set of vertices corresponding each to the entries in the input matrices $\mathbf{A}$ and $\mathbf{B}$. We now present the proof of Lemma 4.

Proof of Lemma 4.

The proof proceeds by induction on the number of Type 1 MSPs $\nu_{1}$ generated by $\mathcal{A}$. In the base case $\nu_{1}=1$ and, by Definition 2 the entire problem is the only, improper Type 1 MSP generated by $\mathcal{A}$. Hence, the sets $\mathcal{Y}{}$ and $\mathcal{X}{}$ coincide and the statement is trivially verified.

Assuming now inductively that the statement holds $\nu_{1}=k\geq 1$, we shall show it also holds for $\nu_{1}=k+1$.

As $\nu_{1}>1$, we have that the algorithm executes at least one recursive step. Let $P^{(j)}$ denote the 7 sub-problems generated at the first recursion step. We distinguish two cases: (a) at least two of of the seven sub-problems $P^{(j)}$ generate each at least one Type 1 MSP; (b) only one of the seven sub-problems generates all $\nu_{1}$ Type 1 MSPs. We first address case (a), as case (b) will follow from a simple extension.

For case (a), let $G^{(j)}$, for $j=1,2,\ldots,7$ denote the seven sub-CDAGs of $G^{\mathcal{A}}$, each corresponding to one of the seven sub-problems generated in the first recursive step of $\mathcal{A}$ according to the chosen Strassen-like scheme as discussed in Section 3. By Definition 2 and as, by assumption, $\mathcal{A}$ generates Type 1 MSPs, we have that each of these seven sub-problems has input size greater or equal to $2\sqrt{M}\times 2\sqrt{M}$. Further, each of the seven sub-problems $P^{(j)}$ generates at most $\nu_{1}-1$ Type 1 MSPs.

Let $Y^{(j)}$ (resp., $\mathcal{Y}{}^{(j)}$) denote the subsets of $Y$ (resp., $\mathcal{Y}{}$) in $G^{(j)}$, for $j=1,2,\ldots,7$. That is $Y^{(j)}=Y\cap\mathcal{Y}^{(j)}$, By Lemma 2 the $G^{(j)}$’s have distinct input values and, hence, are pairwise vertex-disjoint sub-CDAGs of $G^{\mathcal{A}}$. Thus, the $Y^{(1)},Y^{(2)},\ldots,Y^{(7)}$ partition $Y$ (resp., $\mathcal{Y}{}^{(1)},\mathcal{Y}{}^{(2)},\ldots,\mathcal{Y}{}^{(7)}$ partition $\mathcal{Y}{}$) and $\sum_{j=1}^{7}|Y^{(j)}|=|Y|$ (resp., $\sum_{j=1}^{7}|\mathcal{Y}{}^{(j)}|=|\mathcal{Y}{}|$). Let $|Y^{(j)}|=a_{j}/\sqrt{b_{j}}$, we have $|Y|=\sum_{j=1}^{7}a_{j}/\sqrt{b_{j}}$.

By the inductive hypothesis, any dominator set $D^{(j)}$ of $Y^{(j)}$ with respect to the set $\mathcal{K}^{(j)}$ composed by the input vertices of $G^{(j)}$ must be such that $|D^{(j)}|\geq\min\{2M,a_{j}/\sqrt{b}_{j}\}$. By Definition 1 this implies that vertices in $Y^{(j)}$ can be connected to a subset $K^{(j)}\subseteq\mathcal{K}^{(j)}$ of the input vertices of $G^{(j)}$ such that $|K^{(j)}|\geq a_{j}/\sqrt{b}_{j}$, using vertex-disjoint paths. Since the sub-CDAGs $G^{(j)}$ are vertex disjoint, so are the paths connecting vertices in $Y^{(j)}$ to vertices in $K^{(j)}$. In the following we show that it is indeed possible to extended at least $\min\{2M,\sum_{j=1}^{7}a_{j}/\sqrt{\sum_{j=1}^{7}b_{j}}\}$ of these paths to vertices in $\mathcal{X}{}$ while maintaining them vertex disjoint.

According to the construction of $G^{\mathcal{A}}$ as discussed in Section 4, vertices in $\mathcal{X}{}$ corresponding to the entries of input matrix $\mathbf{A}$ (resp., $\mathbf{B}$) are connected to vertices in $\mathcal{K}^{(1)},\mathcal{K}^{(2)},\ldots,\mathcal{K}^{(7)}$ (and, hence, $K^{(1)},K^{(2)},\ldots,K^{(7)}$) by means of $n^{2}$ encoding sub-CDAGs $Enc_{A}$ (resp., $Enc_{B}$). None of these $2n^{2}$ encoding sub-CDAGs share any input or output vertices. No two output vertices of the same encoder sub-CDAG belong to the same sub-CDAG $G^{(j)}$, for $j=1,2,\ldots,7$. This fact ensures that for a single sub-CDAG $G^{(j)}$, for $j=1,2,\ldots,7$, it is possible to connect all the vertices in $\mathcal{K}^{(j)}$ (and, hence, $K^{(j)}$) to a subset of the vertices in $\mathcal{X}{}$ via vertex disjoint paths.

For each of the $2n^{2}$ encoder sub-CDAGs, let us consider the vector $\mathbf{y}_{l}\in\{0,1\}^{7}$ such that $\mathbf{y}_{l}[j]=1$ iff the corresponding $j$-th output vertex of the encoder, which is an input of $G^{(j)}$, is in $K^{(j)}$. Therefore $|\mathbf{y}_{l}|$ equals the number of output vertices of the $l$-th encoder sub-CDAG which are in $K$. From Lemma 3, for each encoder sub-CDAG there exists a subset $X_{l}\in\mathcal{X}$ of the input vertices of the $l$-th encoder sub-CDAG for which it is possible to connect each vertex in $X_{l}$ to a distinct output vertex of the $l$-th encoder sub-CDAG using vertex disjoint paths, each constituted by a singular edge with $\min\{|\mathbf{y}_{l}|,1+\lceil\left(|\mathbf{y}_{l}|-1\right)/2\rceil\}\leq|X_{l}|\leq|\mathbf{y}_{l}|$. The number of vertex disjoint paths connecting vertices in $\mathcal{X}{}$, to vertices in $\cup_{j=1}^{7}K^{(j)}$ is therefore at least $\sum_{l=1}^{2n^{2}}\min\{|\mathbf{y}_{l}|,1+\lceil\left(|\mathbf{y}_{l}|-1\right)/2\rceil\}$, under the constraint that $\sum_{l=1}^{2n^{2}}\mathbf{y}_{l}[j]=a_{j}/\sqrt{b_{j}}$, for $j=1,2,\ldots,7$.

Let us assume w.l.o.g. that $a_{1}/\sqrt{b_{1}}\geq a_{2}/\sqrt{b_{2}}\geq\ldots\geq a_{7}/\sqrt{b_{7}}$. As previously stated, it is possible to connect all vertices in $K_{1}$ to vertices in $\mathcal{X}{}$ through vertex disjoint paths. Consider now all possible dispositions of the vertices in $\cup_{j=2}^{7}K^{(j)}$ over the outputs of the $2n^{2}$ encoder sub-CDAGs. Hence, if $a_{1}/\sqrt{b_{1}}\geq 2M$ we have that there are therefore at least $M$ vertex disjoint paths connecting vertices in $\mathcal{X}{}$ to vertices in $K_{1}$, and, thus, to vertices in $Y$ as desired. In the following we assume $a_{1}/\sqrt{b_{1}}<M$.

Recall that the output vertices of an encoder sub-CDAG belong each to a different sub-CDAG $G^{(j)}$. From Lemma 3 we have that for each encoder there exists a subset $X_{l}\subset{X}$ of the input vertices of the $l$-th encoder sub-CDAG, with |X_{l}|\geq\min\Big{\{}|\mathbf{y}_{l}|,1+\left\lceil\left(|\mathbf{y}_{l}|-1\right)/2\right\rceil\Big{\}}\geq\mathbf{y}_{l}[1]+\left(\sum_{j=2}^{7}\mathbf{y}_{l}[j]\right)/2, for which is possible to connect all vertices in $X_{l}$ to $|X_{l}|$ distinct output vertices of the $l$-th encoder sub-CDAG which are in $\cup_{j=1}^{7}K^{(j)}$ using $|X_{l}|$ vertex disjoint paths. As all the $Enc$ sub-CDAGs are vertex-disjoint, we can add their contributions so that the number of vertex-disjoint paths connecting vertices in $\mathcal{X}{}$ to vertices in $\cup_{j=1}^{7}K^{(j)}$ is at least

[TABLE]

where the last passage follows by applying Lemma 10. There are therefore at least $a/\sqrt{b}$ vertex disjoint paths connecting vertices in $\mathcal{X}{}$ to vertices in $\cup_{j=1}^{7}K^{(j)}$, and, thus, to vertices in $Y$ as desired. This concludes the proof for case (a).

For case (b), only one of the seven sub-problems $P^{(j)}$ generates all $\nu_{1}$ Type 1 MSPs. Without loss of generality, let $P^{(1)}$ denote such sub-problem and let $G^{\mathcal{A}_{P^{(1)}}}$ be the corresponding sub-CDAG. According to the construction of $G^{\mathcal{A}}$ as discussed in Section 4, vertices in $\mathcal{X}{}$ corresponding to the entries of input matrix $\mathbf{A}$ (resp., $\mathbf{B}$) are connected to the input vertices of $G^{\mathcal{A}_{P^{(1)}}}$, by means of $n^{2}$ encoding sub-CDAGs $Enc_{A}$ (resp., $Enc_{B}$). None of these $2n^{2}$ encoding sub-CDAGs share any input or output vertices. No two output vertices of the same encoder sub-CDAG belong to the same sub-CDAG $G^{(j)}$, for $j=1,2,\ldots,7$. This fact ensures that it is possible to connect all the input vertices of $G^{\mathcal{A}_{P^{(1)}}}$ to a subset of the vertices in $\mathcal{X}{}$ via vertex disjoint paths. The proof for case (b) then follows by recursively applying the arguments in the Proof of Lemma 4 to $G^{\mathcal{A}_{P^{(1)}}}$. ∎

A.2 Proof of Lemma 5

Given an algorithm $\mathcal{A}\in\mathfrak{H}{}$ used to multiply input squared matrices $\mathbf{A},\mathbf{B}\in\mathcal{R}^{n\times n}$, let $G^{\mathcal{A}}$ denote the corresponding CDAG constructed according to the description in Section 4. In the following we denote as $\mathcal{X}{}$ the set of input vertices of $G^{\mathcal{A}}$. That is the set of vertices corresponding each to the entries in the input matrices $\mathbf{A}$ and $\mathbf{B}$. Further, for each Type 2 MSP $P_{i}$ we denote as $\mathcal{Y}{}_{i}$ the set of input vertices of the sub-CDAG $G^{\mathcal{A}}_{P_{i}}$ associated with $P_{i}$. That is the set of vertices corresponding each to the entries in the input matrices $\mathbf{A}_{i}$ and $\mathbf{B}_{i}$ of $P_{i}$. Also, we define $\mathcal{Y}{}=\cup_{i=1}^{\nu_{2}}\mathcal{Y}{}_{i}$.

In order to simplify the presentation of the proof of Lemma 5 we first introduce Lemma 11 which is a heavily modified version of a result previously introduced in [11, Lemma 6].

Lemma 11**.**

Consider an algorithm $\mathcal{A}\in\mathfrak{H}{}$ with input matrices $\mathbf{A},\mathbf{B}\in\mathcal{R}{}^{n\times n}$ which generates $\nu_{2}$ Type 2 MSPs. Given the corresponding CDAG $G^{\mathcal{A}}$, let $Q{}$ be a set of internal (i.e., not input) vertices of its $\nu_{2}$ sub-CDAGs corresponding each to one of the generated Type 2 MSPs. For any $Z\subseteq\mathcal{Z}{}$ with $|Z|\geq 2|Q{}|$ there exist $X\subseteq\mathcal{X}{}$ with $|X|\geq 2\sqrt{M\left(|Z|-2|Q{}|\right)}$ such that each vertex in $X$ is connected to some vertex in $Z$ by a directed path with no vertex in $Q{}$.

Proof.

The proof proceeds by induction on the number $\nu_{2}$ of Type 2 MSPs generated by $\mathcal{A}$. In the base case $\nu_{2}=1$ and, by Definition 2 the entire problem is the only, improper, Type 2 MSP generated by $\mathcal{A}$ and, thus, the sets $\mathcal{Y}{}$ and $\mathcal{X}{}$ coincide. Consider now the following lemma:

Lemma 12** ([11, Lemma 5]).**

Let $G^{n\times n}$ denote the CDAG corresponding to the execution of an unspecified algorithm for the square matrix multiplication function with input matrices of size $n\times n$. Let $O^{\prime}\subseteq O$ be a subset of its output vertices $O$. For any subset $D{}$ of the vertices of $G^{n\times n}$ with $|O^{\prime}|\geq 2|D{}|$, there exists a set $I^{\prime}\subseteq I$ of the input vertices $I$ of $G^{n\times n}$ with cardinality $|I^{\prime}|\geq 2n\sqrt{|O^{\prime}|-2|D{}|}$, such that all vertices in $I^{\prime}$ are connected to some vertex in $O^{\prime}$ by directed paths with no vertex in $D{}$.

As, by Definition 2, $G^{\mathcal{A}}$ is a $G^{n\times n}$ CDAG, with $n\geq 2\sqrt{M}$, we can conclude that the statement of Lemma 11 is verified in the base case as a consequence of Lemma 12 and of the fact that, as previously mentioned, $\mathcal{Y}{}$ and $\mathcal{X}{}$ coincide,

Lemma 12 was originally introduced in [11, Lemma 5] and is based on the analysis of the Grigoriev’s flow of the matrix multiplication function. For the sake of completeness we present the proof in Appendix A.4.

Assuming now inductively that the statement holds for $\nu_{2}=k>1$, we shall show it also holds for $\nu_{2}=k+1$. As $\nu_{2}>1$, we have that the algorithm executes at least one recursive step. Let $P^{(j)}$ denote the 7 sub-problems generated at the first recursion step. We distinguish two cases: (a) at least two of of the seven sub-problems $P^{(j)}$ generate each at least one Type 2 MSP; (b) only one of the seven sub-problems generates all $\nu_{2}$ Type 2 MSPs. We first address case (a), as case (b) will follow from a simple extension.

For case (a), let $G^{(j)}$, for $j=1,2,\ldots,7$ denote the seven sub-CDAGs of $G^{\mathcal{A}}$, each corresponding to one of the seven sub-problems generated in the first recursive step of $\mathcal{A}$ according to the chosen Strassen-like scheme as discussed in Section 4. By Definition 2 and as, by assumption, $\mathcal{A}$ generates Type 2 MSPs, we have that $n/2\geq 2\sqrt{M}$. Further, each of the seven sub-problems $P^{(j)}$ generates at most $\nu_{2}-1$ Type 2 MSPs.

Let $Z^{(j)}$, $\mathcal{Y}{}^{(j)}$ and $Q{}^{(j)}$ respectively denote the subsets of $Z$, $\mathcal{Y}{}$ and $Q{}$ in $G^{(j)}$, for $j=1,2,\ldots,7$. By Lemma 2 the $G^{(j)}$’s have distinct input values and, hence, are pairwise vertex-disjoint sub-CDAGs of $G^{\mathcal{A}}$. Thus, $Z_{1},Z_{2},\ldots,Z_{7}$ partition $Z$, $\mathcal{Y}{}_{1},\mathcal{Y}{}_{2},\ldots,\mathcal{Y}{}_{7}$ partition $\mathcal{Y}{}$ and $Q{}_{1},Q{}_{2},\ldots,Q{}_{7}$ partition $Q{}$. This implies $\sum_{j=1}^{7}|Z^{(j)}|=|Z|$ and $\sum_{j=1}^{7}|Q{}^{(j)}|=|Q{}|$. Let $\delta^{(j)}=\max\{0,|Z^{(j)}|-2|Q{}^{(j)}|\}$, we have $\delta=\sum_{j=1}^{7}\delta^{(j)}\geq|Z|-2|Q{}|$.

Applying the inductive hypothesis to each $G^{(j)}$, we have that there is a subset $Y^{(j)}\subseteq\mathcal{Y}{}^{(j)}$ with $|Y^{(j)}|\geq 4\sqrt{M\delta^{(j)}}$ such that vertices of $Y^{(j)}$ are connected to vertices in $Z^{(j)}$ via paths with no vertex in $Q{}^{(j)}$. In the sequel the set $Y$ referred to in the statement will be identified as a suitable subset of $\cup_{i=j}^{7}Y^{(j)}$ so that property (b) will be automatically satisfied. Towards property (a), we observe by the inductive hypothesis that vertices in $Y^{(j)}$ can be connected to a subset $K^{(j)}$ of the input vertices of $G^{(j)}$ with $|K^{(j)}|=|Y^{(j)}|$, using vertex-disjoint paths. Since the sub-CDAGs $G^{(j)}$ are vertex disjoint, so are the paths connecting vertices in $Y^{(j)}$ to vertices in $K^{(j)}$. It remains to show that at least $4\sqrt{M\left(|Z|-2|Q{}|\right)}$ of these paths can be extended to $\mathcal{X}{}$ while maintaining them vertex-disjoint.

According to the construction of $G^{\mathcal{A}}$ as discussed in Section 4, vertices in $\mathcal{X}{}$ corresponding to the entries of input matrix $\mathbf{A}$ (resp., $\mathbf{B}$) are connected to vertices in $K_{1},K_{2},\ldots,K_{7}$ by means of $n^{2}$ encoding sub-CDAGs $Enc_{A}$ (resp., $Enc_{B}$). None of these $2n^{2}$ encoding sub-CDAGs share any input or output vertices. No two output vertices of the same encoder sub-CDAG belong to the same sub-CDAG $G^{(j)}$, for $j=1,2,\ldots,7$. This fact ensures that for a single sub-CDAG $G^{(j)}$, for $j=1,2,\ldots,7$, it is possible to connect all the vertices in $K_{i}$ to a subset of the vertices in $\mathcal{X}{}$ via vertex disjoint paths.

For each of the $2n^{2}$ encoder sub-CDAGs, let us consider the vector $\mathbf{y}_{l}\in\{0,1\}^{7}$ such that $\mathbf{y}_{l}[j]=1$ iff the corresponding $j$-th output vertex of the encoder, which is an input of $G^{(j)}$, is in $K^{(j)}$. Therefore $|\mathbf{y}_{l}|$ equals the number of output vertices of the $l$-th encoder sub-CDAG which are in $K$. From Lemma 3, for each encoder sub-CDAG there exists a subset $X_{l}\in\mathcal{X}$ of the input vertices of the $l$-th encoder sub-CDAG for which it is possible to connect each vertex in $X_{l}$ to a distinct output vertex of the $l$-th encoder sub-CDAG using vertex disjoint paths, each constituted by a singular edge with $\min\{|\mathbf{y}_{l}|,1+\lceil\left(|\mathbf{y}_{l}|-1\right)/2\rceil\}\leq|X_{l}|\leq|\mathbf{y}_{l}|$. The number of vertex disjoint paths connecting vertices in $\mathcal{X}{}$, to vertices in $\cup_{j=1}^{7}K^{(j)}$ is therefore at least $\sum_{l=1}^{2n^{2}}\min\{|\mathbf{y}_{l}|,1+\lceil\left(|\mathbf{y}_{l}|-1\right)/2\rceil\}$, under the constraint that $\sum_{l=1}^{2n^{2}}\mathbf{y}_{l}[j]=4\sqrt{M\delta^{(j)}}$, for $j=1,2,\ldots,7$. Let us assume w.l.o.g. that $\delta^{(1)}\geq\delta^{(20}\geq\ldots\geq\delta^{(7)}$. As previously stated, it is possible to connect all vertices in $K_{1}$ to vertices in $\mathcal{X}{}$ through vertex disjoint paths. Consider now all possible dispositions of the vertices in $\cup_{j=2}^{7}K^{(j)}$ over the outputs of the $2n^{2}$ encoder sub-CDAGs. Recall that the output vertices of an encoder sub-CDAG belong each to a different sub-CDAG $G^{(j)}$. From Lemma 3 we have that for each encoder there exists a subset $X_{l}\subset{X}$ of the input vertices of the $l$-th encoder sub-CDAG, with |X_{l}|\geq\min\Big{\{}|\mathbf{y}_{l}|,1+\left\lceil\left(|\mathbf{y}_{l}|-1\right)/2\right\rceil\Big{\}}\geq\mathbf{y}_{l}[1]+\left(\sum_{j=2}^{7}\mathbf{y}_{l}[j]\right)/2, for which is possible to connect all vertices in $X_{l}$ to $|X_{l}|$ distinct output vertices of the $l$-th encoder sub-CDAG which are in $\cup_{j=1}^{7}K^{(j)}$ using $|X_{l}|$ vertex disjoint paths. As all the $Enc$ sub-CDAGs are vertex disjoint, we can add their contributions so that the number of vertex disjoint paths connecting vertices in $\mathcal{X}{}$ to vertices in $\cup_{j=1}^{7}K^{(j)}$ is at least $|K_{1}|+\frac{1}{2}\sum\limits_{j=2}^{7}|K^{(j)}|=4\sqrt{M}\left(\sqrt{\delta^{(1)}}+\frac{1}{2}\sum\limits_{j=2}^{7}\sqrt{\delta^{(j)}}\right)$. Squaring this quantity leads to:

[TABLE]

As, by assumption, $\delta^{(1)}\geq\ldots\delta^{(7)}$, we have: $\sqrt{\delta^{(1)}}\sqrt{\delta^{(j)}}\geq\delta^{(j)}$ for $j=2,\ldots,7$. Thus:

[TABLE]

There are therefore at least $4\sqrt{M\left(|Z|-2|Q{}|\right)}$ vertex disjoint paths connecting vertices in $\mathcal{X}{}$ to vertices in $\cup_{j=2}^{7}K^{(j)}$ as desired. This concludes the proof for case (a).

For case (b), only one of the seven sub-problems $P^{(j)}$ generates all $\nu_{1}$ Type 1 MSPs. Without loss of generality, let $P^{(1)}$ denote such sub-problem and let $G^{\mathcal{A}_{P^{(1)}}}$ be the corresponding sub-CDAG. According to the construction of $G^{\mathcal{A}}$ as discussed in Section 4, vertices in $\mathcal{X}{}$ corresponding to the entries of input matrix $\mathbf{A}$ (resp., $\mathbf{B}$) are connected to the input vertices of $G^{\mathcal{A}_{P^{(1)}}}$, by means of $n^{2}$ encoding sub-CDAGs $Enc_{A}$ (resp., $Enc_{B}$). None of these $2n^{2}$ encoding sub-CDAGs share any input or output vertices. No two output vertices of the same encoder sub-CDAG belong to the same sub-CDAG $G^{(j)}$, for $j=1,2,\ldots,7$. This fact ensures that it is possible to connect all the input vertices of $G^{\mathcal{A}_{P^{(1)}}}$ to a subset of the vertices in $\mathcal{X}{}$ via vertex disjoint paths. The proof for case (b) then follows by recursively applying the arguments in the Proof of Lemma 11 to $G^{\mathcal{A}_{P^{(1)}}}$. ∎

Lemma 11, provides the base for the proof of Lemma 5, which is itself a reworked version a result from [11] (Lemma 7) modified in order to account for the hybrid nature of algorithms being considered in this work.

Proof of Lemma 5.

Suppose for contradiction that $D$ is a dominator set for $Z$ in $G^{\mathcal{A}}$ such that $|D|\leq 2M-1$. Let $D^{\prime}\subseteq D$ be the subset of the vertices of $D$ composed by vertices which are not internal to the sub-CDAGs corresponding to the Type 2 MSPs generated, by assumption, by $\mathcal{A}$. From Lemma 11, with $Q=D\setminus D^{\prime}$, there exist $X\subseteq\mathcal{X}{}$ and $Y\subseteq\mathcal{Y}$ with $|X|=|Y|\geq 4\sqrt{M\left(|Z|-2\left(|D|-|D^{\prime}|\right)\right)}$ such that vertices in $X$ are connected to vertices in $Y$ by vertex-disjoint paths. Hence, each vertex in $D^{\prime}$ can be on at most one of these paths. Thus, there exists $X^{\prime}\subseteq X$ and $Y^{\prime}\subseteq Y$ with $|X^{\prime}|=|Y^{\prime}|\geq\phi=4\sqrt{M\left(|Z|-2\left(|D|-|D^{\prime}|\right)\right)}-|D^{\prime}|$ paths from $X^{\prime}$ to $Y^{\prime}$ with no vertex in $D^{\prime}$. From Lemma 11, we also have that all vertices in $Y$, and, hence, in $Y^{\prime}$, are connected to some vertex in $Z$ by a path with no vertex in $D\setminus D^{\prime}$. Thus, there are at least $\phi$ paths connecting vertices in $X^{\prime}\subseteq\mathcal{X}$ to vertices in $Z$ with no vertex in $D$. We shall now show that the contradiction assumption $|D|\leq 2M-1$ implies $\phi>0$:

[TABLE]

By $|D|\leq 2M-1$, we have $|Z|-2|D|>4M-2(M-1)>0$. Furthermore, from $D^{\prime}\subseteq D$, we have $32M>2M-1>|D|\geq|D^{\prime}|$. Therefore:

[TABLE]

Again, $|D|\leq 2M-1$ implies $M\left(|Z|-2\left(|D|-|D^{\prime}|\right)\right)>0$. Hence, we can take the square root on both sides of (15) and conclude that $\phi>0$. Therefore, for $|D|\leq 2M-1$ there are at least $\phi>0$ paths connecting a global input vertex to a vertex in $Z$ with no vertex in $D$, contradicting the assumption that $D$ is a dominator of $Z$. ∎

A.3 Proof of Lemma 6

Proof of Lemma 6.

Without loss of generality, let us assume $|\mathcal{Y}{}_{i}^{(\mathbf{A})}|\geq|\mathcal{Y}{}_{i}^{(\mathbf{B})}|$. The proof for the case $|\mathcal{Y}{}_{i}^{(\mathbf{A})}|<|\mathcal{Y}{}_{i}^{(\mathbf{B})}|$ follows an analogous argument. Let $D{}$ be a dominator set for the set of vertices corresponding to $T_{i}^{\prime}$ with respect to $\mathcal{Y}{}_{i}^{(\mathbf{A})}$.

Consider a possible assignment to the values of $B_{i}$ such that all the values corresponding to vertices in $\mathcal{Y}{}_{i}^{(\mathbf{B})}$ are assigned value 1. Under such assignment, for every variable $a\in A_{i}$ corresponding to a vertex in $\mathcal{Y}{}_{i}^{(\mathbf{A})}$ at least one of the elementary products in $T_{i}^{\prime}$ assumes value $a$. The lemma follows combining statements (i) and (ii):

(i) There exists an assignment of the input variables of $P_{i}$ corresponding to vertices in $\mathcal{Y}{}_{i}\setminus\mathcal{Y}{}_{i}^{(\mathbf{A})}$ such that the output variables in $T_{i}^{\prime}$ assume at least $|\mathcal{R}|^{|\mathcal{Y}{}_{i}^{(\mathbf{A})}|}$ distinct values under all possible assignments of the variables corresponding to vertices in $\mathcal{Y}{}_{i}^{(\mathbf{A})}$.

(ii) Since all paths form $\mathcal{Y}{}_{i}^{(\mathbf{A})}$ to the vertices corresponding to the variables in $T_{i}^{\prime}$ intercept $D$, the values of the elementary products in $T_{i}^{\prime}$ are determined by the inputs in $\mathcal{Y}{}_{i}\setminus\mathcal{Y}{}_{i}^{(\mathbf{A})}$, which are fixed, and by the values of the vertices in $D{}$; hence the elementary products in $T_{i}^{\prime}$ can take at most $|\mathcal{R}|^{|D|}$ distinct values. ∎

A.4 Proof of Lemma 12

Before presenting the proof of Lemma 12, we present the concept of “flow of a function” which was originally introduced by Grigoriev [17] and then revised by Savage [30]. In this work, we use the following version:

Definition 3** (Grigoriev’s flow).**

A function $f:\mathcal{\mathcal{R}}^{p}\rightarrow\mathcal{\mathcal{R}}^{q}$ has a $w\left(u,v\right)$ Grigoriev’s flow if for all subsets $X_{1}$ (resp., $Y_{1}$), of its $p$ input (resp., $q$ output) variables, with $|X_{1}|\geq u$ and $|Y_{1}|\geq v$, there is a sub-function $h$ of $f$ obtained by making some assignment to variables of $f$ not in $X_{1}$ and discarding output variables not in $Y_{1}$, such that $h$ has at least $|\mathcal{R}|^{w(u,v)}$ points in the image of its domain.

Grigoriev’s flow is an inherent property of a function, agnostic to the algorithm used to compute the function. It provides a lower bound to the amount of information that suitable sub-sets of outputs encode about suitable sub-sets of inputs. Since any information about inputs that is encoded by outputs must be also be encoded by any dominator of those outputs, we have the following lower bound on the size of dominator sets.

Lemma 13** ([11, Lemma 4]).**

Let $G=(V,E)$ be a CDAG computing a function $f:\mathcal{\mathcal{R}}^{p}\rightarrow\mathcal{\mathcal{R}}^{q}$, defined on the ring $\mathcal{R}$, with $p$ input and $q$ output variables. Let $I$ (resp., $O$) denote the set of input (resp., output) vertices of $G$. Further let $I_{X}$ (resp., $O_{Y}$) denote the subset of $I$ (resp., $O$) which correspond to a subset $X$ (resp., $Y$) of the $p$ input (resp., $q$ output) variables of $f$. For any $O^{\prime}\subseteq O_{Y}$ and any $I^{\prime}\subseteq I_{X}$, any dominator of $O^{\prime}$ with respect to $I^{\prime}$ satisfies $|D|\geq w_{f_{X,Y}}(|I^{\prime}|,|O^{\prime}|)$.

A lower bound on the Grigoriev’s flow for the square matrix multiplication function $f_{n\times n}:\mathcal{R}^{2n^{2}}\rightarrow\mathcal{R}^{n^{2}}$ over the ring $\mathcal{R}$ was presented in [30, Theorem 10.5.1].

Lemma 14** (Grigoriev’s flow of $f_{n\times n}:\mathcal{R}^{2n^{2}}\rightarrow\mathcal{R}^{n^{2}}$  [30]).**

$f_{n\times n}:\mathcal{R}^{2n^{2}}\rightarrow\mathcal{R}^{n^{2}}$* has a $w_{n\times n}\left(u,v\right)$ Grigoriev’s flow, where:*

[TABLE]

We can now state the proof of Lemma 12. Recall that we denote as $G^{n\times n}$ the CDAG corresponding to the execution of a unspecified algorithm for the square matrix multiplication function with input matrices of size $n\times n$.

Proof of Lemma 12.

The statement follows by applying the results in Lemma 13 and Lemma 14 to the CDAG $G^{n\times n}$. Let $I^{\prime\prime}\subseteq I$ denote the set of all input vertices of $G^{n\times n},$ such that all paths connecting these vertices to the output vertices in $O^{\prime}$ include at least a vertex in $D{}$ (i.e., $I^{\prime\prime}$ is the largest subset of $I$ with respect to whom $D{}$ is a dominator set for $O^{\prime}$). From Lemmas 14 and 13 the following must hold:

[TABLE]

Let $I^{\prime}=I\setminus I^{\prime\prime}$. By the definition of $I^{\prime\prime}$, the vertices in $I^{\prime}$ are exactly those that are connected to vertices in $O^{\prime}$ by directed paths with no vertex in $D{}$. Since $|I|=2n^{2}$, from (17) we have $|I^{\prime}|^{2}\geq 4n^{2}\left(|O^{\prime}|-2|D{}|\right)$. ∎

Appendix B Strassen’s fast multiplication algorithm

The original version of Strassen’s fast matrix multiplication [35] is reported in Algorithm 1. We refer the reader to [37] for Winograd’s variant, which reduces the number of additions.