{'fraction_non_alphanumeric': 0.041515472514593525, 'fraction_numerical': 0.008384232234217585, 'mean_word_length': 4.035179300953246, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Automated Lower Bounds on the I/O Complexity of Computation Graphs

We consider the problem of finding lower bounds on the I/O complexity of arbitrary computations. Executions of complex computations can be formalized as an evaluation order over the underlying computation graph. In this paper, we present two novel methods to find I/O lower bounds for an arbitrary computation graph. In the first, we bound the I/O using the eigenvalues of the graph Laplacian. This spectral bound is not only efficiently computable, but also can be computed in closed form for graphs with known spectra. In our second method, we leverage a novel Integer Linear Program that directly solves for the optimal evaluation order; we solve this ILP on constant sized sub-graphs of the original computation graph to find I/O lower bounds. We apply our spectral method to compute closed-form analytical bounds on two computation graphs (hypercube and Fast Fourier Transform). We further empirically validate our methods on four computation graphs, and find that our methods provide tighter bounds than current empirical methods and behave similarly to previously published I/O bounds.


Introduction
Many important applications are bottlenecked not by processing speeds, but by I/O cost: the speed to transfer data items between fast memory (e.g., registers or the CPU cache) and slow memory (e.g., RAM or disk). As a result, there has been considerable interest in designing I/O efficient algorithms and in understanding lower bounds on I/O for a given computation [11,4,15,16].
Past work on I/O lower bounds has largely focused on finding bounds for specific algorithms, such as matrix multiplication or the Fast Fourier Transform [11,4,15,16]. However, these approaches leverage properties specific to the tasks at hand, and thus do not translate across tasks. In this paper, we instead explore methods that can be applied to arbitrary computations and can be computed efficiently in an automatic fashion. Such generic bounds can thus be used to characterize the I/O cost of computations that are too complex to analyze by hand.
We approach the problem of minimizing I/O for an arbitrary computation as finding an optimal evaluation order on the underlying directed computation graph. In a computation graph, each vertex represents a single operation: the parents of the vertex indicate the operands of the operation. We assume a two-level memory architecture with a fixed amount of fast memory and infinite slow memory: I/O is incurred when transferring data between fast and slow memory (Section 3).
We present two novel methods for computing I/O lower bounds on computation graphs in this setting. Our first method bounds the I/O using the eigenvalues of the graph Laplacian (Section 4). This spectral bound is not only efficiently computable, but also can be computed in closed form for graphs with known spectra. We use it to compute closed form bounds for hypercubes and the Fast Fourier Transform (butterfly graph). We find that the resulting bound for the FFT graph is close to the previously published asymptotically tight bound.
Our second method leverages a novel Integer Linear Program that directly solves for the optimal evaluation order to minimize I/O. We solve this ILP on constant sized sub-graphs of the original computation graph to find lower bounds on the total I/O (Section 6). Like our first method, this method gives a bound that can be computed automatically for any given graph.
Finally, we evaluate our methods empirically by using them to compute lower bounds for four types of computation graphs. We compare the asymptotic growth of our results to other proposed automatic methods [11,12] and to previously published analytical bounds on computations such as Fast Fourier Transform, matrix multiplication (naïve and Strassen), and the Discrete Cosine Transform (Section 7). We find that our bounds are tighter than current automatic methods and behave similarly to published analytical bounds. The partitioned ILP method provides sharper bounds for small fast memories with well-behaved graphs, while the spectral method is more efficient to execute and can handle large fast memories and degenerate graphs.

Related Work
Hong and Kung first framed the problem of I/O complexity as the "red-blue pebble game" and used it to prove several bounds [16]. The game represents slow memory as an infinite pool of blue pebbles and fast memory as a finite set of red pebbles. An evaluation then corresponds to pebbling each vertex of the graph according to the game; I/O is incurred when placing a red pebble on top of a blue pebble (reading from slow memory) or vice-versa (writing to slow memory).
Lower bounds on naïve matrix multiplication often use the Loomis-Whitney theorem, which embeds operations in the voxels of a computation cube [15,2]. However, volume based arguments such as Loomis-Whitney do not apply for more general computations. I/O bounding techniques for algorithms beyond matrix multiplication generally focus on the computation graph itself.
Most current work on lower bounds via computation graphs requires manual inspection of the graph. In [16], the authors find a 2S partition of the computation graph to bound I/O-a proof technique that is non-trivial for complex graphs. In [22] and [5], the authors use path routing and dichotomy width respectively to find lower bounds. In [4], the authors reduce the I/O problem to a graph partitioning problem in order to find a lower bound for Strassen's matrix multiplication algorithm using the edge expansion of the graph, which was computed by hand. None of these methods can easily be computed automatically for arbitrary graphs. Instead, lower bounds on each graph must be separately proved by inspecting the specific graph. These approaches thus do not generalize and are difficult to use for highly complex graphs. We instead focus on methods that can automatically compute lower bounds for any input graph, regardless of its structure.
To our knowledge, there are only two works that discuss automated methods for lower bounds that can be applied for any graph. The first work uses an Integer Linear Program (ILP) to solve for the 2S partition of the computation graph [11]. This method is computationally expensive because it uses an exact ILP solver, and only solves for a weak heuristic (the 2S partition). The second work leverages a decomposition relaxation: the I/O cost of a computation graph can be lower bounded by the I/O cost of disjoint sub-graphs. The authors then find convex min s-t cuts on the sub-graphs [12]. We evaluate both of these automated methods in Section 7, and find that they yield looser bounds than our proposed methods.
In this paper, we present two methods for finding automated lower bounds on arbitrary computation graphs. The first is a spectral result that builds on the graph partition technique used for Strassen multiplication in [4]. Rather than manually finding the edge expansion of the computation graph, we lower bound the edge expansion using the spectrum of the graph Laplacian. This spectrum can be computed efficiently via standard numerical methods for finding eigenvalues. As a result, we can extend the edge expansion argument to graphs other than the original Strassen graph in [4]. This method also yields a new technique for closed-form bounds when the graph spectrum can be computed analytically.
Our second method builds on the decomposition theorem in [12]. However, instead of solving for a heuristic such as the min-cut or the 2S partition, we solve the I/O minimization problem on the sub-graph directly. To do this, we develop a novel ILP that can exactly solve for the actual I/O cost of a computation graph. This ILP can then be used on sub-graphs of a fixed size to provide a stronger lower bound for the total I/O cost.
Upper Bounds: Communication optimal variants of computations such as naïve matrix multiplication [15], Strassen's algorithm [3], and Fast Fourier Transform [18] have also been published. To our knowledge, our Integer Linear Program is the first automated method (polynomial time or otherwise) to find an optimal ordering given a computation graph.

Computation Graphs and Memory Model
A computation can be represented by an underlying directed computation graph G. Each operation, including the inputs and outputs, is represented by a vertex. An edge from u to v indicates that the operation v was computed with u as an operand. The graph is acylic, with the inputs as sources and the outputs as sinks. For example, the inner product of x = [x 1 , x 2 ] and y = [y 1 , y 2 ] can be represented as a 7 vertex graph: 4 vertices for each of the elements of x and y, 2 vertices for the intermediate products x 1 y 1 and x 2 y 2 , and a single vertex for the sum x 1 y 1 + x 2 y 2 ( Figure 1).
We assume a two level memory hierarchy on a single processor with infinite slow memory and a limited size fast memory of size M elements, where each operation in the computation graph is a single element. Every operation in the computation graph must be evaluated. When a vertex v is evaluated, the parents of v must be loaded into fast memory from slow memory if they are not already present. As in [4,22,11], we disallow recomputation of the same vertex: therefore, if a computed result is needed elsewhere in the computation graph and is about to be evicted, the result must first be written to slow memory.
I/O can be separated into trivial (reading inputs and writing outputs) and non-trivial I/O. We focus on non-trivial I/O: we thus do not directly include the cost of reading inputs or writing outputs. Instead, we assume that inputs can be read from the user directly into fast memory, and outputs are reported to the user immediately as they are computed. However: if an input is evicted from fast memory and is still needed elsewhere in the computation, it must be written to slow memory. This assumption is inherent in the proof in [4,12]. Because we seek lower bounds, we do not constrain the eviction policy of fast memory. I/O is incurred when, during computation, an element is written to slow memory from fast memory or read from slow memory into fast memory.
An evaluation order is then the order that operations are evaluated in the graph. Since a vertex can only be evaluated after its parents, a valid evaluation order must be topological with respect to the graph. We thus seek lower bounds on the I/O incurred by the optimal evaluation order.

Optimization Task
Formally, let G = (V, E) be a computation graph with vertices V and edges E. Let n = |V | be the number of operations in the graph, and let M be the size of fast memory. Furthermore, let V in , V out ⊆ V be the input and output vertices respectively. Note that each vertex in the graph is evaluated exactly once; therefore, the total computation takes exactly n time-steps.
We thus formalize an evaluation order on G as a permutation matrix X ∈ R n×n , where X ij = I{v j is computed at time-step i}. Let O G be the set of valid topological orders on G. Because vertices must be evaluated after their operands, X ∈ O G .
An I/O is incurred every time an element must be read into fast memory from slow memory or written to slow memory from fast memory. Let J G (X) be the number of I/Os that were incurred by evaluating G in the order specified by X on G. Given a computation graph G, we seek a lower bound on the optimal I/O incurred:

Spectral Bounds via the Graph Laplacian
In this section, we find an automated lower bound based on the eigenvalues of the graph Laplacian.
We first frame the problem as a Quadratic Assignment Problem (QAP) by counting the number of edges that cross boundaries over a graph segmentation. We then use the eigenvalues of the graph Laplacian to find a lower bound on the problem.
, and d(v) be the in-degree, out-degree, and total degree of v respectively. For vectors x, y ∈ R n , we let x, y − be the minimal scalar product of x and y (the smallest dot product of x with any permuted y). Finally, for any subset S ⊆ V , we define ∂S as the edge boundary of S:

Counting Edges over Graph Segmentation
As in [4], we count I/O over a segmentation of the graph. For any evaluation X on G, we can choose a segmentation P ⊆ 2 V that divides V into disjoint subsets of vertices such that each S ∈ P is contiguously ordered by X. P thus defines breakpoints on X. Figure 2 depicts an example of a segmentation on a graph. The numbers on the vertices indicate the evaluation order determined by X. The graph is then segmented into green, yellow, and blue segments. Each segment is contiguous with respect to the order.
Let P X be the set of valid segmentations on X. The following lemma links the segmentation to the I/O cost. Theorem 1. For fast memory size M and graph G, the optimal I/O is lower bounded by: Proof. For each subset S ∈ P , define the following adjacent sets: R S is the vertices not in S with an edge into S, and W S is the vertices in S with an edge outside of S. Ballard et. al in [4] then present the following lemma:  Proof. We summarize their proof here. To evaluate the nodes in S, the vertices in R S must be read into fast memory (or were already in fast memory before beginning computation of S). Similarly, the vertices in W S are freshly computed and needed elsewhere in the evaluation and thus must be written out or left in fast memory at the end of S. (Figure 3). Since the fast memory size is only M , at least |R S | + |W S | − 2M I/O's are incurred by evaluating the nodes in S. Summing over all S ∈ P leads to a bound on the IO incurred by G. Any P is valid so long as P splits V into components contiguous in X. Specifically, if P X is the set of valid segmentations with respect to X: It is easier to compute the number of edges crossing into and out of S rather than the vertex sets R S and W S . Ballard et. al simply bounded |R S | + |W S | by dividing by the maximum undirected degree, thus assuming a regular graph: |R S | + |W S | ≥ 1 dmax(u) |∂S|. For bounds with highly heterogenous degrees, such as matrix multiplication, this bound is relatively weak.
However, because we seek an automated method, we can provide a tighter bound because we retain access to the graph. We thus instead bound |R S | and |W S | as: .
Summing reads and writes, we have: Minimizing over all X, we get the full bound as stated in the Theorem.
Intuitively, an adversary picks some evaluation order X on G. We pick a hard segmentation P on X to maximize the I/O incurred. Ballard et. al explicitly counted the edge crossings of any fixed α sized sub-graph of the Strassen multiplication graph and then simply multiplied the bound by n/α. This relaxation is loose for graphs where the I/O is concentrated in a small portion of the vertices, and their approach cannot be automated.
Instead, we consider the edge crossings incurred by the entire segmentation, rather than a single sub-graph. In the next section, we format our problem as a Quadratic Assignment Problem using the graph Laplacian of a out-degree normalized graph. We then solve this problem in an automated fashion via the eigenvalues of the Laplacian.

Formulation via the Graph Laplacian
In Theorem 1, we solved for the minimum order over a maximum segmentation. However, since any segmentation will give us a lower bound, we can choose to split our graph into evenly sized segments. We pick some k ≤ n as our number of segments: splitting into k subsets of as equally as possible (such that the first n mod k segments have n/k + 1 vertices and the rest have n/k vertices). For an evaluation order X, let P (X,k) ∈ P X be the k-segmentation described above. Further defineŴ (k) ∈ R n×k as (Ŵ (k) ) ij = I{i ∈ P (I,k) j }. Then XŴ (k) ∈ R k×n is the partition matrix for the k-segmentation P (X,k) .
We transform our graph G into the unweighted directed graph as follows: for each directed edge (u, v) ∈ G, we add the undirected edge (u, v) toG with weight 1 dout(u) . Henceforth, we indicate the degree function, degree matrix, and adjacency matrix of the original G as d(v), D, and A respectively; we analogously denoted,D,Ã as the degree function, degree matrix, and adjacency matrix ofG.
Let L =D −Ã be the graph Laplacian ofG. We note that L is positive semi-definite, so all of its eigenvalues are nonnegative. The Laplacian is also convenient for expressing the edge boundaries of vertex subsets. Specifically, for subset S ⊆ V , let x ∈ R n be the one-hot encoding of S (i.e x i = I{v i ∈ S}). Then: Using this property we can bound the edge crossing as: .
Letting W (k) =Ŵ (k)Ŵ (k) T , and rewriting Equation 1 leads to the following theorem: . For a computation graph G and any k ≤ n with L and W (k) defined as above, J * G is lower bounded by the solution of: The optimization problem in Theorem 3 is a Quadratic Assignment Problem (QAP) [8]. We can bound the solution using the eigenvalues of L and W (k) .

Spectral Bounds
We derive the following eigenvalue bound: Proof. We relax the topological constraint X ∈ O G , and instead constrain over orthogonal X. We thus have for any k: Then by [13] for any symmetric L, W and orthogonal matrix X, where λ 1 , ..., λ n and µ 1 , ..., µ n are the eigenvalues of L and W respectively, we have tr(X T LXW ) ≥ λ, µ − . Since W (k) is a block diagonal matrix with n mod k blocks of size J n/k +1 and k − (n mod k) blocks of size J n/k , W (k) has n/k + 1 as an eigenvalue of multiplicity n mod k, n/k as an eigenvalue of multiplicity k − (n mod k), and 0 as an eigenvalue of multiplicity n − k. Therefore, we apply our lower bound as a sum of the first k eigenvalues of L: completing the proof.
The bound in Theorem 4 can be found in O(n 3 ) time. We first find the eigenvalues λ(L) in O(n 3 ). We then iterate over possible values of k which takes constant time per iteration to find the best eigenvalue.
However, we generally only need small number of eigenvalues to find a good k. Since any value of k is a lower bound, it suffices to find the h smallest eigenvalues of L. These values can be found using a method such as Lanczos-Arnoldi with time complexity O(hn 2 ): this complexity decreases even further with sparse L using sparse eigenvalue solvers.
A tighter lower bound on the QAP can also be found via a Semi-definite Program (SDP) by methods such as in [8]. SDP solvers tend to be extremely slow, so we do not cover a SDP solution in this paper.

Analytical Bounds for Specific Graphs
For computation graphs with known eigenvalues, we can compute the bound in Theorem 4 directly.
As illustrative examples, we first perform this analysis for a graph with very simple eigenvalues: a hypercube. We then apply our method to the more complex Fast Fourier Transform graph. We derive the spectrum of the FFT graph in Appendix A; to our knowledge, this is the first closed form of the spectrum of the unwrapped butterfly graph that includes multiplicity. Previously, [16] found an asymptotically tight bound of Ω( l2 l log M ) for an l level FFT through manual inspection of the 2S partition. We find that our bound of Ω( l2 l log 2 M ) is only a factor of 1/ log M off from this previously published tight bound. Examples of the hypercube and FFT graphs are shown in Figure 4.

Hypercube Graph
We first begin with a simple example of a graph with well-known eigenvalues: the hypercube. The l-dimensional hypercube has n = 2 l vertices and Laplacian eigenvalues 2i for i = 0, ..., l with multiplicity l i . If we choose k to encompass the top eigenvalues up to i = α, we have k = α i=0 l i . Letd out < l be the maximal out-degree. Then, for any α < 2 l : While any α < l would be a lower bound, for simplicity we here choose α = 1 (i.e k = l + 1). For a tighter bound we can optimize more specifically over α. Then (assuming the worst cased out = l for the second inequality):

Fast Fourier Transform Graph
The Fast Fourier Transform (FFT) computation graph is a butterfly graph. For a l point FFT, the butterfly graph B l has (l + 1)2 l vertices, which can be arranged into l + 1 columns. The butterfly graph can be defined inductively: allow B 0 to be defined as a single vertex. Then B l can be constructed as two copies of B l−1 that are joined by an extra final column of (l + 1) vertices. We derive the eigenvalues of the Laplacian of B l in Appendix A. To our knowledge, closed forms were only previously known for the wrapped butterfly graph [10]. The eigenvalues are: The smallest eigenvalue is therefore 0, but the next eigenvalues are governed by the second equation with j = 0 as long as 2i + 1 ≥ l + 1. We choose some α and set k = 2 α+1 . We compute the lowest k eigenvalues of the Laplacian of B l . Of these eigenvalues, 2 α have (with i = l − α): To compute our lower bound, we assume the other eigenvalues are 0. We note that n = (l + 1)2 l . Then we have (dividing by our maximal out-degree 2): Suppose that we set α = l − log 2 M , under the assumption that M l. Then: To see how this behaves, we can use the small angle approximation θ 2 /2 ≈ 1 − cos(θ) for small θ to get: Thus, for large M and l where M l, our bound behaves at least as well as Ω( l2 l log 2 M ). This bound is only a 1/ log 2 M factor worse than the tight lower bound for butterfly graphs: Ω( l2 l log M ), which is computed by inspection on the specific graph using S-partitions.

Integer Linear Program and Associated Lower Bounds
The spectral method finds a lower bound by counting the edges crossing graph partitions using the spectrum of the graph Laplacian. In this section, we present an alternate method based on an ILP formulation that can solve for the optimal I/O cost directly. While this exact ILP may be computationally infeasible to solve for a full computation graph, we use the partitioning method in [12] to bound the I/O by applying our ILP on small sub-graphs of a fixed size. We find that, empirically, this method can sometimes give tighter bounds than our spectral method (Section 7), although it is generally more expensive to compute.

Exact Solution via Integer Linear Program
We first present an Integer Linear Program for finding the optimal evaluation order. In addition to our evaluation order X, define new variables Y ∈ R n×n and z ∈ R n . The columns of Y indicate the contents of fast memory after each time-step: since fast memory is initially empty, we define y 0 = 0. z indicates the operations that have been written to slow memory over the course of the evaluation. For ease of notation, let x i , y i be the i'th rows of X and Y respectively. Furthermore, let A be the adjacency matrix of the computation graph G.

Objective
The I/O cost can be split into reads and writes. We define a new variable W ∈ R n×n , where the i'th row w i = y i ∧ ¬y i−1 encodes the elements in fast memory at time step i but not at time-step i − 1.
At every time step, a read is incurred for each new element in fast memory, except for the elements that were directly computed. The total number of reads is then n i=1 n j=1 W ij − X ij . The write cost is the number of elements in slow memory: 1 T (z).
Summing reads and writes, we seek to minimize over the total I/O cost:

Constraints
Beyond X ij , W ij , Y ij , z i ∈ {0, 1}, we further add the following constraints: 1. X is a permutation matrix:

Integer Linear Program
In total then we have the following optimization problem (ILP(G)): and ∀i ∈ {1, ..., n} : This problem has 3n 2 + n boolean variables. Solving this program with an ILP solver retrieves both the optimal order (X) as well as the best I/O (the objective). To our knowledge, this is the first method to explicitly find the optimal order. However, ILPs are hard to solve and can take exponential time. While this method can find the exact solution for smaller graphs using branch-and-bound methods, the problem becomes intractable for larger sizes. One immediate solution would be to relax the integer constraints. However, the resulting linear program is highly fractional and thus a relatively weak lower bound. In the next section, we present an automated lower bound that partitions the graph into smaller sub-graphs and solves the ILP for each sub-graph in isolation.

Lower Bounds via ILP on Sub-graphs
In this section, we seek a method to automatically find lower bounds on a computation graph G using the above ILP on sub-graphs of G. We take our cue from [12] and bound the I/O of G by finding the I/O of disjoint sub-graphs of G. We leverage the following lemma: The removal of an edge from G can only lower the maximum I/O cost, because the optimal evaluation order of G with the edge is a valid solution for G without the edge. Therefore, if G is simply G without an edge, we have: J * G ≥ J * G . Given disjoint subsets V 1 , ..., V k of V , we transform G by removing edges that exist between vertex subsets. The transformed graph has a lower maximum I/O cost and is made up of disconnected components G 1 , ..., G k . Therefore J We leverage this lemma by partitioning the graph G into more manageable sub-graphs and individually solving ILP(G) for each sub-graph. We refer to this method as the partitioned ILP method. Let α be the desired size of the sub-graphs such that ILP(G) takes a reasonable amount of time. The algorithm is shown in Algorithm 1.

Algorithm 1: Partitioned ILP Method for I/O lower bound.
Data: G, α Result: Return a lower bound on Since any partition maintains the validity of the lower bound, we seek partitions with few edges between partitions because these edges are lost in the relaxation. We use METIS [17] to create these partitions. We evaluate the empirical performance of this algorithm in Section 7.
Since we solve fixed-size ILP's, each sub-graph can be solved in constant time for a small enough α. Moreover, each ILP can be solved in parallel. Fixing α, the number of sub-graphs grows linearly with the graph size.

Solver
We evaluate our two lower bounds on four common computation graphs. To facilitate our evaluation, we develop a solver that traces operations during a Python computation and thus extracts a computation graph 1 . The solver inter-operates with standard arithmetic operations and supports the inclusion of custom operations (or execution on a given computation graph). We implemented the exact solution solution via our ILP, as well as lower bounds via the spectral method and the partitioned ILP. We leveraged the Gurobi MIP solver for the ILP problem/sub-problems [14].

Evaluation Computation Graphs
We evaluate the following graphs. Examples of these graphs can be found in Figure 5. This graph follows the butterfly shape discussed in Figure 4. This graph has a published bound [16] of Ω l2 l log M .

Naïve Matrix Multiplication:
We evaluate the graph formed by naïve multiplication of two n × n matrices C = AB. Specifically, we compute C ij as the dot product of the ith row of A and the jth column of B. This graph has a published bound [15] Ω n 3 √ M .

Strassen Multiplication
We evaluate the graph formed by multiplying to n × n matrices C = AB via Strassen's method. Since Strassen's method is a recursive method that splits matrices into quadrants, we evaluate on values of n that are powers of 2. This graph has a published bound [4] of M .

Discrete Cosine Transform:
We evaluate the computation graph of a type II Discrete Cosine Transform. We use the l level Fast Cosine Transform (FCT) described in [19] of a 2 l element input array. To our knowledge, there is no published I/O lower bound for this transform.
When generating our computation graphs, we treat cumulative operations (such as addition of multiple inputs) as a series of binary operations. For example, the operation a + b + c will be executed as (a + b) + c. Thus, the maximum in-degree of any node is 2.

Baselines
The only current methods for creating automatic lower bounds for any arbitrary graph are the convex min-cut method [12] and the 2S partition method [11] (see Section 2).
Convex-Min Cut: The convex min-cut method transforms the graph with respect to a vertex v into a flow problem and then finds the minimum s-t cut of the transformed graph. The method maximizes over all v in the graph. The method decomposes trivial (reading inputs and writing outputs) and non-trivial I/O, and thus fits well with our problem set-up: to include the I/O of writing inputs to slow memory, we include inputs and outputs in the graph. The runtime of this bound is O(n 5 ) where n is the number of vertices in the graph. In order to reduce runtime, the authors first partition the graph into smaller sub-graphs using METIS and run convex min-cut on each sub-graph. If C(v, G) is the minimum convex cut for G transformed with respect to v, and P is the partition reported by METIS, the authors report the bound: More details can be found in their paper. The authors suggest that each sub-graph have at most 2 * M vertices, and evaluate their bounds on a series of small, simple computation graphs with very uniform structure. However, when evaluating on more complex computation graphs such as matrix multiplication or FFT, we found that the above bound gave trivial results (J * ≥ 0) for every one of the four graphs. We hypothesize that the prescribed sub-graph size of 2 * M is too small for more complex graphs . In our evaluation, we display results of the convex min-cut method run over the entire graph (without partitioning): The above bound is linear in M for any graph. In the worst case, this bound can take O(n 5 ) time to compute. 2S-Partition: In the 2S-partition method, the authors propose an ILP for bounding the 2S partition of the entire graph. However, their formulation, unlike for the convex min-cut method above, includes reading inputs and writing outputs within the I/O count. Moreover, we find that, even under this setting without any relaxation due to decomposition, this bound is trivial for our computation graphs: simply counting the inputs and outputs provides a stronger bound than running their method. In Appendix B, we show that the bound is trivial when run over the FFT graph. We thus do not include this method in our evaluation results.

Bound Behavior with Various Graph Sizes
We first examine graph behavior with M = 4, 5 for varying graph sizes. We plot the the spectral method, the partitioned ILP method, and the convex min-cut method against the graph size parameter (the level l for FFT/DCT and the matrix side length n for matrix multiplication). For the partitioned ILP methods we choose α = 16 (the maximum size of our partitioned sub-graphs).  Larger values of α will result in tighter bounds but longer runtimes; we chose this level of α to ensure a reasonable runtime for solving each partition.
To compare against the published bounds, we also plot the computed I/O against the graph parameter term in the analytical bounds in Section 7.2. For example, we plot the computed I/O for the FFT graph against l2 l , where l is the level of the FFT graph. If our bounds follow the growth patterns of the analytical bounds, then these plots should be roughly linear. Since DCT does not have a published analytical bound, we plot against the FFT bound l2 l . Furthermore, because the spectral bounds are far more tractable for large graphs, we plot an extra few points in some plots for the spectral bound to further illustrate the growth.
For the transform graphs (FFT, DCT), the partitioned ILP provides tighter bounds than the spectral method for both the M = 4 and M = 5 case (Figures 6, 7). However, both bounds roughly grow with the analytical growth of Ω(l2 l ), since the I/O vs the published bound is linear.
For Strassen matrix multiplication, the partitioned ILP bound is tighter than the spectral bound for M = 4 but not M = 5 ( Figure 8). Specifically, the partitioned ILP bound deteriorates dramatically for M = 5 while the spectral bound does not. This deterioration is likely a result of our partition size α = 16: to properly count I/Os for a large cache size, we need to compute over a larger sub-graph. As an extreme example, if M = 16, then we could memorize the entirety of the partitioned sub-graph and would thus count zero I/O for each sub-ILP. The spectral bound increases the size of each segment as the graph grows and thus does not suffer this problem. A possible fix for this deterioration would be to increase α with M : however, an increase in α also significantly increases runtime. We note that the bounds do grow with the analytical growth of Ω(n log 2 7 ) for fixed M as the plot against the published bound is linear.
For Naïve Matrix Multiplication, the spectral bounds behave smoothly while the partitioned ILP bounds are erratic: the bound for M = 4 is not monotonically increasing and the bound for M = 5 is trivial. This is likely due to the large number of "addition trees" in the computation graph. As stated in Section 7.2, we evaluate a + b + c as (a + b) + c, resulting in a tree with 5 nodes. These trees add intermediate vertices that do not add additional I/O. When the graph is partitioned, only α = 16 nodes can exist per partition: thus the partitions with vertices contributing genuine I/O are placed with irrelevant intermediate nodes from these addition trees. Since edges are deleted between partitions, the interactions contributing to actual I/O cost are obscured. The spectral bound, which specifies the number of segments rather than the size of segments, can absorb these irrelevant within each partition, and thus behaves more similarly to the published bound.
We note that both the spectral and partitioned ILP bounds provide tighter bounds than the convex min-cut bound in all cases except for the DCT graph, where the baseline's bound is slightly tighter than the spectral bound but not the partitioned ILP.
For all four graphs, the partitioned ILP degrades significantly with increased M , while the spectral bound decreases only marginally. As explained above, this deterioration in the partitioned ILP can be attributed to the fixed size of the partition: as M grows, the solution to each sub-ILP becomes increasingly trivial. However, the spectral bound, which does not rely on fixed size partitions, does not suffer from this same deterioration. In the next section, we evaluate the behavior of the spectral bound as M increases.

Bound Behavior with Varying M
For each of the above computation graphs, we evaluate the spectral bound for a fixed size graph with increasing M (Figures 10 -13). The spectral bound does not exhibit the same deterioration for increasing M the partitioned ILP Bound. Thus, the spectral method is a more useful bound for large M . In Appendix C, we plot against the analytical bound with respect to M , again assuming the analytical bound of FFT for DCT. While the spectral bound matches the published bound for DCT and Strassen, the bound degrades more quickly than the published bounds for FFT and naïve multiplication.

Discussion
From our evaluation, we find that both the partitioned ILP method and the spectral method can effectively find lower bounds on arbitrary computation graphs. The partitioned ILP bound is relatively strong for small values of M with graphs that have consistent structure, such as the FFT or DCT graphs. However, the bound deteriorates for large M or for graphs with intermediate components that are not well connected, such as the addition trees in naïve matrix multiplication. We found that the bound matched the published analytical bounds in terms of graph size for every graph except for naïve matrix multiplication.
The spectral bound, while weaker for small M , is more versatile and is faster to compute than either the partitioned ILP or the convex min-cut; moreover, the bound does not deteriorate for large M and can handle a wider variety of graphs. For graphs with known spectra, the spectral bound can be found in closed form. We found that the bound grows with the published analytical bounds in terms of graph size for each of the graphs, and in terms of M for 2 of the 4 graphs. Finally, the spectral method and the partitioned ILP method both provided tighter bounds than the min-cut method in 3 of the 4 graphs.

Conclusion
Finding I/O bounds for general computations remains a challenging problem. In this paper, we presented two novel methods to find I/O bounds for computation graphs, which can be computed automatically to find lower bounds for arbitrary graphs or used in analytical proofs. Our first method uses the spectrum of the graph Laplacian to provide I/O lower bounds. This spectrum can be computed efficiently even for large graphs and can also be computed in closed form for some graphs, yielding a proof technique to find new closed-form bounds. Our second method develops an ILP that can solve for an I/O-optimal execution order on a graph directly and applies it on fixed-size partitions of a graph to bound total I/O in a computationally tractable manner. We evaluated these metrics empirically on four computation graphs and showed that they find tighter bounds than previous automated methods and behave similarly to published analytical bounds. We also used the spectral method to derive closed-form bounds for several graphs, including the butterfly graph for the Fast Fourier Transform.  Proof. We approach this problem by successively decomposing the butterfly graph into smaller, weighted graphs. We allow weights to exist both on edges and vertices via the weight functions ω : V × V → R and φ : V → R respectively. We define ω(u, v) = 0 if (u, v) / ∈ E. Unless otherwise stated, an edge (u, v) has weight 1. We further let φ add additional weight to a vertex unassociated with an edge: unless otherwise stated, a vertex v has weight φ(v) = 0. Then the Laplacian L of a vertex/edge weighted graph G = (V, E) is Let L(G) indicate the laplacian of G with eigenvalues λ(L(G)). Let indicate the multiset union. We define the following augmented butterfly graph B k (q). Let q ∈ {0, 1} m be a sequence of m ones or zeros for m ≤ k. Then B k (q) has k + 1 columns of 2 k−m nodes. The first k − m + 1 columns are the butterfly graph B k−m . The next m columns are connected with the following rule: for i ∈ {1, ..., m}, if q i = 1, then create edges of weight 2 (V k−m+i,j , V k−m+i+1,j ) and (V k−m+i,j , V k−m+i+1,j ) for j ∈ {1, ..., 2 k−m−1 }. Otherwise if q i = 0, then add 2 to the vertex weights of V k−m+i,j , V k−m+i+1,j , V k−m+i,j and V k−m+i+1,j . Note that B k ([]) = B k , and any augmentation of B k maintains the degrees of each of the vertices. We note that, under our formulation, B k ( [1]) is equivalent to the augmented butterfly graph mentioned in [1].
Proof. Let V 1 be the vertices V i,j and V 2 be the vertices V i,j for i ∈ {1, ..., k + 1}, j ∈ {1, ..., 2 k−m−1 } in B k (q). In Figure 15, V 1 is the top half and V 2 is the bottom half of the vertices in each graph. Note that the sub-graphs induced by V 1 and V 2 are identical. Let D , A be the degree and adjacency matrices of the sub-graph, and letÂ be the adjacency matrix between V 1 and V 2 .
Due to the symmetry of the augmented butterfly, we can decompose the graph Laplacian into quadrants as completes the proof of this lemma.
We can then inductively apply the above lemma to decompose B k into a series of flat graphs. To this end, define the graph K(q) for q ∈ {0, 1} k as follows: let K(q) be the path graph of length k + 1 with vertex weights 0 and edge weights 1. Then for each i ∈ {1, ..., k}, if q i = 0, then delete the edge (v i , v i+1 ) and add 2 to the weights of the endpoints of that edge. Examples of K(q) for B 3 can be found in Figure 16. Proof. The proof follows from induction on Corollary 7. We make the following inductive statement: for permutation q , let Q = {0, 1} k−|q | . Then we claim: λ(g(B k (q ))) = q∈Q λ(g(B k (q + q ))).
For our base case, suppose that |q | = k. Then P is the emptyset, so the inductive claim is trivially true. For our inductive step, assume that the claim is true for all q so that |q | = m (assuming m > 0). Then we prove that the claim is true for q so that |q | = m − 1. Note that according to Corollary 7, λ(g(B k (q ))) = λ(g(B k ([0] + q ))) λ(g(B k ([1] + q ))).
We note that each K(q) consists of disconnected path graphs. Since the spectra of connected components is simple the multi-set union of the spectra of each component, we can view K(q) as a multi-set of paths. In Figure 16, for example, K([1, 0, 1]) has three connected components. Define P i as the unweighted path graph with i vertices, P i as P i where one end vertex has weight φ(v i ) = 2, and path P i where both end vertices have weight φ(v 1 ) = φ(v i ) = 2. Examples of these paths can be found in 17.   where P consists of: Proof. Define K as the multiset of path components of each K(q) built from all q ∈ Q = {0, 1} k : K = q∈Q K(q). We need to prove that P = K.
(i.) The only instance of P k+1 is created by (and equivalent to) K(q) for q as all ones.
(ii.) We first examine instances of P i such that the weighted vertex is on the right. Then, since the edge immediately after the first i − 1 edges must be deleted, there are k − i edges to the left that can be either 0 or 1. Then there are 2 k−i paths of P i that appear such that the weighted vertex is on the right. Since the number of paths with the weighted vertex as the left or the right is symmetric, there are a total of 2 k−i+1 instances of P i .
iii. We examine an instance of P i that starts on a specific index. The edges on both the left and the right of the path must be deleted. Therefore, there are k − i − 1 edges that are free to be either 0 or 1. There are thus 2 k−i−1 instances of P i which begin at a specific instance. P i can begin at k − i possible valid indices, thus resulting in (k − i)2 k−i−1 instances.
Thus K = P. Lemma 7 completes the proof.
Thus, to finish our proof of Theorem 6, we simply need to find the eigenvalues of P i , P i , and P i .
We note that L(P i ) is a tri-diagonal matrix: it has 4 on the diagonal, −2 on the off-diagonals, and 2 in the lower right corner. We show that λ(L(P i )) are the odd eigenvalues of λ(L(P 2i+1 )). We can decompose L(P 2i+1 ) as: where L(P i ) R is just L(P i ) R with reversed index order such that the first vertex has extra weight. Suppose we found an eigenvector x of L(P 2i+1 ) corresponding to eigenvalue λ with the form: where y R is the reversed vector of y. Then, Therefore L(P i )y = λy so y is an eigenvector of L(P i ). Thus, if we can find i orthogonal eigenvectors of L(P 2i+1 ) of the form Equation 5, we can reduce each of those eigenvectors to find all the eigenvectors of P i .
We show that only the odd eigenvectors of P 2i+1 fall into the form of Equation 5. The eigenvector formulas for L(P l ) are (for j = 0, ..., l − 1) [21]: L(P l ) : x j (h) = 2 cos( (2h−1)jπ 2l ), h = 1, ..., l. Let l = 2i + 1, which means that l is odd. By inspection, we note that if h is odd and l is odd, then both x j (h) andx j (h) follow the form in Equation 5. That means that the eigenvalues of L(P i ) are the odd eigenvalues of L(P 2i+1 ). Plugging into the formulas for λ(L(P 2i+1 )) finds the above closed for for λ(L(P i )).
Finally, L(P i ) is Toeplitz matrix with 4 on the diagonals and −2 on the off-diagonals. The eigenvalues of tridiagonal Toeplitz matrices have a closed form [20], thus finishing the proof.

B 2S Partition on FFT Graph
In the 2S-partition method in [11], the authors propose an ILP for bounding the 2S partition of the sub-graph. Their formulation, unlike ours, combines trivial and non-trivial I/O: they include reading inputs and writing outputs in their I/O count. Under this setting, a trivial solution is to return the number of inputs and outputs. In Figure 18, we plot the bound returned by the 2S partition method and the trivial solution (the number of inputs and outputs) for the FFT graph with M = 4. The trivial solution gives a tighter bound than the 2S partition solution.

C Spectral Bound Evaluation vs. Published Bound
In Section 7.5, we evaluated the performance of the Spectral bound for increasing M . In Figure C, we plot the I/O against the published bound with respect to M . We find that the spectral bound shrinks as the published bound for DCT and Strassen, but the bound degrades more slowly than the published bound for FFT and Naïve Multiplication.