Tightening I/O Lower Bounds through the Hourglass Dependency Pattern
Lionel Eyraud-Dubois Guillaume Iooss Julien Langou Fabrice Rastello

























On Software Ageing Indicators in OpenStack
    
	
		Yevhen Yazvinskyi1, 
		Jasmin Bogatinovski 1, 
		Jorge Cardoso2,3, 
		Odej Kao1
	
	
	
		Technical University Berlin1, Berlin, Germany, email: {jasmin.bogatinovski, odej.kao}@tu-berlin.de
	
		2Huawei Munich Research, Munich, Germany 	
	3Departamento de Engenharia Informatica, University of Coimbra, Portugal, email: jorge.cardoso@huawei.com
		

    
=================================================================================================================================================================================================================================================================================================================================================================











When designing an algorithm, one cares about arithmetic/computational complexity, but data movement (I/O) complexity plays an increasingly important role that highly impacts performance and energy consumption.
For a given algorithm and a given I/O model, scheduling strategies such as loop tiling can reduce the required I/O down to a limit, called the I/O complexity, inherent to the algorithm itself. 

The objective of I/O complexity analysis is to compute, for a given program, its minimal I/O requirement among all valid schedules. 
We consider a sequential execution model with two memories, an infinite one, and a small one of size S on which the computations retrieve and produce data.
The I/O is the number of reads and writes between the two memories.

We identify a common “hourglass pattern” in the dependency graphs of several common linear algebra kernels.

Using the properties of this pattern, we mathematically prove tighter lower bounds on their I/O complexity, which improves the previous state-of-the-art bound by a parametric ratio.
This proof was integrated inside the IOLB automatic lower bound derivation tool.


Keywords: I/O complexity, Data Movement Lower Bound



§ INTRODUCTION



    
    We are concerned with understanding the I/O requirements 
    associated with a specific dependency pattern, that we named wide hourglass pattern and that exists in the dependency graphs 
    of some algorithms. This pattern can be found, for example, in 
    three numerical linear algebra kernels: “MGS”, “A2V” and “V2Q”.
    This paper is mostly focused on these three kernels, although the proof techniques can be adapted to any algorithm which contains
    a similar pattern in its dependency graph.
    
    The MGS kernel performs the QR factorization of an M-by-N matrix A using the Modified Gram-Schmidt (MGS) algorithm.
    MGS is a fundamental algorithm in Numerical Linear Algebra: (1) it is
    used in many applications, typically it can be used in Krylov iterative methods such as GMRES <cit.> or Arnoldi; 
    (2) it is taught to students early on in any Numerical Linear Algebra class <cit.>;
    (3) in addition, MGS appears in the popular PolyBench benchmark suite <cit.>;
    finally (4) due to its integration in Polybench but also due to its relevance for scientific computing, MGS is a standard testbed for experimental studies
    and is regularly used as a case study when demonstrating the contributions of a scientific paper in an application,
    e.g. <cit.>.
    
    
    
    
    A2V and V2Q are also related to the QR factorization of a matrix but with an approach 
    that uses unitary Householder elementary reflectors instead of orthogonal projections for MGS.
    A2V and V2Q are provided as part of the LAPACK library <cit.> as the “https://github.com/Reference-LAPACK/lapack/blob/master/SRC/dgeqr2.fGEQR2” and “https://github.com/Reference-LAPACK/lapack/blob/master/SRC/dorg2r.fORG2R” subroutines, respectively.
    “https://github.com/Reference-LAPACK/lapack/blob/master/SRC/dgeqrf.fGEQRF” and “https://github.com/Reference-LAPACK/lapack/blob/master/SRC/dorgqr.fORGQR”,
    the LAPACK optimized implementations of 
    the subroutines https://github.com/Reference-LAPACK/lapack/blob/master/SRC/dgeqr2.fGEQR2 and https://github.com/Reference-LAPACK/lapack/blob/master/SRC/dorg2r.fORG2R, are used, for example, by 
    Matlab, Numpy, R, Octave, julia, etc.
    during a command of the type, (using the Matlab syntax,)  or, if M > N, .
    Most implementations of https://github.com/Reference-LAPACK/lapack/blob/master/SRC/dgeqrf.fGEQRF and https://github.com/Reference-LAPACK/lapack/blob/master/SRC/dorgqr.fORGQR use https://github.com/Reference-LAPACK/lapack/blob/master/SRC/dgeqr2.fGEQR2 (i.e. A2V) and https://github.com/Reference-LAPACK/lapack/blob/master/SRC/dorg2r.fORG2R (i.e. V2Q) as their “panel” factorization.
    
    For these three kernels, (MGS, A2V and V2Q), 
    it is commonly accepted and  understood as common knowledge 
    that these kernels are “memory-bound” in the sequential case.
    This has been stated, for example in <cit.>.
    Julien CHECK references ; rephrase, maybe move some to abstract?
    However, no formal “proof” of this “common
    knowledge” exists.
    In addition, two recent attempts <cit.>
    have failed to provide a descriptive I/O lower bound.
    Indeed, given a cache of size S,
    the authors of these two papers provide a lower bound on the 
    amount of communication required as
    Ω( MN^2/√(S) ),
    which is a correct lower bound, but does not explain the  
    “accepted truth” that these kernels are I/O limited.
    Indeed, since the number of arithmetic operations of these three kernels
    is Ω( MN^2 ), these papers imply an upper bound
    on the arithmetic intensity in Ω(√(S)); this is correct
    but does not reflect a memory-bound kernel. 
    The best known ordering requires an amount of I/O of
    either (MN^2) or, if M<S, 
    (M^2N^2/S) <cit.>, that is, an arithmetic intensity of  Ω(max(1,S/M).
    There is a significant discrepancy between the 
    greatest known lower bound 
    and the least known upper bound.
    This piqued our interest, and the objective of this paper is to improve the state of the art. 
    
    For MGS, we show that, when the total memory footprint does not fit in cache, 
    a lower bound on the I/O complexity is Ω(MN^2 ×M/S+M).
    The proof is then generalized, and the results are extended to A2V and V2Q. These results are the main contributions of this paper.
    
    As already mentioned, in <cit.>, 
    the authors present 
    an ordering for MGS that requires an amount of I/O operations in (min(M N^2 ×M/S, M N^2 )).
    The ratio S / M represents the number of columns that fit in cache. Therefore, this tiling reduces the I/O cost by the factor “number of columns that fit in cache”: if 10 columns fit in cache, the I/O cost is divided by 10.
    We feel that this result has been overlooked by the community. As a second contribution of this paper, 
    we re-explain the tiling, generalize it to A2V and V2Q, and present  
    the first experimental performance evaluation of the tiling. We show that it can, as expected, considerably reduce cache miss and improve performance.
    In addition, 
    our Ω(M^2N^2 / (S+M)) lower bound on the I/O complexity can be written as 
    Ω(MN^2) when M>S, and as Ω(M^2N^2 / S ) when M<S.
    This proves that this tiling is optimal and that our I/O lower bound is optimal as well.
    
    
    
    
    

  
Contributions In this paper, we improve the lower bounds on the minimal data movement required by the Modified Gram-Schmidt computation.
    Then, we generalize our proof to a pattern of dependencies and apply it to derive an improved bound for both parts of the QR Householder factorization (A2V and V2Q). We also present algorithms for MGS, A2V and V2Q that match our lower bounds in terms of data movement. Finally, we present an experimental evaluation to assess the benefits of these algorithms in terms of reduction of communications and improved performance.
    
    

  
Summary of results For a matrix of size M× N and a cache of size S, we obtain the following bounds on the communication volume Q:
    
        
  * For MGS, Q≥M^2N(N-1)/4(S+M). Furthermore, if S≤ M, we also have Q≥(M-S)N(N-1)/2.
        
  * For A2V, Q ≥(3M-N)N^2(M-N)^2/24(MS+(M-N)^2). When M ≫ N, this becomes Q ≥M^2N^2/8(S+M).
        
  * For V2Q, Q ≥N (N-1)(3M-N-1)(M-N)^2/24 ( (M-N)^2 + SM ). When M ≫ N, this becomes Q ≥N(N-1)M^2/8 (M + S).
    
    
    There exist tiled orderings for each of these algorithms, based on the ideas in  <cit.>, which show that these bounds are asymptotically optimal.
    Their communication volumes are:
    
        
  * For MGS, assuming 2M < S, Q^TILED = 1/2M^2N^2/S; otherwise Q = 2MN^2; 
        
  * For A2V and V2Q, assuming 2M < S, Q^TILED = 1/2M^2N^2 - N^3/3/S;
        otherwise Q = 2MN^2 - 2N^3/3.
    





When designing an algorithm, we usually reason about its computational complexity, to estimate the increase in its execution time, when its problem sizes increase.
However, the amount of computation is not the only factor when estimating the performance of an algorithm.
The data movement (I/O) also plays an increasingly important role in both performance and the energy consumed by an algorithm.



In order to minimize the amount of data movement of an algorithm, we modify its schedule by using program transformations, such as the tiling transformation <cit.>.
However, finding an optimal schedule is not trivial, due to the size and complexity of the optimization space.
The notion of I/O complexity of an algorithm indicates the minimal amount of data movement required for any valid schedule of an algorithm.
It provides an indication on how far it is possible to optimize the I/O of an algorithm.





However,  due to the number of possible valid schedules, it is not possible to directly compute the I/O complexity.
Instead, we search for a lower and upper bound on the minimal amount of data movement.
To find an upper bound, one only needs to exhibit a valid schedule and compute its I/O cost.
However, finding a lower bound requires us to reason over all possible schedules, which requires a mathematical proof based on the data dependency graph of the program.



Several proof techniques exist for deriving lower bounds, such as the wavefront <cit.> or the K-partitioning technique <cit.>.
Depending on the shape of the dependencies, some of these techniques might work better than the others: for example, the wavefront technique is usually the most effective for stencil-like computation, while the K-partitioning technique is usually better for linear-algebraic computation.
Other algorithms might require a specialized proof to obtain an asymptotically tight bound (e.g., for the SYRK kernel <cit.>).


However, there are still some kernels for which the ratio between their proven lower bound and their best upper bound is parametric. For example, such a ratio is described in Table 1 of  <cit.> for the kernels of the Polybench benchmark suite <cit.>.


Some of the proof techniques mentioned above have been automatized and implemented in automatic data movement lower bound derivation tools, such as IOLB <cit.>, so that they can be applied to any kernel given as an input.




  
Contributions In this paper, we consider several important linear algebra kernels whose dependence graph exhibits a dependency pattern called the hourglass pattern.
We propose a new proof technique that uses the properties of this pattern to improve the lower bound on the minimal data movement required by these kernels.

In more detail:

    
  * We define the hourglass pattern, a pattern of the dependencies of a program, and present its properties.

    
  * We provide a lower bound derivation proof, based on an adaptation of the K-partitioning technique, that tighten the lower bound of a program exhibiting an hourglass pattern.
    This proof has been fully automatized, inside the tool IOLB <cit.>.
    
    
  * We present new data movement lower bounds for several linear algebra kernels: Modified GramSchmidt (MGS), QR Householder (GEQR2 and ORG2R in the LAPACK <cit.> library; also called A2V and V2Q), bidiagonal matrix reduction (GEBD2) and Hessenberg matrix reduction (GEHD2).
    
    For all of these kernels, their asymptotic bound was improved by a parametric factor, compared to the bound obtainable by the classical K-partitioning technique.
    Figure <ref> and Figure <ref> summarize the new lower bounds found for several linear algebra kernels.

    
  * We also provide tiled orderings for MGS and Householder, resulting in upper bounds that asymptotically match these new lower bounds.
    This proves the optimality of the new I/O lower bounds.







  
Outline
In Section <ref>, we provide some background on the I/O complexity and introduce the K-partitioning method for deriving a lower bound on the minimal amount of data movement of a computation.
In Section <ref>, we present the hourglass pattern, a pattern of dependencies whose properties can be used to improve the derived lower bound.
In Section <ref>, we show how to exploit the hourglass pattern to adapt the K-partitioning method, in order to obtain a tighter bound.
In Section <ref>, we list different linear algebra kernels exhibiting an hourglass pattern, and their associated improved lower bound.
Additionally, Annex <ref> contains the tiled algorithm for two of these kernels. The data movement of these algorithms provides an upper bound to the minimal amount of data movement, which matches asymptotically their new lower bound.








§ BACKGROUND - I/O COMPLEXITY AND THE K-PARTITIONING METHOD




In this section, we present a state-of-the-art method of proof – the K-partitioning method – which infers a lower bound on the amount of data movement needed by a computation.




  
Memory model and I/O complexity
We consider a simple two-level memory model, composed of (1) a slow memory of unbounded size, and (2) a fast memory of size S.
Both memories can transmit data from one to another, as long as the constraint on the size of the small memory is satisfied.
When we perform an operation, the data used must be present in the small memory, and the data produced must be committed in the small memory.


We consider a program, which performs a collection of operations organized in statements.
Each statement SX has multiple instances SX[i⃗], where i⃗ is a vector of the surrounding loop indexes.


In this paper, we consider a subclass of programs called polyhedron (or affine) programs.
These programs are combinations of nested loops and statements such as:
(i) the loop bounds are affine constraints using the surrounding loop indexes and the program parameters (e.g., the sizes of an input array);
(ii) the array accesses are affine expression of the surrounding loop indexes and program parameters. All the programs presented in this paper, such as Figure <ref>, are polyhedral programs.
Also, we will call a quantity “parametric” when it is a function of the parameters of the program, which are considered as symbolic constants.


A statement instance might depend on the data produced by another instance, for example when the first instance uses a value produced by the second instance.
We call this a dependency between these two statement instances.
These dependencies impose constraints on the order of execution of the program.
This order of execution is called the schedule or the ordering.
Frequently, the data consumed/produced by the statements of a program are too big to fit all at once in the small memory, it it thus necessary to spill, i.e., to transfer some data back into the slow memory and retrieve it later when needed.
However, doing so increases the amount of data movement between both memories.


Given a valid schedule for a program, the I/O cost for this program and for this ordering is the amount of data movement required, i.e. the number of data transfers between both memories.
The I/O complexity of a program is the minimal I/O cost that can be reached by any valid schedule.


This quantity is interesting, in particular in the context of an architecture where the transfer of data is the limiting factor for performance.
Knowing the minimal amount of data transfers is a good algorithmic indicator to know if it could be theoretically optimized further.


However, because we need to find the minimal I/O cost for all possible orderings, the exact I/O complexity is hard to evaluate.
Instead, we rely on bounds on the I/O complexity of a program: we can provide a mathematical proof for the lower bound, and exhibit an ordering (i.e., an implementation of a program) that reaches an I/O cost and provides an upper bound.




  
CDAG and red-white pebble game
When trying to prove a lower bound, one should decide whether redundant computation is allowed or not. 
The red-white pebble game, a variation of the red-blue pebble game of Hong and Kung <cit.>, was introduced by Olivry et al. <cit.> in order to model the state of the memories during the execution of a program without recomputation, which matched the assumptions we make in this paper.


This game is played on the Computational Directed Acyclic Graph (CDAG) of the program. This is a directed graph G, where

    
  * the nodes V represent the computation (statement instance) of a program, and
    
  * the edges represent the flow dependencies between the computations of the program.

Notice that the inputs of a program are nodes that do not have incoming edges. The outputs of a program are a subset of nodes O ⊂ V; they might have outgoing edges.


During a red-white pebble game, red and white pebbles are placed on the nodes of a CDAG.
A white pebble represents a computation that was performed, and a red pebble represents a computation whose output is currently stored in the small memory.
A game follows this set of rules:

    
    
  * At the start, the only pebbles in the CDAG are white pebbles, placed on the inputs of the program.

    
    
  * At most S red pebbles can be simultaneously present on the nodes of the CDAG.

    
    
  * Spill: a red pebble can be removed from a node.

    
  * Compute: When a node does not have a white pebble, but all its predecessors have red pebbles, then we can place both a white and a red pebble on it.

    
  * Load: A red pebble can be added to nodes with a white pebble.

    
    
  * The game ends when each node has a white pebble on it.

Notice that once a white pebble is placed on a node, it cannot be removed. This prevents recomputation.
In order to compute the amount of data movements, we focus on the number of red pebbles added with the Load rule during a game.

This means that we only focus on the “Load” portion of the data movements and ignore its “Store” part.
The resulting bounds are still valid, and because the number of “Load” often dominates the number of “Store”, their tightness should not be strongly impacted.
This assumption is identical to the one made in <cit.>.




  
K-partitioning method The K-partitioning method introduced in the seminal paper of Hong and Kunk <cit.> is a proof technique that allows to derive a lower bound.

The first idea is to consider a partition of the CDAG and games that play on each set of the partition one by one.

Then, we consider the notion of K-bounded set.
An inset of a set E of nodes of the CDAG, noted InSet(E), is the set of data used by E but not produced by a computation of E.
A K-bounded set is a set of nodes E of the CDAG whose inset has a size at most K: |InSet(E)| ≤ K.

This notion is interesting because an (S+T)-bounded set requires at least T additional data to fit in the small memory (of size S). Thus, even if the small memory is filled with interesting data, we will need at least T load operations to perform the computations of this set.


In addition, we assume that our K-bounded sets are convex: if there is a dependency chain between two points of a K-bounded set E, then all the intermediate points must belong to E.

Finally, a K-partition is a partition into convex K-bounded sets.
The idea is to consider all K-partitions of a CDAG and to count how many sets are in this partition.
By choosing K=S+T, we know that there will be at least T loads per set of the partition.
Therefore, a lower bound on the number of loads is T times the minimal number of sets in such a partition.


Let S be the size of the small memory, and for any T>0 let U be the maximal size of a (S+T)-partition. Let V be the set of nodes of the CDAG of the program.


Then, a lower bound on the number Q of data movement of the program is:

    T.⌊|V|/U⌋≤ Q






Then, we pick a value of T (which means a value of K) that leads to the tightest lower bound.

To estimate the minimal number of sets in a K-partition, we can estimate the maximum size of a set inside this partition.
In other words, an upper bound on the size of a K-bounded set can be transformed into a lower bound on the amount of data movement required.



  
Upper bound on the size of a K-bounded set
An upper bound on the size of a K-bounded set E can be obtained by analyzing the dependencies of the program.

Indeed, for a polyhedral program, dependencies between its statement instances are associated with affine relations, matching the loop indices of the data producing instance with the data consuming instance.
When examining the path of affine dependencies starting from any node of E to a node of the inset of E, we can either obtain a projection or a translation. In both cases, the image of E through these affine functions ϕ can be mapped to geometrical borders or projections of E, and can be associated with parts of InSet(E).
This is the key geometrical intuition that leads us to use the Brascamp-Lieb theorem.

The Brascamp-Lieb theorem is a geometrical way to bound the volume of a set by the volume of its projections, which can be bounded by K.

Let d and d_j be non-negative integers and ϕ_j: ℤ^d ↦ℤ^d_j be a collection of group homomorphisms for all 1 ≤ j ≤ m.

If we have a collection of coefficients s_j ∈ [0,1] such that, for any subgroup ℋ⊂ℤ^d:

    rank(ℋ) ≤∑_j=1^m s_j × rank(ϕ_j(ℋ)).


Then, for any non-empty finite set E ⊂ℤ^d:

    |E| ≤∏_j=1^m |ϕ_j(E)|^s_j.







For example, if we consider a 3D set and the 3 canonical projections ϕ_j,k(i,j,k)=(j,k), ϕ_i,k(i,j,k)=(i,k) and ϕ_i,j(i,j,k)=(i,j), this theorem gives us the following inequality between the volume of E and the area of its faces ϕ_x(E):


    |E| ≤ |ϕ_j,k(E)|^1/2× |ϕ_i,k(E)|^1/2× |ϕ_i,j(E)|^1/2.


As another example, we could consider instead the projections ϕ_i(i,j,k) = (i), ϕ_j(i,j,k) = (j) and ϕ_k(i,j,k) = (k), to obtain the following inequality:

    |E| ≤ |ϕ_i(E)| × |ϕ_j(E)| × |ϕ_k(E)|.



In the case of a K-bounded set, we consider the path of dependencies to automatically derive these projections ϕ∈Φ.
So, ϕ(E) can be mapped to one part of the inset of E, and its size is bounded by K.






§ THE HOURGLASS PATTERN




In this section, we describe the intuition of our core contribution.
We consider a specific pattern of dependencies, called the hourglass pattern, that forces a convex K-bounded set to have a specific shape.
We can exploit this property to significantly improve the derived I/O complexity lower bound of programs that exhibit such a pattern.

In the whole section, we use the Modified Gram-Schmidt algorithm as an illustrative example, whose right-looking variant is provided in Figure <ref>.
Using the automatic tool IOLB <cit.> to apply the K-partitioning method (described in Section <ref>) to the MGS computation results in a lower bound in Ω(MN^2/√(S)).


By using the hourglass pattern, we obtain a more precise lower bound:

    M^2N(N-1)/8(S+M)≤ Q(MGS)














 §.§ Intuition of the hourglass pattern










  
Intuition Figure <ref> presents the main idea of the hourglass pattern.
It is a repeating succession of reduction and broadcast statements, such that the number of elements reduced/broadcasted is parametric, thus greater than the cache size S.
There are 3 categories of dimensions in this pattern:
(a) the dimensions over which the reduction and the broadcast are performed (horizontal axis of Figure <ref>),
(b) the “temporal” dimensions over which the hourglass pattern is repeated (vertical axis of Figure <ref>), and
(c) the neutral dimensions that do not interact with the hourglass pattern.





  
Running example The hourglass pattern appears on several linear algebra kernels, including the Modified Gram-Schmidt kernel (Figure <ref>).
The pattern appears between the last two statements: statement SR which updates R[k][j], and statement SU which updates A[i][j].
The statement SR is a reduction along the i dimension and uses, in particular, all the values of A[·][j] produced during the previous iteration of k.
The statement SU broadcasts R[k][j] across the i dimension to update all the A[·][j] of the current iteration of k.
Therefore, dimension k is a temporal dimension, dimension i is the reduction/broadcast dimension and dimension j is a neutral dimension. There is exactly one dimension in each category in this example, but in general, there might be several.




  
Consequences of the hourglass pattern
When considering a K-bounded set over this pattern, we notice that if it spans over several iterations of the temporal dimension t, then the set must include all the nodes of the broadcast/reduction in between, due to the convexity property of the set.
Therefore, we have two situations:

    
  * Either the K-bounded set spans over several iterations of t and includes all the nodes over the reduction/broadcast dimension. Notice that this is not always possible, depending on the size of the broadcast/reduction dimension and the value of K.
    
  * Or, the K-bounded set is “flat” along the t dimension.

For both situations, we can deduce much stronger constraints on the sizes of the projection |ϕ(E)| used in the Brascamp-Lieb theorem (Theorem <ref>).
This provides the derivation of an improved lower bound compared to the classical methodology.






 §.§ Hourglass pattern - formal definition








In this section, we provide a formal definition of the hourglass pattern.



  
Preliminary notations

Given an instance for a statement SX, we call an iteration vector the tuple of the (integral) values of its surrounding loop indices.
The iteration domain 𝒟_SX of the statement SX is the set of iteration vectors that respect the conditions on the indices of the surrounding loops.


As mentioned above, in general, there might be several reduction or temporal dimensions. Thus, we will consider sets of dimensions, and represent their iteration as a vector i⃗ = (i_1, i_2, …).

Given an iteration k, we write  k+1 to represent the next valid lexicographic value of k. We extend this notation to k+n where n is an integer.



  
The hourglass pattern
Considering a statement S of the CDAG of a program, the hourglass pattern is a pattern of dependencies with the following properties:

    
    
  * Partitioning of the dimensions. The dimensions of the statement SX can be partitioned into 3 groups: (i) the temporal dimensions k⃗, (ii) the reduction/broadcast dimensions i⃗, and (iii) the neutral dimensions j⃗.
    For simplicity of the presentation, we assume that k⃗ are the first/outer dimensions and that i⃗ are the last/inner dimensions.
    
    
    
  * Path in the dependence graph. For any valid value of i⃗ and i⃗'⃗, there is a dependency chain between the instances SX[k⃗,j⃗,i⃗] and SX[k+1, j⃗, i⃗'⃗].
    
    
    
  * Large width of the hourglass. Let us consider W, the number of statement instances on all the dependency chains between SX[k⃗,j⃗, i⃗] and SX[k+2,j⃗, i⃗]. This expression depends on the parameters of the program, and cannot be bounded by a constant value.



Notice that, to have such a chain of dependencies, it must include a reduction and a broadcast.
The dimensions i⃗ are the dimensions which are reduced on and broadcasted over.
The dimensions k⃗ are the dimensions that are incremented by a constant factor when looping along this loop.
So, once a path is found, partitioning the dimensions should be unambiguous, if this condition is also satisfied.

The automatic detection of such an hourglass pattern has been implemented inside the IOLB tool, using a polyhedral library <cit.>.





  
Examples
For the MGS computation (Figure <ref>), we consider the statement SU and the cycle of dependencies going through the SR.
We confirm that there are 2M statement instances inside a dependency chain between two instances SU[k,j,i] and SU[k+2,j,i]: SR[k+1,j,·] and SU[k+1,j,·].

The same reasoning would hold if we considered the statement SR instead of SU.

For the A2V QR Householder computation (Figure <ref>), we consider the statement SU, and the cycle of dependencies going through the SR statement.

There are (M-k) statement instances SR[k,j,·], inside a dependency chain between two instances SU[k,j,i] and SU[k+1,j,i].








§ LOWER BOUND PROOF USING THE HOURGLASS PATTERN




In this section, we show how to exploit the hourglass pattern of a program to derive a tighter lower bound on the data movement.
We will use the Modified Gram-Schmidt program (Figure <ref>) as a running example to illustrate our proof.

Once again, this proof has been integrated inside the automatic data movement lower bound derivation tool, IOLB <cit.>, and is applied when an hourglass pattern is detected.



  
Preliminary notations
Let us denote with E a set of integral iteration vectors. The cardinality of such a set |E| is the number of integer points inside this set.


Given a dimension k, ϕ_k : (i,j,k) ↦ (k) is the projection to the dimension k. This notation can be extended to several dimensions. For example, ϕ_j,k: (i,j,k) ↦ (j,k).


Given a set E, a dimension k and a value k_0, E_k=k_0 is a slice of E, i.e., the set of points of E whose value along the dimension k is k_0: E_k=k_0 = { (i,j,k) ∈ E  |  k=k_0 }.
We extend this notation to several dimensions, as in E_k=k_0,j=j_0.
In the rest of the paper, we use a compact notation for slices, by omitting the dimension when it is not ambiguous, like E_k_0,j_0.



  
Starting point and intuition
The goal of the proof is to find an upper bound on the size of a K-bounded set, to transform it into a lower bound on the data volume required by a portion of the program.
So, we start the proof with an arbitrary convex K-bounded set E containing instances of the broadcast statement SX[k⃗, j⃗, i⃗], which we assume is part of an hourglass pattern. Just like in the classical proof, an analysis of the dependencies of the program yields a set Φ of projections ϕ_x⃗ for some dimensions x⃗, which project on the inset of E.


The classical proof involves applying the Brascamp-Lieb theorem on the K-bounded set E, by bounding the size |ϕ_x⃗(E)| of each of these projections by K.

In our proof, we split E into two parts (Section <ref>) and we adapt the set of projections when we apply the Brascamp-Lieb theorem on each of these parts (Sections <ref> and <ref>) to obtain more precise bounds.



  
Running example In the case of MGS, by following the chain of dependencies associated with the accesses A[i][j], Q[i][k] and R[k][j] of statement SU, we infer that the projections are ϕ_i,j, ϕ_i,k and ϕ_k,j.

When following the classical proof, the application of the Brascamp-Lieb theorem results in the following inequality:

    |E| ≤ |ϕ_i,j(E)|^1/2 . |ϕ_i,k(E)|^1/2 . |ϕ_k,j(E)|^1/2≤ K^3/2

We prove in Section <ref> a tighter upper bound: |E| ≤K^2/M + 2K.





 §.§ Part 1 - Decomposition of the K-bounded set




In the first part of the proof, we decompose E into the union of two fragments: (i) I' which has volume along the k⃗ dimensions, and (ii) F which is flat along the k⃗ dimensions.
Both parts have their own upper bound, obtained with different reasoning, that is described in Section <ref> and Section <ref>.





  
Decomposition of E
We consider the number of different values of k⃗ for a given value of j⃗:

    Tick_j⃗ = {k⃗|∃i⃗,   SX[k⃗,j⃗,i⃗] ∈ E_j⃗}.

We split E into two sets:

    
  * E' = ∪_j⃗∈ J_3+ E_j where J_3+={j⃗| 3≤ Tick_j⃗}, and
    
  * E” = ∪_j⃗∈ J_12 E_j where J_12={j⃗| 1 ≤ Tick_j⃗≤ 2 }.


The intuition between this separation is that (i) E' contains the connected components which need to include an entire line of statement instances along the i⃗ dimensions, and (ii) E” contains the connected components which are “flat” along the k⃗ dimensions.

We start by focusing on E', to show that these components need to include a parametric number of iterations along the i⃗ dimension, using the hourglass pattern.



    Given some j⃗∈ϕ_j⃗(E'), let us consider the slice E'_j⃗ along the dimensions j⃗.
    Let us consider:
    
        
  * (j⃗) = min_k⃗{k⃗| SX[k⃗, j⃗, i⃗] ∈ E'_j⃗},
        
  * (j⃗) = max_k⃗{k⃗| SX[k⃗, j⃗, i⃗] ∈ E'_j⃗},
        
  * a⃗, any subset of the indices of i⃗,
        
  * |ϕ_a⃗(𝒟_S)| ≥ W_a⃗,
        a lower bound on the size of the projection of 𝒟_S, the iteration domain of the statement SX.
        
    

    Then:
    
        
  * E'_j⃗ is a connected component.
        
  * For all (j⃗) <_lexk⃗ <_lex(j⃗), |ϕ_a⃗(E'_j⃗, k⃗)| ≥ W_a⃗.
    
    




    (1) Because SX satisfies the hourglass pattern properties, for any two instances SX[k⃗,j⃗, i⃗] and SX[k⃗'⃗,j⃗, i⃗'⃗] in E'_j⃗ where k⃗<_lexk⃗'⃗, we can prove that there is a chain of dependencies from one statement instance to another.
    Therefore, using the convexity property of E, we conclude that E'_j is a connected component.
    
    
    (2) By considering one instance of index k⃗ = and another of index k⃗ =, we can show that all the SX[i⃗, j⃗, k⃗] are in the middle of a dependency chain between the two of them.
    Therefore, by projecting the i⃗ on the subset of dimension a⃗, we conclude:
    
    ∀(j⃗) <_lexk⃗ <_lex(j⃗),  |ϕ_a⃗(E'_j⃗, k⃗)| ≥ |ϕ_a⃗(𝒟_S)| ≥ W_a⃗.




We consider E' = I' ⊎ B' defined by:

    
  * I' = ∪_j⃗∈ J_3+ (E_j⃗ - E_j⃗,(j⃗) - E_j⃗,(j⃗)): the “inside” of the connected components of E', according to dimensions k⃗.
    
  * B' = ∪_j⃗∈ J_3+ (E_j⃗,(j⃗)∪ E_j⃗,(j⃗)): the “boundaries” of the connected components of E', according to dimensions k⃗.


We define F = B' ⊎ E” the flat parts of E, and we adapt our previous decomposition of E to obtain the desired decomposition:

    E = E' ⊎ E” = I' ⊎ F.





 §.§ Part 2 - Bound on the size of I'




Let us focus on the upper bound of the size of I', using Lemma <ref>.
The goal is to have tighter bounds on the projections involving some of the reduce/broadcast dimensions i⃗, to use them with the Brascamp-Lieb theorem, instead of the classical “≤ K” bound.





Let us consider a projection ϕ_x⃗, a⃗∈Φ, where a⃗ is a subset of i⃗, and x⃗ is a subset of (j⃗, k⃗).

Assume that we have a lower bound |ϕ_a⃗(𝒟_S)| ≥ W_a⃗ on the size of the projection of the iteration domain 𝒟_S of statement S.
Then:

    |ϕ_x⃗(I')| ≤K/W_a⃗.






    Because I' is a subset of a K-partition, then we have |ϕ_a⃗,x⃗(I')| ≤ K.
    
    By slicing I' along the dimensions of x⃗, we also have:
    
    |ϕ_a⃗,x⃗(I')| = |ϕ_a⃗,x⃗( ∪_x” I'_x” )| 
                = ∑_x” |ϕ_a⃗,x⃗(I'_x”)|.

    
    Furthermore:
    
    |ϕ_a⃗,x⃗(I'_x”)|    ≥ |ϕ_a⃗ I'_x”)|       (a slice along x⃗ only has a single
             point to be projected along x⃗)
       ≥ |ϕ_a⃗(I'_j”, k”)|         (I'_j”, k”⊂ I'_x”, for some (j”,k”)
             matching x” on its dimensions)
       ≥ W_a⃗.       (by Lemma <ref>)

        
    Therefore: |ϕ_a⃗,x⃗(I')| ≥ W_a⃗× |ϕ_x⃗(I')|.
    In other words:
    
    |ϕ_x⃗(I')| ≤|ϕ_a⃗,x⃗(I')|/W_a⃗≤K/W_a⃗.




To obtain a bound on |I'|, we use the set of projections Φ, modified in the following way:

    
  * We add a projection on i⃗ whose bound is:
    |ϕ_i⃗(I')| ≤ W, where W is the width of the hourglass.
    
  * When a projection on (x⃗, a⃗) shares some of its dimensions a⃗ with i⃗, we use the projection on x⃗ instead, and the bound given by Lemma <ref>.
    
  * The rest of the projections, which do not involve dimensions of i⃗, are unchanged and associated with a classical upper bound |ϕ_x⃗(I')| ≤ K.

Then, we apply the Brascamp-Lieb theorem, while optimizing the values of the power s of the size of the projections |ϕ(E)|.
Notice that some projections can be filtered out through this optimization process (s=0), to obtain a tighter bound on I'.






  
Running example

We notice that Φ contains two projections ϕ_i,j and ϕ_i,k,
both of them using the dimension i and another dimension.
Therefore, from Lemma <ref>, we have the following bounds:

    |ϕ_j(I')| ≤K/M  and   |ϕ_k(I')| ≤K/M.


We apply the Brascamp-Lieb theorem on I', using the projections on i, on j instead of a projection on (i,j), and on k instead of a projection on (i,k), each with coefficient 1:

    |I'| ≤ |ϕ_i(I')| × |ϕ_j(I')| × |ϕ_k(I')|.

By using our specialized bounds, we have:

    |I'| ≤ M ×K/M×K/M = K^2/M.









 §.§ Part 3 - Bound on the size of F




We now focus on the bound of the remaining part of E, i.e., F = (B' ⊎ E”). 
Our goal is to exploit the fact that this set is “flat” along the k⃗ dimensions, to improve the upper bounds on the projections used in the Brascamp-Lieb application to F.

In particular, our proof will consider each F_j⃗ separately, the slice of F for a given value j⃗. Then, by picking a list of well-chosen projections for the Brascamp-Lieb theorem, we obtain interesting bounds on the size of F_j⃗, that are finally summed to obtain our bound on F.


We recall that E' and E” contains the j⃗-slices of E for the values of j⃗ such that Tick_j⃗ are respectively at least 3 and at most 2. Since B' ⊂ E', B' and E” do not share the same set of j⃗.
So we get a “flatness bound”: ∀j⃗∈ϕ_j⃗(F),   |ϕ_k⃗(F_j⃗)| ≤ 2.

In addition, we notice that because F_j⃗⊂ E, for any projection ϕ∈Φ, |ϕ(F_j⃗)| ≤ |ϕ(E)| ≤ K.






In this part of the proof, instead of focusing on F, we apply the Brascamp-Lieb theorem to the slice F_j⃗.

As in Section <ref>, we start with the same list of projections Φ to the in-set of E, obtained by inspecting the chain of dependencies of the considered statement.

Then, we customize this list of projections to exploit the properties of F_j⃗:

    
  * We add a projection on k⃗ and the flatness bound |ϕ_k⃗(F_j⃗)| ≤ 2.
    
  * We identify a projection that involves a non-empty subset of the dimensions of j⃗: after applying the Brascamp-Lieb theorem, we will not try to immediately use an upper bound for this projection. Let us call ϕ_w⃗ this projection.
    
  * The remaining projections are left alone, and associated with the classical upper bound |ϕ_x⃗(F_j⃗)| ≤ |ϕ_x⃗(F)| ≤ K.


At that point, we have obtained an upper bound on |F_j⃗| of the following form:

    |F_j⃗| ≤ e × |ϕ_w⃗( F_j⃗ )| .

where e is a parametric expression using the parameter K and independent of the value of j⃗.
Notice that the Brascamp-Lieb power above |ϕ_w⃗| must be equal to 1, due to the fact that this is the only projection involving the j⃗ dimensions.






Let us consider the collection of projected sets ϕ_w⃗(F_j⃗).
Two of these projected sets are either (i) identical, along the dimensions of j⃗ which are not present in w⃗, or (ii) disjoint.
Because the distinct values of these projections are always a subset of the inset of E, of size K, the sum of the size of the disjoint union of these distinct values is also bounded by K.
Let us call R the number of values that can be taken by the dimensions of j⃗ which are not in w⃗. This is also the maximum number of times a ϕ_w⃗(E) projects on the same value.
Notice that R=1 when w⃗ covers all the dimensions of j⃗, and that R can be a parametric expression in general.






So, we have:

    ∑_j⃗∈ϕ_j⃗(F) |ϕ_w⃗( F_j⃗ )| ≤ R × K.



Combining all these observations, we obtain the following upper bound on the size of F:

    |F|     = ∑_j⃗∈ϕ_j⃗(F) |F_j⃗|   ≤∑_j⃗∈ϕ_j⃗(F) e × |ϕ_w⃗( F_j⃗ )|  
       ≤ e ×∑_j⃗∈ϕ_j⃗(F) |ϕ_w⃗( F_j⃗ )| ≤ e × R × K .







  
Running example
In the MGS case, the flatness bound yields: ∀ j, |ϕ_k(F_j)| ≤ 2.
We can apply the Brascamp-Lieb theorem on F_j, with the projection ϕ_i,j, together with the projection ϕ_k:

    |F_j| ≤ |ϕ_k(F_j)| × |ϕ_i,j(F_j)| ≤ 2 × |ϕ_i,j(F_j)|.


And since dimension j is included in the dimensions (i,j) of the projection, we have R=1 and therefore:

    |F| = ∑_j ∈ϕ_j(F) |F_j| ≤ 2 ×∑_j ∈ϕ_j(F) |ϕ_i,j(F_j)| ≤ 2K.






 §.§ Part 4 - Wrapping things up



We have found in Section <ref> and Section <ref> an upper bound of the size of the two parts of E.
We simply sum them together to obtain an upper bound of the size of E, then apply Theorem <ref> to deduce a lower bound on the data movement.




  
Running example
Thanks to the above results, we can obtain lower bounds on the data movement of MGS.

    The communication volume Q for the MGS algorithm on a M× N matrix can be bounded as follows:
    
    M^2N(N-1)/8(S+M)≤ Q


    Furthermore, if S≤ M, we also have: 
    
    (M-S)N(N-1)/4≤ Q

    



From the previous results, we have:

    |E| = |I'| + |F| ≤K^2/M + 2K .


Then, by using Theorem <ref> with K=2S:

    (K-S) ×MN(N-1)/2·(K^2/M + 2K) = M^2N(N-1)/8(S+M)≤ Q



Due to Lemma <ref> and because we have at least an input dependency involving the dimension i, |InSet(E')| > M.
So, if we have S ≤ M, then E' must be empty. Therefore, E = F, and we can use only the second part of the bound: |E|≤ 2K.

Using Theorem <ref> again, but this time with K = M, we obtain:

    (K-S)·MN(N-1)/2 × 2K = (M-S)·N(N-1)/4≤ Q











§ EXPERIMENTAL RESULTS - NEW LOWER BOUNDS




In this section, we report the data movement lower bounds generated by IOLB for four kernels exhibiting an hourglass pattern. We compare the results using our technique (new bound) with those obtained without it (old bound).
These kernels are:

    
  * Modified Gram-Schmidt (Figure <ref>), already used as a running example in Section <ref>.
    
  * QR Householder algorithm: both its A2V (Figure <ref>) and V2Q parts (left-looking variants of respectively the GEQR2 and ORG2R subroutines in LAPACK <cit.>).
    
  * Reduction to a bidiagonal matrix (GEBD2 subroutine)
    
  * Reduction to a Hessenberg matrix (GEHD2 subroutine)
    
    






Figure <ref> summarizes all the newly found full lower bounds, and Figure <ref> focuses on their leading term, to emphasize their improvements.
More precisely:

    
  * Section <ref> contains an asymptotic analysis of the MGS lower bound.
    
  * Section <ref> presents the lower bound for both QR Householder parts, and the GEBD2 computation.
    
  * Section <ref> presents the lower bounds for the GEHD2 computation.
    



Annex <ref> contains a description of a tiled algorithm for MGS and for HH A2V and the computation of their amount of data movement.
This provides an upper bound to the minimal amount of data movement required by these algorithms, which matches asymptotically the provided lower bound.









































 §.§ MGS - Asymptotic analysis




In this section, we analyze the bound obtained in Theorem <ref> for MGS, by specializing it for different ordering of M and S:

    
  * If S≤ M/2, we have M/2≤ M-S, so that the second bound yields:
    
    MN^2/8 = Ω(MN^2)≤ Q

    
  * If M/2≤ S, we have S+M≤ 3S, so that the first bound becomes:
    
    M^2N^2/24S = Ω(M^2N^2/S)≤ Q




The first result matches the amount of data movements obtained by the classical ordering of the MGS algorithm, as presented at the end of Section <ref>. Both the algorithm and the bound are thus asymptotically optimal when S is small.

Demmel et al. <cit.> propose a tiled ordering for the case 2M≤ S (in Section F.2 of their paper), that achieves an amount of data movement O(M^2N^2/S), thus matching asymptotically the second upper bound. Both of these bounds are thus optimal up to a constant factor. For reference, the tiled ordering is detailed in Appendix <ref>, together with the proof on its amount of data movement.


We can compare these bounds with the lower bound returned by the classical hourglass-less proof, whose asymptotic bound is Ω(MN^2/√(S)). Our first bound for small values of S is stronger by a factor of Theta(√(S)).
By writing our second bound as Ω(M/√(S)·MN^2/√(S)), we see that our bound is stronger by a factor of Θ(M/√(S)). Since the input matrix has size M × N, we can assume that S< MN, otherwise, the whole matrix fits in the cache and there is no need to minimize data transfers. Besides, we know that N≤ M, which leads to S<M^2.
We can conclude that M/√(S) > 1: our bound is asymptotically at least as strong as the previous bound.

The results of Theorem <ref> can also be presented differently, to improve the constant in front of the dominant term. Indeed:

    If S≪ M, the first bound yields   MN^2/4   ≤ Q,
       Similarly, if M≪ S, the second result becomes   M^2N^2/8S   ≤ Q.








 §.§ Householder QR factorization and GEBD2




In this section, we present the bounds to both parts of the Householder QR factorization (LAPACK routines GEQR2 et ORG2R). The computation of the first part (A2V) is given in Figure <ref> and the computation of the second part (V2Q) in Figure <ref>.








By using the hourglass reasoning, we obtain the following new lower bounds for both kernels.



    The communication volume Q for the A2V part of the HH algorithm on a M × N matrix, with M>N, can be bounded as follows:
    
    (3M-N)N^2(M-N)^2/24(MS+(M-N)^2)≤ Q

    If M ≫ N, then the bound becomes:
    
    M^2 N(N-1)/8(S+M)≤ Q






    The communication volume Q for the HH V2Q algorithm on a M × N algorithm, with M > N, can be bounded as follows:
    
    N (N-1)(3M-N-1)(M-N)^2/24 ( (M-N)^2 + SM )≤ Q

    When M ≫ N, this bound becomes:
    
    N(N-1)M^2/8 (S+M)≤ Q







One interesting detail of the proof concerns the detection of the hourglass, and the criteria on its size (third criterion introduced in Section <ref>).
For the MGS computation, the size of its hourglass was constant and equal to M.
In the case of both Householder computations, the size of their hourglass is parametrized by the outer loop iteration value k, and is equal to (M-1-k).
Its smallest value happens for k=N-1, and is (M-N).
By using this value in Lemma <ref>, we can derive the announced lower bound.


The lower bound proof of the GEBD2 subroutine is similar to both Householder proofs.



    The communication volume Q for the GEBD2 subroutine on a M × N matrix, with M≥ N, can be bounded as follows:
    
    MN^2(M-N+1)/8(S+M-N+1)≤ Q

    
    If M ≫ N, then the bound becomes:
    
    M^2 N^2/8(S+M)≤ Q









 §.§ Hessenberg matrix reduction







In this section, we present the bound for the Hessenberg matrix factorization kernel (GEHD2 kernel in LAPACK), whose code is in Figure <ref>.



By using the hourglass reasoning, we obtain the following new lower bound for the GEHD2 kernel.


    The communication volume Q for the Hessenberg matrix factorization algorithm on a N × N algorithm can be bounded as follows:
    
    1/12 . N^4/N+2S≤ Q

    
    When N ≫ S,
    
    N^3/24≤ Q






As in Section <ref>, the size of the hourglass depends on the value of iteration j of the outermost loop (temporal dimension), and is equal to (N-2-j), where 0 ≤ j < N-2.
This means that its minimal value is 1, which causes an issue with our proof.

A way to solve this problem is to introduce a new parameter M < N-2 and to consider a loop splitting of the outermost temporal dimension j at iteration M.
First, notice that a loop splitting does not change the dependencies of a program, thus the lower bound of the split version will apply to the original one.

The first half of the split program satisfies the hourglass pattern and has a minimum size of its hourglass of (N-M-1). Thus, as long as this quantity is not bounded by a constant, the hourglass reasoning should apply.
The second half does not satisfy the last condition of the hourglass pattern, so the classical derivation is used there.
Notice that its bound is asymptotically worse than the bound on the first part, so the bound of the first part will dominate.
We consider M=N/2-1 to obtain the first bound, and M=N-S-2 to obtain the second bound when N ≫ S.






§ RELATED WORK









The work presented here (especially the blocked algorithms) is relevant when all N columns of the matrix A are known at once. There are applications where the columns of A are only known one column at a time, for example, for a Krylov expansion <cit.>. Our work cannot be used in this context. It can be used for Block Krylov methods though which are useful to solve 
linear system of equations with multiple right-hand
sides or eigenvalue problems, for example.



  
Variations of QR factorization algorithms
While no formal proof on a tight lower bound for the 
I/O requirements of MGS, A2V and V2Q 
existed before our work,
these kernels were experimentally identified as memory-bound bottlenecks for applications and therefore significant work
has been done to circumvent this issue. There are three main techniques, all of which compute the same result, a QR factorization. These techniques 
 are not reorderings of the operations, 
they involve mathematical transformations that fundamentally alter the graph of dependencies.


We give below
three numerical transformations that improve the I/O and concurrency bottlenecks of MGS, A2V and V2Q:
(1) by changing the algorithmic variant
of the orthogonalization scheme, (2) by “blocking” the “update” phase using the compact WY representation, (3) 
by slicing in the row dimension using the TSQR and the CAQR algorithms.

A first approach is changing the orthogonalization scheme. 
For example, instead of MGS, one can use Classical Gram-Schmidt (CGS), 
or Classical Gram-Schmit with reorthogonalization (CGS2), 
or CholeskyQR (CholQR) <cit.>, 
or Cholesky QR with reorthogonalization (CholQR2) <cit.>.
These methods have graphs of dependencies that enable lower I/O than MGS or A2V.
In these cases, the numerics, the stability, the theory, and the 
fundamental of the method 
are completely changed. But one still obtains the QR factorization.

Secondly, for MGS and A2V, it is possible to “block” the “update” phase using the Compact WY representation <cit.>. This was proposed initially for A2V, but this can easily be generalized to MGS <cit.>.
A blocking factor of size B is chosen.
The factorization proceeds with a “panel” factorization, 
then constructs the “T” matrix, 
and then the operator (projection for MGS or reflection for HH)
with T is used in a “blocked” format. 
A standard strategy is to use a block size B
so that the M-by-B panel fits in cache. The panel factorization step
of this blocked algorithm
is the algorithm (MGS or A2V) studied in this paper. 
While this slightly increases the number of floating point operations, this technique enables to 
make efficient use of (compute-bound) matrix-matrix multiplications.
A compact
WY representation of A2V is implemented 
in LAPACK in subroutine 
“https://github.com/Reference-LAPACK/lapack/blob/master/SRC/dgeqrf.fGEQRF.”
We note that https://github.com/Reference-LAPACK/lapack/blob/master/SRC/dgeqrf.fGEQRF does rely on many smaller A2V kernels for its panel factorization, 
therefore optimizing A2V is relevant for LAPACK even if the main algorithm for a QR factorization is https://github.com/Reference-LAPACK/lapack/blob/master/SRC/dgeqrf.fGEQRF.

Thirdly, CAQR <cit.> has been proposed to break the many reductions in the M dimension on the critical path of A2V.
This is particularly interesting when M is so large that no columns, 
or just a few columns can fit in cache.
The problems of A2V are two-fold: long reductions in the column dimension M, and sequential updates. 
Compact WY representation enables to remove the “sequential updates” problem.
CAQR enable to remove the long reduction in the column dimension M problem. This leads to so-called “tiled algorithm” where
the matrix A is sliced in smaller tiles, and the CAQR algorithm 
is run on the tiles. This approach enables to get I/O in
ω( MN^2/√(S). Similarly, as A2V, for CAQR, 
the Q-factor is not computed explicitly, and it is actually stored in a different format than A2V.

While these three classes of algorithms enable, in practice and in theory, better performance
than MGS, A2V or V2Q, these are not the focus of the paper. MGS, A2V and V2Q are used
in practical applications. They are relevant in practice.
They are  used in benchmarks. They are  studied in
theoretical and scientific work. They are taught in our undergraduate classes. So we judge that
any scientific contribution increasing our knowledge of these fundamental algorithms is important.

To provide an idea of what kind of performance more “efficient” algorithms can provide, 
in the experimental section,
we perform comparisons with 
MKL https://github.com/Reference-LAPACK/lapack/blob/master/SRC/dgeqrf.fGEQRF 
which performs the QR factorization using 
the compact WY representation.





  
I/O complexity and automation of the proof

The seminal work of Hong and Kung <cit.> introduced the red-blue pebble game and used it as a formalism to manually prove the I/O lower bound of programs.

Following this contribution, many papers <cit.> focused on the manual proof of lower bounds for various (classes of) programs.
A large portion of these papers consider a formalism that forbid recomputation. This assumption is necessary to decompose complex CDAGs into simpler subregions, and to be able to recombine the bounds of each region.


Irony et al. <cit.> used the Loomis-Whitney bound, a specialization of the Brascamp-Lieb theorem with canonical projections, on the  algorithm.
This was later extended by the work of Christ et al. <cit.> for the class of affine programs, by using the Brascamp-Lieb theorem.
Then, Elango et al. <cit.> added the idea of considering paths of dependencies, completing the K-partitioning method as introduced in Section <ref> of our paper.



This leads to the work of Olivry et al. <cit.> that introduced the first automatic lower bound derivation tool, based on combining bounds obtained from the K-partitioning and the wavefront methods of proof.
This tool includes several refinements to the bound derivation, by taking advantage of some specific properties.
For example, if the projections use disjoint parts of the inset, the derived bound could be improved by a constant factor.
A later extension <cit.> exploits the fact that some dimensions have very small sizes and applies this reasoning to convolutions.












§ CONCLUSION




In this paper, we introduced a new proof reasoning to derive data movement lower bounds, based on a pattern of dependencies, called the hourglass pattern.
This pattern appears on several linear algebra kernels, such as the Modified Gram-Schmidt, the QR Householder decomposition (A2V and V2Q parts), the bidiagonal matrix reduction, and the Hessenberg matrix reduction.
We managed to derive asymptotically tighter lower bounds for these kernels, compared to the previous classical methods.
We also provided tiled algorithms for MGS and QR A2V, whose amount of data movement matches asymptotically the new lower bounds, up to a constant factor.

We integrated the hourglass detection and its associated proof in the automatic data movement lower bound derivation tool of Olivry et al., IOLB <cit.>, whose implementation is freely available.
















  
Individual contributions
Fabrice identified the original problem on MGS and the intuition on how to attack this problem.
The implementation in IOLB and the proof were written by Guillaume. The proof was refined by Lionel, who also provided a detailed study of the bound, and an experimental performance analysis.
Julien identified the interesting kernels exhibiting an hourglass pattern and provided valuable insights on the linear algebra side. Lionel and Julien are also at the source of the upper bound proof in the appendix.



  
Acknowledgments Julien was partially supported for this work by NSF award #2004850.








abbrv














§ DATA MOVEMENT UPPER BOUND - TILED ALGORITHMS




In this Appendix, we present the tiled algorithm for the MGS and Householder A2V kernels.
We compute the data movement of these algorithms to obtain an upper bound on the minimal amount of data movement needed, and we show how it relates to their corresponding lower bound.





 §.§ Tiled algorithm for MGS



As shown in Section <ref>, we have an improved asymptotic lower bounds for MGS when M ≪ S.
For reference and completeness, we show in Figure <ref> a tiled left-looking ordering, based on the ideas from <cit.>. Then, we provide the proof that this ordering asymptotically matches this bound.


Let us show that, if we choose the block size B such that (M+1)B < S, then this algorithm satisfies the following properties:

    
  * it does not spill out of memory
    
  * the amount of “read” operations is 1/2 M^2  N^2 / S (leading term).
    
  * the amount of “write” is MN + 1/2N^2 (leading terms).
    
  * the total amount of I/O (read+write) is 1/2 M^2  N^2 / S (leading term).



The input of this algorithm is an M-by-N matrix A. The outputs are the M-by-N matrix Q (stored in the array A) and the upper triangular matrix R.
The first for-loop (line 1) moves along the columns of A by block. At step j0, we work with the block A(0:M,j0:j0+B).
We consider the left-looking variant of the algorithm. So, at the start of the j0 iteration, the block A(0:M,j0:j0+B) contains the initial data of the matrix A, and at the end of the j0 iteration, the same block contains the final data of the matrix Q.
The block A(0:M,j0:j0+B) stays in memory during the j0 iteration, so we read it at the start of the j0 iteration and write it at the end of this iteration.





For each vector column on the left of this block of columns (the A(0:M,i) where i<j0), we first need to project our block of A, A(0:M,j0:j0+B), onto the orthogonal complement of this column.
This is done by the first loop nest (line 3), which performs this projection one column at a time.
This requires the program to load the column A(0:M,i).
Since we are using the left-looking version of the algorithm, we only need to “read” A(0:M,i) and we do not need to “write” A(0:M,i) back in memory.

Now, let us count the number of I/O accesses of this program.
To fit the (B+1) columns of A in cache (block A(0:M,j0:j0+B) and column A(0:M,i)), we assume that B is chosen such that MB + M < S.
Over the whole course of the algorithm, moving the blocks A(1:M,j0:j0+B) into cache only happens once per block, so the I/O cost for these movements is only MN read and MN write operations.
The main I/O cost happens when reading the columns A(0:M,i) one at a time, inside the very first two for-loops (lines 1 and 3, on iterations j0 and i).
    
Notice that the I/O benefit of the blocked algorithm occurs here: every time we move the column A(1:M,i) in memory, we reuse it to project over B columns at once, instead of projecting over a single column in the unblocked version.
This leads to a reduction of the I/O by a factor B.

More precisely, the total amount of data movement involved in the reading of column A(0:M,i) is:

    [ M ×|{ j0,i | 0≤ j0 <N  and  j0  mod  B = 0  and  0≤ i< j0 }|;                = M ×|{ j0,i | 0≤ j0 <N/B  and  0≤ i< j0.B }|;        = M ×∑_j0=0^N/B j0 × B  ≈  MB(N^2/2B^2)  = 1/2MN^2/B. ]

The total data movement for the algorithm reading the columns A(1:M,i), 
(i = 0; i < j0; i++), (j0 = 0; j0 < N; j0+=B) is:

    [ ∑_(j0 = 0; j0 < N; j0+=B)∑_(i = 0; i < j0; i++) M;                  = M ∑_(j0 = 0; j0 < N; j0+=B) j0;               ≈ M B ∑_(j0 = 0; j0 < N/B; j0++) j0;                     ≈1/2 M B (N/B)^2 = 1/2MN^2/B. ]









This is indeed a factor of B lower than the number of reads of the unblocked algorithm.
    

Because we have chosen B such that M ( B  + 1)< S, we can take its maximal value to minimize the number of reads: B = ⌊S/M⌋ - 1 ≈S/M.
Thus, the I/O cost of this tiled algorithm is, assuming B ≥ 1:

    I/O ≈1/2M^2N^2/S.





  
Additional remarks
This algorithm works for any block size B, any M, any N (N < M), and any S. When we make the analysis of data movement, we further assume that M (B  + 1)< S.
So, for the I/O analysis to be correct, we need at least S > 2M. In other words, we assume that at least two columns fit in cache.

Also, we do not consider R in the data movement or in the cache constraint.
Once a value of R is computed, this value is used only once, and right after its computation.
We can thus write these values as soon as we are done with them. Hence, (1) the total number of writes due to R is N (N+1) / 2 ≈ N^2/2.  This N^2/2 “write” I/O is negligible for the algorithm (and cannot be optimized anyway), so we will  neglect it; and (2) the requirement of cache for R is 1 for the whole computation, so we do not consider R at all for the cache usage.


Note that it is also possible to follow the same reasoning to write a tiled version of the right looking variant of the algorithm.
We get a volume of I/O with a similar order of magnitude, but the constant is higher.
The right-looking variant performs more I/O than the left-looking variant, but, in addition, writes are more expensive than reads, and the I/Os of the right-looking variant are dominated by writes, whereas, the I/Os of the left-looking variant are dominated by reads. 
So, in the sequential context, we prefer to use the left-looking variant. There are some  advantages in using the right-looking variants in the sense that the updates can be done in parallel, so in a parallel multicore setting, for example, the right-looking variants would be interesting.
     




 §.§ Tiled algorithm for Householder A2V






Figure <ref> shows an implementation that matches asymptotically this bound. It is based on the same idea as the tiled algorithm for MGS in Section <ref>, but adapted to the A2V algorithm.

In particular, we show that, if we chose B such that (M+1)B < S, then our algorithm satisfies the following properties:

    
  * it does not spill out of memory.
    
  * the amount of “read” is 1/2M^2N^2 -MN^3/3/S.
    
  * the amount of “write” is MN.
    
  * the total amount of I/O (read+write) is 1/2M^2N^2 - MN^3/3/S.

We have only show the leading terms for these previous I/O costs.



The input of our algorithm is the M-by-N matrix A. The outputs are the N-by-N upper triangular matrix R and the  M-by-N unit lower triangular matrix V, that are computed in place in the M-by-N array A.

The first for-loop (line 1), using iteration k0, moves along the column of A by block.
At step k0, we work on the block A(0:M,k0:k0+B).

We consider the left-looking variant of the algorithm.
So, at the start of the k0 iteration, the block A(0:M,k0:k0+B) contains the initial data of the matrix A; and at the end of the k0 iteration, 
it contains the final data (V and R).
This block remains in cache during the k0 iteration, so we read it at the start of the k0 iteration, and we write it at the end of this iteration.

For each vector column on the left of this block (the A(0:M,j) where j<k0), we need to reflect our current block A(j:M,k0:k0+B) with the elementary Householder reflectors defined by A(j+1:M,j)
and tau[j], one column at a time.

This requires the algorithm to load A(j+1:M,j).
Since we are using the left-looking version of the algorithm, we do not need to “write” it back in memory.

Now, let us count the number of I/O accesses of this program.
To fit the (B+1) columns of A in cache, we assume that B is chosen such that M(B + 1) < S.
Over the whole course of the algorithm, moving the blocks A(1:M,k0:k0+B) in cache happens only once per block, so the I/O for this is 2MN and is negligible.
The main I/O cost happens when reading the columns A(j+1:M,j), one at a time, inside the two very first for-loops (using iteration k0 and j).

Notice that the I/O gain of this blocked algorithm against the unblocked version occurs there:
every time we load a column A(j+1:M,j), we reuse it to reflect over B columns at once, instead of a single column for the unblocked version.
This leads to a reduction of the I/O cost by a factor of B.
More precisely, the total reads of columns A(j+1:M,j) in the algorithm is:

    [ ∑_k0 = 0,  k0+=B^N-1∑_j = 0^k0-1 (M - j)  ≈ ∑_k0 = 0^N/B∑_j = 0^B.k0 (M - j);                                         ≈ ∑_k0 = 0^N/B( B.k0.M - (B.k0)^2/2);                      =  B.M. (∑_k0 = 0^N/B k0) - B^2/2 . (∑_k0 = 0^N/B k0^2);                       ≈  B.M.N^2/2.B^2 - B^2/2 . N^3/3B^3 = MN^2/2B - N^3/6B;                                                         ≈ 1/2MN^2 - N^3/3/B. ]



    ∑_(k0 = 0; k0 < N; k0+=B)∑_(j = 0; j < k0; j++) M - j 
        =     M ∑_(k0 = 0; k0 < N; k0+=B) k0 
    - 1/2∑_(k0 = 0; k0 < N; k0+=B) k0^2
       ≈    M B ∑_(j0 = 0; j0 < N/B; j0++) j0 
    - 1/2 B^2 ∑_(k0 = 0; k0 < N/B; k0++) k0^2
       ≈   1/2 M B (N/B)^2- 1/6 B^2 (N/B)^3
       ≈   1/2MN^2 -1/3 N^3/B.


This is indeed a factor of B less than the unblocked algorithm.

Because we require B to satisfy M ( B  + 1)< S, we choose B = ⌊S/M⌋ - 1 ≈S/M.
So the I/O cost of our algorithm is, assuming B ≥ 1,

    I/O ≈1/2M^2N^2 - N^3/3/S.




  
Additional remark

Just like for the MGS algorithm, the I/O analysis requires S>2M, which means that we assume that at least two columns fit in cache.




