Auguste Olivry

Guillaume Iooss

Nicolas Tollenaere

Atanas Rountev

P Sadayappan

Fabrice Rastello

IOOpt: Automatic Derivation of I/O Complexity Bounds for Affine Programs *

Keywords: Compilation, I/O complexity, Polyhedral model, Convolution, Tensor Contraction

Evaluating the complexity of an algorithm is an important step when developing applications, as it impacts both its time and energy performance. Computational complexity, which is the number of dynamic operations regardless of the execution order, is easy to characterize for affine programs. Data movement (or, I/O) complexity is more complex to evaluate as it refers, when considering all possible valid schedules, to the minimum required number of I/O between a slow (e.g. main memory) and a fast (e.g. local scratchpad) storage location.

This paper presents IOOpt, a fully automated tool that automatically bounds the data movement of an affine (tilable) program. Given a tilable program described in a DSL, it automatically computes: 1. a lower bound of the I/O complexity as a symbolic expression of the cache size and program parameters; 2. an upper bound that allows one to assess the tightness of the lower bound; 3. a tiling recommendation (loop permutation and tile sizes) that matches the upper bound. For the lower bound algorithm which can be applied to any affine program, a substantial effort has been made to provide bounds that are as tight as possible for neural networks: In particular, it extends the previous work of Olivry et al. to handle multi-dimensional reductions and expose the constraints associated with small dimensions that are present in convolutions. For the upper bound algorithm that reasons on the tile band of the program (e.g. output of a * This work was

polyhedral compiler such as PluTo), the algebraic computations involved have been tuned to behave well on tensor computations such as direct tensor contractions or direct convolutions. As a bonus, the upper bound algorithm that has been extended to multi-level cache can provide the programmer with a useful tiling recommendation.

We demonstrate the effectiveness of our tool by deriving the symbolic lower and upper bounds for several tensor contraction and convolution kernels. Then we evaluate numerically the tightness of our bound using the convolution layers of Yolo9000 and representative tensor contractions from the TCCG benchmark suite. Finally, we show the pertinence of our I/O complexity model by reporting the running time of the recommended tiled code for the convolution layers of Yolo9000.

CCS Concepts: • Theory of computation → Design and analysis of algorithms; • Software and its engineering → Automated static analysis.

Introduction

Over the last few decades, the rate of improvement of the peak arithmetic (computational) performance (GFLOPs) of processors has exceeded the rate of improvement of the peak data movement bandwidth from memory. The machine balance of systems, the ratio of peak computational performance to memory performance [START_REF] Callahan | Estimating Interlock and Improving Balance for Pipelined Architectures[END_REF], has steadily increased. Unless a computation has an operational intensity (ratio of number of executed arithmetic operations to number of data for(i = 0; i < Ni; i++) for(j = 0; j < Nj; j++) for(k = 0; k < Nk; k++)

C[i][j] += A[i][k] * B[k][j];
for(i1 = 0; i1 < Ni; i1+=Ti) for(j1 = 0; j1 < Nj; j1+=Tj) for(k = 0; k < Nk; k++) for(i = i1; i < i1+Ti; i++) for(j = j1; j < j1+Tj; j++)

C[i][j] += A[i][k] * B[k][j];
Listing 1. Matrix-matrix multiplication (untiled and tiled versions)

elements moved from/to main memory) that exceeds the machine balance, the performance achieved will be limited by the memory bandwidth bottleneck.

While the computational complexity of algorithms, that is, the number of arithmetic/logic operations, is well understood, their data movement complexity, that is, the minimum required number of data movement operations to/from memory, is very challenging to characterize. The early work of Hong and Kung [START_REF] Hong | I/O complexity: The Red-Blue Pebble Game[END_REF] devised a methodology for (manually) deriving asymptotic data movement lower bounds (they used the term I/O lower bounds, and we will do so too). Recent work in program analysis [START_REF] Elango | On Characterizing the Data Access Complexity of Programs[END_REF][START_REF] Olivry | Automated Derivation of Parametric Data Movement Lower Bounds for Affine Programs[END_REF], using novel mathematical ideas [START_REF] Christ | Communication Lower Bounds and Optimal Algorithms for Programs that Reference Arrays -Part 1[END_REF], has resulted in tools for the automated derivation of I/O lower bounds for affine computations, as expressions of symbolic problem size parameters and the size of fast memory. Thus, for specific problem sizes and cache size, a lower bound on the minimum required I/O can be computed.

However, an open problem has been the assessment of the tightness of the derived lower bounds for programs. Indeed, zero is a valid lower bound on the I/O for any computation and any cache size, but it is a useless and trivial lower bound.

The only way of establishing the tightness of an I/O lower bound is to find a matching upper bound on I/O for the computation, by exhibiting an implementation that attains it. Consider the example of matrix-matrix multiplication. Hong and Kung established that for a fast memory of size 𝑆, at least Ω 𝑁 3 √ 𝑆 I/O operations would be required for any valid schedule for the standard 𝑂 (𝑁 3 ) algorithm for matrixmatrix multiplication. More recent analysis [START_REF] Michael | A Tight I/O Lower Bound for Matrix Multiplication[END_REF] established a lower bound with scaling constant as 2 𝑁 3 √ 𝑆 -3𝑆, which exactly matches the highest order term for the data movement of an optimally tiled implementation of matrix-matrix multiplication (see Listing 1).

In this paper, we present the first automated static analysis tool (IOUB) that derives parametric expressions for I/O upper bounds for affine computations. In addition, we present an improvement to the prior state of the art in automated I/O lower bounds (IOLB), of particular relevance to modeling the I/O complexity of convolutional neural networks. IOOpt, which combines the IOLB and IOUB tools, has been used to tightly bound the inherent I/O complexity of two important classes of affine computations that we highlight as case studies in this paper: the family of tensor contraction expressions, and convolutional neural networks. IOOpt can also recommend a tiled code (with appropriate loop permutation and tiles sizes) that minimizes the data movement as modeled in our formulation. The complete flow of our analysis is depicted in Figure 1.

The paper makes the following contributions:

1. Design of the first algorithm for computing a symbolic over-approximation of the data movement for a parametric (multi-dimensional) tiled version of an affine code; 2. Design of the first fully automated scheme for expressing as an operations research problem the minimization of this data movement expression; 3. Extension of the state of the art [START_REF] Olivry | Automated Derivation of Parametric Data Movement Lower Bounds for Affine Programs[END_REF] for the derivation of tight I/O complexity lower bounds in the presence of small dimensions; 4. Integration of these techniques into a tool that computes, for a class of affine computations: 1. an arithmetic complexity; 2. proved lower and upper bounds on I/O complexity; 3. a suggested tiled code that minimizes data movement.

We evaluate our tool by considering both the layers (convolutions) of Yolo9000 [START_REF] Redmon | YOLO9000: Better, Faster, Stronger[END_REF] and the main benchmarks (tensor contractions) of the TCCG suite [START_REF] Springer | Design of a High-Performance GEMM-like Tensor-Tensor Multiplication[END_REF]. To assess the tightness of our lower and upper data movement bounds, we compute those bounds for a range of cache sizes.

The rest of the paper is organized as follows: Sec. 2 presents a high-level overview of our contributions through an example. Sec. 3 provides some useful formalism along with the required background. Sec. 4 describes our upper bound algorithm. Sec. 5 describes our contributions to the lower bound algorithm. Sec. 6 reports the results of our experiments. Sec. 7 discusses related work, before Sec 8 concludes.

Overview

I/O lower and upper bounds. Let us consider a simple memory model with two levels: a small and fast memory (of size 𝑆) and a slow but infinite memory. We assume that a piece of data has to be in the fast memory in order to be used in a computation, but we have the possibility of transferring data between the two memory levels. Now, given a program, supposing we can reorder operations and have complete control over data transfers, how many transfers between these memories do we need to perform at the very least in order to compute it? This quantity is called the I/O complexity of a program. To compute it exactly is usually intractable, as one would have to simulate every possible valid schedule of operations. Hence the need to find lower and upper bounds on this complexity.

It is possible to derive a lower bound of the I/O complexity by studying the structure of the program through a graph representation called a CDAG, which represents operations and data dependencies between computations [START_REF] Ballard | Graph Expansion and Communication Costs of Fast Matrix Multiplication[END_REF][START_REF] Hong | I/O complexity: The Red-Blue Pebble Game[END_REF][START_REF] Olivry | Automated Derivation of Parametric Data Movement Lower Bounds for Affine Programs[END_REF]. However, we also wish to evaluate how tight this bound is. Thus, we seek an upper bound of the smallest volume of I/O required by a computation. This upper bound can be obtained by simply exhibiting a schedule of the computation that has this volume of I/O. Our contribution is a novel algorithm which finds such a schedule, then compares the corresponding upper bound on the I/O to the lower bound, in order to evaluate its tightness. We also improve the stateof-the-art automated lower bound method, allowing us to find matching bounds in more cases.

In the rest of this section, we apply our algorithms informally to the matrix multiplication example described in Listing 1, in order to provide some intuitions on how they work. This program has three parameters 𝑁 𝑖 , 𝑁 𝑗 and 𝑁 𝑘 , corresponding to matrix sizes. Its arithmetic complexity is 𝑁 𝑖 𝑁 𝑗 𝑁 𝑘 .

Derivation of an upper bound on the I/O complexity. Any valid schedule of a program provides an upper bound on its I/O complexity. For a given schedule, a memory allocation strategy, and fixed numerical values of memory size 𝑆 and loop bounds, it is sufficient to count loads and stores. However, the goal is to get a symbolic bound.

For the category of programs we consider, the key transformation to get an I/O-efficient schedule is tiling [START_REF] Xue | Loop Tiling for Parallelism[END_REF]. Thus, we first designed a method to compute an upper bound of the required I/O of a given tiled program: For a given tiled program with parametric tile sizes, IOUB automatically computes the corresponding I/O cost (a symbolic expression as a function of program parameters and tile sizes) along with the footprint constraint (an analytical inequality that bounds the tile footprint with the cache size 𝑆). This expression can then be used by IOUB to derive an upper bound on the actual I/O complexity, that is, to find the best tile sizes that minimize the I/O cost under the footprint constraint.

Another key component is the selection of a subset of relevant tiling schemes (loop permutation). The bottom code in Listing. 1 corresponds to the best tiling scheme found by IOUB: here, 𝑇 𝑖 and 𝑇 𝑗 are tile sizes along loop dimensions 𝑖 and 𝑗, and tiles have size 1 along dimension 𝑘. The corresponding symbolic expressions of the I/O cost and footprint constraint found by IOUB are IO = 𝑁 𝑖 𝑁 𝑗 𝑁 𝑘 1

𝑇 𝑖 + 1 𝑇 𝑗 + 1
𝑁 𝑘 , and 𝑇 𝑖 𝑇 𝑗 +𝑇 𝑖 +𝑇 𝑗 ≤ 𝑆. We would also like to derive a symbolic expression of the I/O complexity upper bound which only depends on program parameters and memory size 𝑆, but not on tile sizes. To do so, we designed a procedure to fix tile sizes as functions of 𝑆, using computer algebra (see Sec. 6 for more details). For matrix multiplication, the final expression is

𝑈 𝐵 = 𝑁 𝑖 𝑁 𝑗 2𝑁 𝑘 √ 𝑆 + 1 -1 + 1 .
As this underlying optimization problem is difficult in general, TileOpt (see Fig. 1) can look for a numerical solution by using fixed given values for 𝑁 {𝑖,𝑗,𝑘 } and 𝑆. As an example, for 𝑁 𝑖 = 2000, 𝑁 𝑗 = 𝑁 𝑘 = 1500, and 𝑆 = 1024, IO is minimized for 𝑇 𝑖 = 𝑇 𝑗 = 31. For those values, UB = 296322580.

Lower bound and tightness. Now that we have automatically derived a symbolic upper bound on the I/O cost for matrix multiplication, we can compare it with the lower bound obtained from the IOLB tool [START_REF] Olivry | Automated Derivation of Parametric Data Movement Lower Bounds for Affine Programs[END_REF] to evaluate its tightness. The lower bound we find is:

𝐿𝐵 = 2𝑁 𝑖 + 𝑁 𝑘 + 2𝑁 𝑗 + 2𝑁 𝑖 𝑁 𝑗 𝑁 𝑘
√ 𝑆 where 𝑆 is the size of the fast memory. We can check that this bound is asymptotically tight (see Sec. 6).

However, these bounds are not always tight. For example, on a convolution (as shown in Sec. 5.4), the lower bound derived by the currently published version of IOLB is several asymptotic orders below the derived upper bound. This lower bound can be improved by taking into account (i) the reductions across multiple dimensions, which refines the dependence analysis performed in the first steps of the algorithm, and (ii) small dimensions. Small dimensions are dimensions with a much smaller number of iterations than the size of the fast memory 𝑆. This property can be used to improve the power factor of the asymptotic bound. Both improvements to the IOLB algorithm are described in Section 5, and allow us to obtain asymptotically tight bounds for tensor contraction and convolution computations.

Lower vs. upper bound method. Both methods are based on polyhedral compilation tools, for abstract program representation and to compute cardinalities of polyhedra with Barvinok's algorithm [START_REF] Verdoolaege | Counting Integer Points in Parametric Polytopes Using Barvinok's Rational Functions[END_REF]

Background

In this section, we present the class of programs, its representations, and the memory model we consider in our analysis.

Class of Programs

Imperfectly nested affine loop programs. We consider imperfectly nested loop programs, which can be recursively defined as a sequence of statements and for loops, each being itself composed of imperfectly nested loops. A program parameter is a symbolic constant during the compilation, whose value will be known during the execution of the program (for example, array sizes). The access functions made by statements to arrays are affine expressions of surrounding loop indices and program parameters. We also assume that conditions on loop indices are affine expressions of surrounding loop indices and program parameters [START_REF] Bastoul | Putting Polyhedral Loop Transformations to Work[END_REF].

For example, the 1D-convolution described in Listing 2 matches these criteria. This program is even perfectly nested: all loops are nested with a single statement inside. Constants Nc, Nf, Nx and Nw are the program parameters. All loop indices are bounded by affine constraints (e.g. 0 ≤ 𝑐 < Nc), and all array subscripts are affine expressions (e.g. (𝑥 + 𝑤, 𝑐) for array Image).

Fully-tilable program. The tiling transformation [START_REF] Xue | Loop Tiling for Parallelism[END_REF] is a key loop transformation to improve the data locality of a program. Given a list of consecutive dimensions (called loop band), this transformation groups their iterations into tiles, which are executed atomically. Listing 1 shows an example of tiling transformation for the matrix multiplication example: the tiled dimensions are i and j and the iterations are grouped into rectangles of size Ti × Tj (tile shape). Note that when the tiling shape is rectangular, each tiled dimension (e.g. i) is strip-mined, giving rise to the tile dimension that will iterate over the tiles (e.g. i1), and the local dimension that iterates inside a tile (e.g. i in the transformed program).

A tiling is legal when there is no cycle of dependencies between the computation of different tiles, that is, when the atomicity condition between tiles can be respected by the schedule of the program. The algorithms presented in this paper are described on the assumption that the input program is fully tilable using rectangular tiles, which means that all its dimensions can be legally tiled by using rectangular tiles.

In theory, any affine program could be pre-processed using a polyhedral compiler to provide a fully permutable (that is, tilable using rectangular tiles) loop band, but the reality is more complex, as analyzing a non-regular iteration domain involves being able to cope with the simplification and solving of a complex system of symbolic expressions. Nevertheless, several important classes of computation fit our simplified hypothesis, such as a convolution, all kinds of tensor contractions and many linear algebraic kernels, including matrix multiplication. Listing 2 shows another example of such computation, which is a simplified convolution. It will be used as our running example.

Program Representation

In this paper, we will mainly use two representations of a program. The CDAG representation is an unrolled graph representation of the program, and is used as the base formalism for proving lower bounds. The polyhedral representation is a mathematical and concise representation of the structure of the program, and corresponds to the intermediate representation manipulated by our algorithms.

Computational Directed Acyclic Graph (CDAG). This representation will be used by our data movement model. Definition 3.1. (CDAG) A computational directed acyclic graph G = (𝑉 , 𝐸, 𝑂) is a graph formed by:

• a finite set of nodes 𝑉 , each one representing an elementary computation, producing a piece of data, • a finite set of edges 𝐸, representing data dependencies:

there is an edge from a node 𝑣 to a node 𝑤 if the computation at 𝑤 requires the data from 𝑣. The input data of the program is represented by the nodes with no incoming edge. The output data is represented by a set 𝑂 of nodes, potentially with no outgoing edges.

We will consider a subset of a CDAG, which is itself a CDAG: Given a sub-graph

(𝑉 ′ , 𝐸 ′ ), the corresponding sub- CDAG is (𝑉 ′ , 𝐸 ′ , 𝑂 ′ ) where 𝑂 ′ = 𝑂 ∩ 𝑉 ′ .
Polyhedral representation. The iteration domain I 𝑆 of a statement 𝑆 is the set of integral values that the surrounding loop indices take during execution. Each point of this space is associated with an execution instance of the statement. Due to the hypotheses made on the loop bounds, the iteration domain of a statement is a Z-polyhedron, that is, a set of integral points whose coordinates satisfy a set of affine constraints. The computation of the cardinality of a set 𝐸 (denoted |𝐸| -symbolic expression as a function of the program parameters), can be efficiently done using the Barvinok algorithm 1 [START_REF] Barvinok | A Polynomial Time Algorithm for Counting Integral Points in Polyhedra when the Dimension is Fixed[END_REF] available in any polyhedral framework library such as ISL [START_REF] Verdoolaege | ISL: An Integer Set Library for the Polyhedral Model[END_REF].

Each array 𝐴 is associated with a memory domain M 𝐴 , which is the multidimensional set of valid array indices for the array 𝐴. Given a statement 𝑆 containing a read or a write to an array 𝐴, this occurrence is associated with a memory access function 𝑓 𝐴 : D 𝑆 ↦ → M 𝐴 , which is an affine function. For any given sub-domain 𝐷 𝑆 ⊂ I 𝑆 defined as a Z-polyhedron, the associated data footprint 𝑓 𝐴 (𝐷 𝑆 ) is also a Z-polyhedron that can be easily computed by polyhedral compilation tools.

For the matrix-matrix multiplication example (Listing 1, untiled version), the dimensions are D = {𝑖, 𝑗, 𝑘 } and the iteration domain is

I = {(𝑖, 𝑗, 𝑘) | 0 ≤ 𝑖 < 𝑁 𝑖 ∧ 0 ≤ 𝑗 < 𝑁 𝑗 ∧ 0 ≤ 𝑘 < 𝑁 𝑘 }. A 3-dimensional rectangular subset of this domain will be defined as [𝑙 𝑖 , 𝑢 𝑖 ] × [𝑙 𝑗 , 𝑢 𝑗 ] × [𝑙 𝑘 , 𝑢 𝑘 ],
where, for each dimension 𝑑 ∈ D, 0 ≤ 𝑙 𝑑 ≤ 𝑢 𝑑 < 𝑁 𝑑 . The example code also contains three array accesses:

• for array 𝐴, with domain

M 𝐴 = {𝑥, 𝑦 | 0 ≤ 𝑥 < 𝑁 𝑖 ∧ 0 ≤ 𝑦 < 𝑁 𝑘 }

Memory Model

In order to model the amount of data movement needed, we use the red-white pebble game, which is a variation of Hong & Kung's red-blue pebble game [START_REF] Hong | I/O complexity: The Red-Blue Pebble Game[END_REF]. It was introduced in Olivry et al. 's work [START_REF] Olivry | Automated Derivation of Parametric Data Movement Lower Bounds for Affine Programs[END_REF] to provide a convenient framework for deriving I/O lower bounds, and in particular to model the no recomputation assumption. Since this work also deals with lower bounds, we use this formalism for both upper and lower bounds. We consider a hierarchy of two memories:

• An infinite slow memory which contains the input data at the beginning of the execution. • A fast scratchpad of limited capacity 𝑆, which must contain the data used by a computation.

In the red-white pebble game, pebbles can be placed on the nodes of the CDAG. A white pebble represents a data which has been computed, and a red pebble represents a data which is present in the fast memory. To account for the limited size of the fast memory, there are only 𝑆 red pebbles. There is an unlimited number of white pebbles.

We can place or remove a pebble of a given color according to the following set of rules, assuming that the limit constraint on the number of simultaneously used red pebbles is fulfilled: 1 In theory, Barvinok's algorithm has an exponential complexity in the number of dimensions of the set. In practice, its execution time is usually less than a few seconds, when applied to the polyhedral sets commonly encountered during program analysis.

for(c1 = 0; c1 < Nc; c+=Tc)

for(f1 = 0; f1 < Nf; f+=Tf) for(x = 0; x < Nx; x++) for(c = c1; c < c1+Tc; c++) for(w = 0; w < Nw; w++)

for(f = f1; f < f1+Tf; f++) Out[f][x] += Image[x+w][c] * Filter[f][w][c];
Tile limit Listing 3. Tiled code for 1D-convolution for tiling schedule (𝑤, 𝑐, 𝑓 , 𝑥) , 𝑇 𝑐 = Tc,𝑇 𝑓 = Tf,𝑇 𝑥 = 1,𝑇 𝑤 = Nw

• Fetch rule: A red pebble can be placed on any node that has a white pebble. • Spill rule: A red pebble can be removed from any node. • Computation rule: If a node has no white pebble and all its immediate predecessors have red pebbles, then a red and a white pebble can be placed on this node. A valid game sequence is a list of applications of these rules (moves) on specific nodes and edges of the CDAG, such that:

• At the start of the game, there is a white pebble placed on every input node. • At the end of the game, all the nodes of the graph are covered by white pebbles, • At any point of the game, there are at most 𝑆 nodes with red pebbles. A game sequence mirrors the data movement that can be performed by a computation. Thus, the amount of loads from the slow memory done by a computation corresponds to the number of times the fetch rule is used in a valid game sequence.

Upper Bound on Data Movement

In this section, we first present an algorithm that computes an over-approximation of the required I/O for a given parametric tiled program. We then show how to use this symbolic expression to select a permutation and tiles sizes that minimize the amount of I/O. As an actual valid schedule is associated with the computed I/O, this provides a valid upper bound for the data movement complexity. An interesting side effect of the method is that the computed optimizing schedule can be exposed to the user as a suggested loop transformation.

Loop Permutation and Tiling

Let us consider the example of Listing 1 that reports a tiled and non-tiled code for matrix-matrix multiplication. The original code contains three dimensions 𝑖, 𝑗, and 𝑘. The tiled code contains five loops on dimensions 𝑖, 𝑗, 𝑘, 𝑖, 𝑗, from outer to inner. The three outermost loops span the entire iteration domain of size Ni×Nj×Nk, while the two innermost ones span a polyhedron of size Ti × Tj. We refer to the innermost (resp. outermost) part as the intra-tile (resp. inter-tile) dimensions. As we will see later, unlike the permutation of the inter-tile dimensions, that of the intra-tile dimensions is not relevant in our model. In other words, while a different schedule with loops ordered as (from outer to inner) (i1,k,j1,i,j) could lead to a different I/O estimation than for the schedule with loops ordered as (i1,j1,k,i,j) (as in Listing. 1), our cost model will not be affected if we permute 𝑖 and 𝑗 as in loop order (i1,j1,k,j,i). This motivates the following notation to represent the tiling schedule of Listing. 1: ((𝑖, 𝑗, 𝑘), 𝑇 𝑖 = Ti,𝑇 𝑗 = Tj,𝑇 𝑘 = 1 where the ordered list P = (𝑖, 𝑗, 𝑘) describes the permutation of the intertile loop dimensions, and T = {𝑇 𝑖 = Ti,𝑇 𝑗 = Tj,𝑇 𝑘 = 1} describes the tile sizes. Similarly, for the example in Listing 2, the tiling shown in Listing 3 is represented as ((𝑤, 𝑐, 𝑓 , 𝑥), 𝑇 𝑐 = Tc,𝑇 𝑓 = Tf , 𝑇 𝑥 = 1,𝑇 𝑤 = Nw}). Note that the loop on 𝑤 in the inter-tile loops (outermost loops), and the loop on 𝑥 in the intra-tile loops are both omitted in the code as they would have only one iteration. The permutation of intra-tile loops in the code is arbitrary.

More formally, a tiling schedule is defined as a tuple (P, T ) where P = 𝑑 𝑗 | D | ≥ 𝑗 ≥1 is a permutation of dimensions D representing the inter-tile loop order, and T = {𝑇 𝑑 } 𝑑 ∈D represent the tile dimensions. In P, 𝑑 | D | represents the outermost loop dimension, while 𝑑 1 represents the innermost loop dimension that encloses the tile: the permutation order is outer (leftmost) to inner (rightmost). The tile size for dimension 𝑑 is 𝑇 𝑑 .

For an inter-tile loop dimension 𝑑 𝑗 , we define its subdomain as the set of points in the iteration domain such that the indices of the enclosing loops 𝑑 𝑘 (|D| ≥ 𝑘 ≥ 𝑗) have fixed values. Taking the example of the 1D-convolution of Listing 3, the sub-domain SD 𝑑 2 (c1,f1) for loop dimension

𝑑 𝑗 = 𝑑 2 = 𝑓 (recall that 𝑑 4 = 𝑤, 𝑑 3 = 𝑐, 𝑑 2 = 𝑓 , 𝑑 1 = 𝑥) is {(𝑤, 𝑐, 𝑓 , 𝑥) | 0 ≤ 𝑥 < Nx ∧ c1 ≤ 𝑐 < c1+Tc ∧ 0 ≤ 𝑤 < Nw ∧ f1 ≤ 𝑓 < f1+Tf}.
More generally, if we denote by ik the fixed index value at level 𝑘:

SD 𝑑 𝑗 (iD, . . . , ij) = {𝑖 | D | , . . . , 𝑖 1 ∈ I | ∀𝑘 ≥ 𝑗, ik ≤ 𝑖 𝑘 < ik + 𝑇 𝑑 𝑘 }
A Sub-Domain data Footprint for an array 𝐴 at level 𝑗 is defined as:

SDF 𝐴,𝑗 (iD, . . . , ij) = 𝑓 𝐴 SD 𝑑 𝑗 (iD, . . . , ij)
As an example, for the 1D-convolution running example we have:

SDF Image,2 (c1,f1) = (Nx+Nw-1) × Tc
Whenever there is a non-empty overlap of data used between two consecutive sub-domains, estimating it allows us to refine our cost model. We define the inter-Sub-Domain Reuse for an array 𝐴 at level 𝑗 as:

SDR 𝐴,𝑗 (iD, . . . , ij) = 𝑓 𝐴 SD 𝑑 𝑗 (iD, . . . , ij) ∩ 𝑓 𝐴 SD 𝑑 𝑗 (iD, . . . , ij -𝑇 𝑑 𝑗 )
A sharp eye would have observed that the notion of reuse from a "previous" sub-domain is meaningful for all subdomains but the first one in the loop. To handle this subtlety, for a given loop dimension 𝑑 𝑗 (assuming its loop index starts at 0), we split the iteration space into two sub-domains: the front domain

𝐼 front = {(𝑖 | D | , . . . , 𝑖 1 ) ∈ I | 𝑖 𝑗 = 0}
and the back domain

𝐼 back = {(𝑖 | D | , . . . , 𝑖 1 ) ∈ I | 𝑖 𝑗 ≠ 0}
This allows us to define, for a given sub-domain and an array 𝐴, its inverse density as the ratio:

ID 𝐴,𝑗 (iD, . . . , ij ≠ 0) = SDF 𝐴,𝑗 (iD, . . . , ij) -SDR 𝐴,𝑗 (iD, . . . , ij) SD 𝑑 𝑗 (iD, . . . , ij) ID 𝐴,𝑗 (iD,. . . ,0) = SDF 𝐴,𝑗 (iD, . . . , ij) SD 𝑑 𝑗 (iD, . . . , ij)
As we will see later, the inverse density is an over-approximation of the best attainable inverse operational intensity (data movement per computation unit). As an example, for the 1D-convolution example,

SD 𝑥 (c1,f1,x) = {𝑐, 𝑓 , 𝑤, 𝑥 | 0 ≤ 𝑤 < Nw ∧ 𝑥 = x ∧ c1 ≤ 𝑐 < c1+Tc ∧ f1 ≤ 𝑓 < f1+Tf} |SD 𝑥 (c1,f1,x)| = Nw × Tc × Tf SDF Image,1 (c1,f1,x) = Nw × Tc SDR Image,1 (c1,f1,x) = Tc × (Nw -1) ID Image,1 (c1,f1,x ≠ 0) = 1/(Nw × Tf) ID Image,1 (c1,f1,0) = 1/Tf
To avoid complicated expressions for non-rectangular domains, we consider the maximum value for the inverse density that we specialize for the front and the back:

ID front 𝐴,𝑗 = max iD,...,ij/ij=0
ID 𝐴,𝑗 (iD, . . . , ij)

ID back 𝐴,𝑗 = max iD,...,ij/ij≠0
ID 𝐴,𝑗 (iD, . . . , ij)

In many cases, the sub-domain sizes and intersections do not depend on loop indices, and there is no need to take the maximum.

Cost Model

Let us start with the case where there is a single array 𝐴. We consider a tiling (P = (𝑑 | D | , . . . , 𝑑 1 ), T = {𝑇 𝑑 } 𝑑 ∈ D ).

First, the data footprint of a tile has to be smaller than the cache capacity 𝑆:

SDF 𝐴,1 ≤ 𝑆. Second, let us define the outermost reuse dimension as the "first" (smallest possible 𝑙) dimension 𝑑 𝑙 in P such that the sub-domain data footprint of 𝐴 is less than the cache capacity 𝑆. That is, 𝑙 is the leftmost index in P such that SDF 𝐴,𝑙 ≤ 𝑆.

The total I/O cost for array 𝐴 derived from the inverse density as follows constitutes an upper bound on the cost of an optimal red-white pebble game:

𝐼𝑂 𝐴 = ID front 𝐴,𝑙 × |𝐼 front | + ID back 𝐴,𝑙 × |𝐼 back | Indeed,
by construction, each sub-domain can be executed by bringing only once (before starting execution) all the required data. Every sub-domain in the back can reuse the data used in the previous sub-domain and the corresponding I/O can thus be saved. This means that the optimal inverse operational intensity for any sub-domain in the front (resp. in the back) is no more than ID front 𝐴,𝑙 (resp. ID back 𝐴,𝑙 ). When there are multiple arrays 𝐴 1 , . . . , 𝐴 𝑠 , we "cut" the cache into array-specific regions of sizes 𝑆 1 , . . . , 𝑆 𝑠 , such that

𝑆 1 + • • • + 𝑆 𝑠 = 𝑆.
The conditions are the same as above, replacing 𝐴 and 𝑆 by 𝐴 𝑗 and 𝑆 𝑗 .

On our running example, for array Image, supposing the outermost reuse dimension is 𝑥: The red-white pebble game and the corresponding lower bound reasoning frameworks can be extended to multi-level memory hierarchies [START_REF] Savage | Extending the Hong-Kung Model to Memory Hierarchies[END_REF]. The above computation of I/O can also be extended by simply considering one tiling band per cache level and independently applying the previous reasoning to each level.

𝐼𝑂 Image = ID front Image,1 × |𝐼 front | + ID back Image,1 × |𝐼 back | = 1

Loop Permutation Selection

We have only focused on computing the required I/O cost for a given permutation of the inter-tile loops. Having a symbolic expression as a function of the tile sizes allows the use of a Non-Linear optimization Problem (NLP) solver such as IPOPT [START_REF] Wächter | On the Implementation of an Interior-Point Filter Line-Search Algorithm for largescale nonlinear programming[END_REF] to find optimizing tile sizes for this particular permutation. Unfortunately, the search space of all possible permutations grows exponentially with the number of dimensions and with the memory depth.

However, it is straightforward to see that many permutations are equivalent in terms of I/O cost: for example, switching the two outer loops c1 and f1 in the conv-1D example in Listing 3 would not change the I/O costs. Moreover, some permutations are better than others: for example, it is not worth exploring the permutation whose innermost dimensions do not allow any data reuse between two successive tiles. In this section, we provide some insights into how to select a subset of relevant permutations for a given code.

First, we define a notion of reuse: for a given dimension 𝑑 and an array 𝐴, there is reuse for array 𝐴 on dimension 𝑑 if, when putting dimension 𝑑 innermost (that is, setting 𝑑 1 = 𝑑), the sub-domain data footprint at level 2 for 𝐴 does not increase much compared to level 1:

SDF 𝐴,2 -SDF 𝐴,1 ≪ SDF 𝐴,1
This allows us to decide which dimensions are worth putting at the innermost levels.

In simple cases, this notion is not ambiguous: for array Out in the conv-1D example, 𝑐 is a reuse dimension as setting it innermost would lead to SDF Out 2 = SDF Out 1 = Tf × Tx, while 𝑓 is not a reuse dimension as setting it innermost would lead to SDF Out 2 -SDF Out 1 = (Nf -Tf) × Tx (unless Nf = Tf but in that case it means that 𝑓 iterates only once). For more complex subscript expressions, such as for array Image, the criterion is less clear: setting 𝑤 as innermost and assuming it is not part of the tile would lead to comparing Tw -1 with Tx (as SDF Image 2 = (Tx + Tw -1) × Tc, SDF Image 1 = Tx × Tc). In other words, a reuse criterion requires an "oracle" (such as the user) to provide information about small and large dimensions.

Assuming such an oracle and a reuse criterion, Algorithm 1 selects a set of representative permutations. This algorithm builds a permutation from the innermost to the outermost dimensions while keeping track of the set of remaining dimensions (variable D ′ ), and of the set 𝑆 of arrays having a potential reuse along each dimension 𝑑 (variable 𝑅). When a dimension is chosen, the set of potential reuse is updated.

When there is no more potential reuse, an arbitrary permutation of the remaining dimensions is chosen.

Figure 2 illustrates the steps of the algorithm on the 1D convolution (Listing 2), and how it narrows the loop permutation space to three permutations (one of them is the one used for the running example in Listing 3). Note that some permutations are pruned during the selection process (red cross in Figure 2), because they have strictly less reuse potential than others. For example, selecting dimension 𝑐 as the innermost dimension allows us to have reuse on array Out, but dimension 𝑤 allows us to have reuse on one more array (Out and Image).

Putting It All Together

The IOUB tool combines the permutation selection and I/O cost computation as follows: it first selects a list of permutations and generates the corresponding set of tiled versions (with tile sizes as parameters) using Algorithm 1. For each permutation, it computes a symbolic expression (as a function of program parameters, tile sizes, and cache size) that represents the required I/O, using the method described in Sec. 4.1 and 4.2, as well as constraints on tile sizes. This is implemented using the ISL [START_REF] Verdoolaege | ISL: An Integer Set Library for the Polyhedral Model[END_REF] and Sympy [START_REF] Meurer | SymPy: Symbolic Computing in Python[END_REF] libraries.

Then, this can be fed to TileOpt, which uses an NLP solver (IPOPT [START_REF] Wächter | On the Implementation of an Interior-Point Filter Line-Search Algorithm for largescale nonlinear programming[END_REF]) to find the best tile sizes for each permutation. The final I/O cost is the minimum over all versions, and TileOpt also provides a basic tiled code that implements the corresponding tiling scheme. It can also be useful to have a symbolic bound that does not depend on tile sizes, but only on program parameters and cache size. A discussion on how to derive such an expression for the examples we used in experiments can be found in Sec. 6.

Lower Bound on Data Movement

In this section, we consider the problem of finding a symbolic lower bound on the volume of loads needed to perform an affine computation. We first present the main intuitions behind the partitioning method, which is one of the state-of-theart techniques to derive a symbolic lower bound [START_REF] Christ | Communication Lower Bounds and Optimal Algorithms for Programs that Reference Arrays -Part 1[END_REF][START_REF] Hong | I/O complexity: The Red-Blue Pebble Game[END_REF][START_REF] Olivry | Automated Derivation of Parametric Data Movement Lower Bounds for Affine Programs[END_REF].

We then provide two improvements on this method, namely reductions and small dimensions. These are needed to find asymptotically tight bounds for some programs, in particular convolutions. Indeed, let us consider the 2D convolution from Figure 3a. For the sake of examining the asymptotic behavior of the bound, let us assume that all program parameters are equal to 𝑁 and greater than the small memory size 𝑆. With the state-of-the-art method, the lower bound that is found is 𝑂 (𝑁 4 ). If we detect the reduction and adapt the dependence analysis before applying the partitioning method (Section 5.3), the bound is improved to 𝑂 (𝑁 7 /𝑆). We can further improve this bound by exploiting the fact that 𝐻 and 𝑊 are usually small parameters that are much smaller than 𝑆. By refining the partitioning method (Section 5.2), we obtain 𝑂 ( √ 𝐻𝑊 𝑁 5 / √ 𝑆) as a lower bound, and this bound is asymptotically tight.

Background: The Partitioning Method for IOLB

The partitioning method considers a red-white pebble game over the CDAG of a computation (both notions were introduced in Section 3). A complete mathematical formalization of the content of this subsection can be found in [START_REF] Olivry | Automated Derivation of Parametric Data Movement Lower Bounds for Affine Programs[END_REF].

Another key notion is a 𝐾-bounded set, which is a set of nodes of the CDAG with at most 𝐾 predecessors. If we consider a (𝑆 + 𝑇 )-bounded set (where 𝑆 is still the size of the fast memory and 𝑇 > 0), then we need to load at least 𝑇 pieces of data in the fast memory to be able to perform the computation. This gives a lower bound on the amount of communication needed to compute this set. Any valid sequence of moves in the red-white pebble game can be decomposed into a (𝑆 + 𝑇 )-partition, i.e., a partition of the CDAG into (𝑆 + 𝑇 )-bounded sets. Therefore, a lower bound on the total number of loads can be deduced from an upper bound on the number of sets in a (𝑆 + 𝑇 )-partition.

When the dependencies of the computation are affine, we can find an upper bound on the size of a (𝑆 +𝑇 )-bounded set by using a geometrical argument called the Brascamp-Lieb inequality. Its mathematical formulation is: Theorem 5.1 (Brascamp-Lieb inequality). Let 𝑑 and 𝑑 𝑗 be non-negative integers and 𝜙 𝑗 : Z 𝑑 ↦ → Z 𝑑 𝑗 be group homomorphisms for all 1 ≤ 𝑗 ≤ 𝑚.

If we have a set of 𝑠 𝑗 ∈ [0, 1] such that, for all subgroups H of Z 𝑑 :

𝑟𝑎𝑛𝑘 (H ) ≤ 𝑚 ∑ 𝑗=1 𝑠 𝑗 • 𝑟𝑎𝑛𝑘 (𝜙 𝑗 (H ))
Then, for all nonempty finite sets 𝐸 ∈ Z 𝑑 :

|𝐸| ≤ 𝑚 𝑗=1 |𝜙 𝑗 (𝐸)| 𝑠 𝑗
Here, 𝑟𝑎𝑛𝑘 (𝐺) is the subgroup rank, defined as the smallest cardinality of a generating set for 𝐺. This can be seen as an equivalent to the dimension of a vector space in a discrete setting.

This result can intuitively be seen as a generalization of the following inequality for 3-dimensional vector spaces:

If the surfaces of all three projections of a 3-dimensional volume 𝑉 on planes 𝑥 = 0, 𝑦 = 0, 𝑧 = 0 are bounded by some constant 𝐶, then |𝑉 | ≤ 𝐶 3/2 .

Vertices in the CDAG are mapped to points in a multidimensional geometric space E ≃ Z 𝑑 through some mapping 𝜌 (where dimensions are typically loop indices), and regular data dependencies in the CDAG are represented as projections on a lower-dimensional space. The condition "set of vertices 𝑃 ⊂ 𝑉 is (𝑆 +𝑇 )-bounded" in the CDAG corresponds to a condition of the form "the size of the projections of 𝜌 (𝑃) in E is bounded by 𝑆 + 𝑇 ". Finding a bound on the size of a (𝑆 + 𝑇 )-bounded set in a CDAG can thus be reduced to finding a bound on the size of a set 𝐸 in a geometric space, given the cardinality bounds on some of its projections.

This theorem is applied to a 𝐾-bounded set 𝐸 (where 𝐾 = 𝑆 + 𝑇 ). We select the group homomorphisms 𝜙 𝑗 such that their images are included in the input set of 𝐸. Indeed, by definition of a 𝐾-bounded set, we will have |𝜙 𝑗 (𝐸)| ≤ 𝐾, for all 1 ≤ 𝑗 ≤ 𝑚. Therefore, we will obtain the upper bound: |𝐸| ≤ 𝐾 𝜎 where 𝜎 = 𝑗 𝑠 𝑗 .

To select such homomorphisms 𝜙 𝑗 , we study the paths of affine dependence of the CDAG involving 𝐸 that originate with a vertex in the input set of 𝐸. These paths are associated by a path relation, which is the composition of the affine functions of the dependencies composing a path. We build a homomorphism from each path relation.

Then, we need to find a set of rationals 𝑠 𝑗 that satisfy the linear constraints from Theorem 5.1, and minimize 𝜎. In order to find the subgroups H that lead to the tightest constraints on 𝑠 𝑗 , we use bases of 𝐾𝑒𝑟 (𝜙 𝑗 ). We select a combination of vectors of these bases, and consider the subgroup generated by these vectors to build a constraint on 𝑠 𝑗 .

Once we have a symbolic upper bound of the size of a (𝑆 + 𝑇 )-bounded set, we obtain a symbolic lower bound on the volume of data communications by the program, for any schedule.

Adapting the Partitioning Method to Small Dimensions

We assume that some of the dimensions of the iteration domain are small, i.e., their size is orders of magnitude smaller than 𝐶. We note the product of all the sizes of the small dimensions 𝑁 𝑠𝑑 . We adapt the partitioning method by considering an additional group homomorphism 𝜙 𝑠𝑑 to the application of Theorem 5.1. This group homomorphism is a projection of the space on the small dimensions, which gives us a much tighter bound |𝜙 𝑠𝑑 (𝐸)| ≤ 𝑁 𝑠𝑑 . Assuming a set of 𝑠 𝑗 and 𝑠 𝑠𝑑 (coefficient associated with 𝜙 𝑠𝑑 ) satisfying the constraints of Theorem 5.1, we obtain the following upper bound:

|𝐸| ≤ 𝑗 |𝜙 𝑗 (𝐸)| 𝑠 𝑗 • |𝜙 𝑠𝑑 (𝐸)| 𝑠 𝑠𝑑 ≤ 𝐾 𝜎 • 𝑁 𝑠 𝑠𝑑 𝑠𝑑
where 𝜎 = 𝑗 𝑠 𝑗 , not including 𝑠 𝑠𝑑 .

Note that the constraints on the 𝑠 𝑗 are relaxed due to the contribution of 𝑠 𝑠𝑑 , giving us an opportunity to find smaller values of 𝑠 𝑗 , compared to the usual method. Moreover, in order to tighten this upper bound, we should minimize 𝜎 first, then 𝑠 𝑠𝑑 next.

Using Reductions

A reduction is the successive application of an associative and commutative binary operator over a set of values. For example, in the convolution kernel described in Figure 3, we have a reduction over 3 dimensions (𝑐, ℎ and 𝑤) with the addition operation. We only consider reductions that sum all the values along some dimensions, called reduced dimensions. This property can be used to sum the values of a set in any order, instead of having a fixed sequence of summation.

Because the program is expressed as a loop nest, its reduction is sequential. Thus, when we examine the CDAG, we have a chain of dependencies linking the nodes along the reduced dimensions. When there are more than 2 reduced dimensions, it impacts negatively our study of the paths of dependencies of a CDAG, thus the quality of the extracted homomorphisms 𝜙 𝑗 and the upper bound found on the size of 𝐸.

We detect a simple class of reduction by performing a pattern-matching on the affine dependencies of the CDAG. We impose that the reduced dimensions are defined over a rectangular domain and that the summation is performed in a lexicographic order over a permutation of these reduced dimensions. Once a reduction is detected, we replace the sequential chain of dependencies (between the nodes of the reduction) by two sets of broadcast dependencies: (i) starting from the computations using the final result of the reduction to all the nodes of the reduction, and (ii) starting from the initialization of the reduction to the nodes of the reduction. ( ì 𝑓 )

1 ≤ 𝑠 1 + 𝑠 3 ( ì 𝑏) 1 ≤ 𝑠 1 + 𝑠 2 (ì 𝑐) 1 ≤ 𝑠 2 + 𝑠 3 ( ì 𝑥) 1 ≤ 𝑠 1 + 𝑠 2 ( ì 𝑤) 1 ≤ 𝑠 2 + 𝑠 3 +𝑠 𝑠𝑑 ( ì 𝑤 -ì 𝑥) 1 ≤ 𝑠 1 + 𝑠 3 +𝑠 𝑠𝑑 ( ì 𝑤, ì 𝑥) 2 ≤ 𝑠 1 + 𝑠 2 + 𝑠 3 +𝑠 𝑠𝑑 ( ì 𝑦) 1 ≤ 𝑠 1 + 𝑠 2 ( ì ℎ) 1 ≤ 𝑠 2 + 𝑠 3 +𝑠 𝑠𝑑 ( ì ℎ -ì 𝑦) 1 ≤ 𝑠 1 + 𝑠 3 +𝑠 𝑠𝑑 ( ì ℎ, ì 𝑦) 2 ≤ 𝑠 1 + 𝑠 2 + 𝑠 3 +𝑠 𝑠𝑑 (c) Constraints on 𝑠 𝑗 from Brascamp-Lieb inequality.
• No small dimensions:

𝑠 𝑗 = 2 3 , 𝜎 = 2. • Small dimensions: 𝑠 𝑗 = 1 2 , 𝑠 𝑠𝑑 = 1 2 , 𝜎 = 3 2 . (d) Solution obtained for the 𝑠 𝑗 , minimizing 𝜎 = 𝑗 𝑠 𝑗 .
Figure 3. Partitioning method -Brascamp-Lieb inequality applied to a convolution, without and with small dimensions (H and W).

On the convolution example of Figure 3, if we do detect the reduction, the dependence on Out is along the dimension 𝑠, and we will find as a corresponding homomorphism 𝜙 1 (𝑏, 𝑐, 𝑓 , 𝑥, 𝑦, ℎ, 𝑤) = (𝑏, 𝑐, 𝑓 , 𝑥, 𝑦, ℎ, 0). If we detect the reduction, the homomorphism of the new path becomes 𝜙 1 (𝑏, 𝑐, 𝑓 , 𝑥, 𝑦, ℎ, 𝑤) = (𝑏, 0, 𝑓 , 𝑥, 𝑦, 0, 0), which leads to a bigger kernel and better constraints on 𝑠 𝑗 .

Example -Convolution

To illustrate this method, we consider a convolution computation, described in Figure 3a.

Derivation with no small dimensions. This derivation corresponds to the one found in Demmel and Dihn [START_REF] Demmel | Communication-Optimal Convolutional Neural Nets[END_REF]. By examining the affine dependencies of the program, we find three homomorphisms 𝜙 1 , 𝜙 2 and 𝜙 3 , described in Figure 3b.

As stated by Theorem 5.1, we consider the subgroups H generated from the combination of kernel vectors of the 𝜙 𝑗 . We obtain the constraints shown in Figure 3c, ignoring the rightmost column. Their corresponding subgroups are shown in the leftmost column.

Finally, we solve this system of linear inequality to obtain a solution for the 𝑠 𝑗 , which corresponds to the upper bound on the size of a 𝐾-bounded set: |𝐸| ≤ 𝐾 2 .

Derivation with small dimensions. Now, we assume that the product of the problem size parameters 𝐻 and 𝑊 is small compared to the size of the small memory (𝑆). Compared to the first derivation, we have an additional homomorphism 𝜙 𝑠𝑑 , shown in Figure 3b. The presence of this new homomorphism adds a new coefficient 𝑠 𝑠𝑑 to the constraints on the 𝑠 𝑗 , as shown in Figure 3c. We obtain a better solution for the 𝑠 𝑗 , which corresponds to the upper bound on the size of a 𝐾-bounded set:

|𝐸| ≤ 𝐾 3 2 • (𝐻 • 𝑊 ) 1 2
. This is a tighter upper bound.

Experiments

Using our new I/O cost computation method, and our improvements to the lower bound algorithm for computing I/O complexity, IOOpt is able to provide tight bounds for both tensor contractions (even with small dimensions) and convolutions.

Benchmarks. We illustrate this by running it on representative benchmarks: For convolutions, we considered 11 different layers of Yolo9000 [START_REF] Redmon | YOLO9000: Better, Faster, Stronger[END_REF]. The parameter values for each layer are shown in Fig. 4. For tensor contractions, we considered the ones from TCCG [START_REF] Springer | Design of a High-Performance GEMM-like Tensor-Tensor Multiplication[END_REF].

The benchmark Python script from TCCG source code gathers 73 tensor contraction kernels (originating from various sources), that can be reduced to 49 different kernels, once synonyms are identified. For each kernel, TCCG selects the sizes of every dimension so that they are multiple of 8, roughly equal, and so that the product of all sizes is around 200 × 2 20 .

We can further reduce the number of relevant kernels that need to be considered in our analysis, by grouping them according to the number of dimensions of each array, and the number of dimensions shared between them. Indeed, the array layouts (dimension order) do not impact our analysis. This yields eight classes of tensor contraction kernels, described in Figure 5. Parametric lower bound expressions. Figure 6 shows the parametric lower bound found for the tensor contraction kernels and the 2D convolution kernels. For each line of this array, the first term of the maximum function is the sum of the volumes of the input arrays used by the computation, which is trivially a lower bound of the number of I/O. Because the lower bound derived is still sound, even if the small dimension hypothesis is not satisfied, we have to combine various lower bounds for different small dimension scenarios together.

Layer

For the tensor contraction kernels, we considered 8 = 2 3 small dimension scenarios. Dimensions shared between two arrays are grouped together, and every combination of small/regular dimensions for those three groups is examined. For convolution kernels, we considered 5 small dimension scenarios, based on array sizes: (i) no small parameters, (ii) 𝐻 and 𝑊 small parameters, (iii) 𝐻 , 𝑊 and 𝐵 small parameters, (iv) 𝐻 , 𝑊 , 𝑋 , 𝑌 and 𝐵 small parameters, and (v) 𝐶, 𝐻 , 𝑊 and 𝐵 small parameters. Among those bounds, only the first three scenarios lead to interesting bounds. They correspond to the last three rows in the expression. The previous version of the algorithm (without reduction management or small dimensions) fails to find an interesting bound, and returns the sum of array sizes (first row in the expression).

Symbolic upper bound expressions. The method presented in Sec. 4 provides a symbolic expression for the I/O of a program, as a function of program parameters and tile sizes. To compare them with the parametric lower bound expressions, we would like to remove tile sizes from the expression, and express them as functions of the cache size instead. In general, this is too complicated to solve. However, for the benchmarks we consider, this is possible by using only a few assumptions.

Let us start with the matrix multiplication example from Listing 1, with tiling (𝑖, 𝑗, 𝑘), 𝑇 𝑖 = Ti,𝑇 𝑗 = Tj,𝑇 𝑘 = 1 . The IO cost expression and the constraints on tile sizes are:

𝐼𝑂 = 𝐼𝑂 𝐴 + 𝐼𝑂 𝐵 + 𝐼𝑂 𝐶 = Ni • Nj • Nk • 1 Ti + 1 Tj + 1 Nk (1)
𝑆𝐷𝐹 𝐴,1 + 𝑆𝐷𝐹 𝐵,1 + 𝑆𝐷𝐹 𝐶,1 = Ti + Tj + Ti • Tj ≤ 𝑆 (2)
Then, we consider square tiles, by assuming that Ti and Tj are equal to the same value 𝑇 . We also assume that the tile completely fills the cache, which means that we consider that inequality (2) is actually an equality. While this may not be the best solution, these hypotheses will still provide a valid bound. We obtain the following equality:

2𝑇 + 𝑇 2 = 𝑆.
It has a unique positive solution:

𝑇 = √ 𝑆 + 1 -1.
We can then plug this value back into expression (1) to get our symbolic bound:

𝐼𝑂 = Ni • Nj • 2Nk √ 𝑆 + 1 -1 + 1
Note that the dominant term matches the dominant term of the I/O lower bound for matrix multiplication.

Here, we chose the tiling scheme manually, but this can be automated using TileOpt with some relevant numerical values of parameters. It will return a permutation and numerical tile sizes, and we can use the information on which dimensions have 𝑇 𝑑 = 1 or 𝑇 𝑑 = 𝑁 𝑑 to guide our symbolic solving. If the problem sizes are significantly bigger than √ 𝑆, the best tiling found by the solver will precisely be this one.

To generalize this reasoning to tensor contraction kernels, we only need to add some extra information on the expected order of magnitude of tile sizes.

Dimensions in a tensor contraction computation can be divided into three groups such that after "merging" the dimensions in each group, the computation is equivalent to a matrix multiplication. This corresponds to the "shared dimensions" in Fig. 5. The condition we impose is that the

+ 𝐵𝐶𝐹𝐺 LB = max 𝐴𝐵𝐶𝐷𝐸𝐹 + 𝐷𝐸𝐺𝐴 + 𝐺𝐹 𝐵𝐶, -2 + 𝐺 -𝑆 + 2𝐴𝐷𝐸 + 2𝐵𝐶𝐹 + 2𝐴𝐵𝐶𝐷𝐸𝐹 (𝐺-1) √ 𝑆 TC abc-adec-ebd UB = 2𝐴𝐵𝐶𝐷𝐸 √ 𝑆+1-1 + 𝐵𝐷𝐸 LB = max 𝐴𝐵𝐶 + 𝐴𝐷𝐸𝐶 + 𝐸𝐵𝐷, -3 + 𝐷𝐸 -𝑆 + 3𝐴𝐶 + 3𝐵 -𝐴𝐵𝐶 + 2𝐴𝐵𝐶 (𝐷𝐸-1) √ 𝑆 TC ab-cad-dcb UB = 2𝐴𝐵𝐶𝐷 √ 𝑆+1-1 + 𝐴𝐵 LB = max 𝐴𝐵 + 𝐶𝐴𝐷 + 𝐷𝐶𝐵, -3 + 𝐶𝐷 -𝑆 + 3𝐴 + 3𝐵 -𝐴𝐵 + 2𝐴𝐵 (𝐶𝐷-1) √ 𝑆 TC ab-ac-cb UB = 2𝐴𝐵𝐶 √ 𝑆+1-1 + 𝐵𝐶 LB = max 𝐴𝐵 + 𝐴𝐶 + 𝐶𝐵, -2 + 𝐶 -𝑆 + 2𝐴 + 2𝐵 + 2𝐴𝐵 (𝐶-1) √ 𝑆 TC abcd-aebf-fdec UB = 2𝐴𝐵𝐶𝐷𝐸𝐹 √ 𝑆+1-1 + 𝐶𝐷𝐸𝐹 LB = max 𝐴𝐵𝐶𝐷 + 𝐴𝐸𝐵𝐹 + 𝐹 𝐷𝐸𝐶, -3 + 𝐸𝐹 -𝑆 + 3𝐴𝐵 + 3𝐶𝐷 -𝐴𝐵𝐶𝐷 + 2𝐴𝐵𝐶𝐷 (𝐸𝐹 -1) √ 𝑆 2D Convolution UB = 𝐶𝐹 𝐻𝑊 𝑋𝑌 1 𝑋𝑌 + 1 𝐻 Δ𝑊 + (𝐻 +Δ-1) (𝑊 +𝑋 -1) 𝐻 Δ 2 𝑊 𝑋 where Δ = -𝐻𝑊 +𝑊 + √ 𝐻 2 𝑊 2 +4𝐻𝑆𝑊 -2𝐻𝑊 2 +4𝑆𝑊 +4𝑆+𝑊 2 2(𝐻𝑊 +𝑊 +1) LB = max 𝐵𝐶 (𝑌 + 𝐻 -1) (𝑋 + 𝑊 -1) + 𝐵𝐹𝑋𝑌 + 𝐹𝐶𝐻𝑊 , -2 -𝑆 + 𝐶 + 4𝐹 + 𝐵𝑌 + 𝐵𝑋 + 2𝐵𝑋𝑌 -2𝐵𝑋𝑌 𝐹 + 𝐵𝐹𝑋𝑌 (𝑊 𝐻𝐶-1) 𝑆 , -2 -𝑆 + 𝐶 + 4𝐹 + 𝐵𝑌 + 𝐵𝑋 + 2𝐵𝑋𝑌 -2𝐵𝑋 𝐹𝑌 + 2𝐵𝑋𝑌𝐶𝐹 √ 𝐻𝑊 √ 𝑆 -2𝐵𝑋𝑌 𝐹 √ 𝐻𝑊 𝑆 , -2 -𝑆 + 𝐶 + 4𝐹 + 𝐵𝑌 + 𝐵𝑋 + 2𝐵𝑋𝑌 -2𝐵𝑋 𝐹𝑌 + 2𝑋𝑌𝐶𝐹 √ 𝐵𝐻𝑊 √ 𝑆 -2𝑋𝑌 𝐹 √ 𝐵 √ 𝐻𝑊 𝑆 Figure 6
. Combined parametric I/O bounds of tensor contraction (TC) and 2D convolution kernels. 𝑆 is the small memory size and other uppercase letters are problem sizes (for TC kernels A. . . F are tensor dimensions, for convolution see Fig. 3a).

products of tile sizes inside each group are equal. For example abc-adec-ebd, the three groups are {𝑎, 𝑐}, {𝑏} and {𝑑, 𝑒}, and the condition is

𝑇 𝑎 𝑇 𝑐 = 𝑇 𝑏 = 𝑇 𝑑 𝑇 𝑒 .
The expressions are shown in Figure 6. They correspond to the "general case" where the input size is not the bottleneck, i.e., when parameters are sufficiently large compared to 𝑆.

The expression for convolution is quite complex, but an asymptotic analysis shows that the highest order term is

2𝐶𝐹𝑋𝑌 √ 𝐻𝑊 √ 𝑆
when 𝐶, 𝐹, 𝑋 and 𝑌 are sufficiently large, which matches the third term in the lower bound (the 𝐵 factor does not appear since it is equal to 1 in all our benchmarks).

Comparison of upper and lower bounds for different cache sizes. Figure 7 shows a comparison between lower and upper bounds for the considered benchmarks. For several cache sizes (from 16 kB to 4MB), the lower bound is computed by plugging the actual values in the symbolic expressions shown above, and the upper bound is obtained by running TileOpt on the optimization problems generated by IOUB.

As a sanity check, we confirm that the upper bound found is always above the lower bound, which was not trivial for the convolution bounds, and that both curves are decreasing or are constant along 𝑆. Both bounds are close to each other, ranging from at most a factor of 3 between them for Yolo9000-2 on the smallest value of S, to almost the same values for the ab-ac-cb TC kernel, which is also known as matrix multiplication. Moreover, the upper bound and the lower bound come closer for the largest values of 𝑆, which shows a match between the asymptotic dominant term of both bounds. When the cache size becomes very large, the dominant cost becomes the initial loading of the input data, Evaluation of the tiling recommendation. We evaluate the efficiency of the tiling recommendation we produce on the Yolo-9000 benchmark (Figure 4). The tiling recommendation has been obtained using a simple multi-level extension of the cache model where the objective function minimizes the weighted (with available measured inverse bandwidth) sum of data movements between the different cache levels. We fix the innermost dimension of the permutation in order to force vectorization on dimension 𝑓 . Note that fixing the innermost dimensions of our tiling recommendation does not impact the rest of our model.

In Figure 8, we compare the original code (no tiling), the tiled code based on our tiling recommendation, with the convolution from OneDNN [START_REF]OneAPI Deep Neural Network Library (oneDNN)[END_REF], a state-of-the-art library. We ran our experiment on an Intel i9-7940X Skylake-X machine with AVX-512 vector instructions, a 32kB L1 cache, a 1MB L2 cache and a shared 20MB L3 cache. Each code is run 200 times, and the median is taken.

As expected, our implementation does not outperform OneDNN (with the exception of the first layer), but is still reasonably fast. Indeed, many optimizations such as register tiling that expose instruction level parallelism without stressing too much register pressure, packing, versioning, etc., which are implemented by OneDNN, are not present in our "naive" implementation. Thus, our tiling strategy can be improved in order to obtain a good performance.

Kernel

No Tiling OneDNN Limitations. While the methodology presented here is very general on the theoretical side, practical implementation has some limitations. When iteration domains get complex, for instance when loop counters depend on one another, or loop bounds are more involved functions of parameters, cardinalities computed by Barvinok's algorithm get more complex, with several possible expressions depending on relations between parameters. This can greatly increase the size of the optimization problem, making it currently intractable. Another limit is the method to generate symbolic upper bounds: currently we rely on solving a polynomial equation, which is not generally doable when the degree exceeds four. This could be addressed by relaxing the problem of finding the precise expression of a tile size that completely fills the cache to only finding a size that does not exceed the cache capacity.

Related Work

I/O complexity lower bound. Most papers on computing lower bounds are handmade proofs for specific algorithms. While early works including the seminal paper of Hong and Kung [START_REF] Hong | I/O complexity: The Red-Blue Pebble Game[END_REF] focused on deriving asymptotic order lower bounds without bothering about the constant factor, several more recent contributions have finally derived tight (up to lower order terms) lower bounds with scaling constants for a few specific classes of algorithms [1-3, 13, 22, 23, 30, 31, 36]. The early work of Chris et al. [START_REF] Christ | Communication Lower Bounds and Optimal Algorithms for Programs that Reference Arrays -Part 1[END_REF], followed by that of Elango et al. [START_REF] Elango | On Characterizing the Data Access Complexity of Programs[END_REF], paved the way to algorithmic approaches for automatic derivation of I/O complexity lower bounds for the subclass of affine programs. The first tool able to actually automate the process and apply it to a full benchmark suite is IOLB from Olivry et al. [START_REF] Olivry | Automated Derivation of Parametric Data Movement Lower Bounds for Affine Programs[END_REF]. One important issue of IOLB with regard to neural networks is the looseness of the obtained bound for convolutions. This is due to its inability to exploit the presence of small dimensions, as opposed to the recent work of Demmel et al. [START_REF] Demmel | Communication-Optimal Convolutional Neural Nets[END_REF] which uses a similar trick to that of Elango et al. in [START_REF] Elango | On Characterizing the Data Access Complexity of Programs[END_REF] for deriving tight bounds for CNN. The main limitation of [START_REF] Demmel | Communication-Optimal Convolutional Neural Nets[END_REF][START_REF] Elango | On Characterizing the Data Access Complexity of Programs[END_REF] resides in the inability to provide scaling constants. The algorithm presented in Sec. 5 extends the work of Olivry et al. for exploiting the presence of small dimensions without losing its main advantage: being able to provide an exact (not restricted to asymptotic order inequality) bound with scaling constant. I/O complexity upper bound, cache miss prediction & tile size selection . Computing an I/O complexity upper bound for an algorithm is the most reasonable way to assess the tightness of a lower bound. While this computation is usually done by hand using ad hoc techniques specific to each studied algorithm [START_REF] Aggarwal | The Input/Output Complexity of Sorting and Related Problems[END_REF][START_REF] Demmel | Communication-Optimal Convolutional Neural Nets[END_REF][START_REF] Kwasniewski | Red-blue Pebbling Revisited: Near Optimal Parallel Matrix-Matrix Multiplication[END_REF][START_REF] Olivry | Automated Derivation of Parametric Data Movement Lower Bounds for Affine Programs[END_REF][START_REF] Ranjan | Upper and Lower I/O Bounds for Pebbling r-Pyramids[END_REF][START_REF] Michael | A Tight I/O Lower Bound for Matrix Multiplication[END_REF], Fauzia et al. [START_REF] Fauzia | Beyond Reuse Distance Analysis: Dynamic Analysis for Characterization of Data Locality Potential[END_REF] proposed a heuristic that directly reasons on the CDAG, which unfortunately does not scale to real programs. Finding an upper bound for a fixed architecture can also be viewed as finding an optimized program transformation that minimizes data movement costs, which also implies being able to evaluate this cost. Thus, restricting the analysis to affine programs and using the polyhedral framework appears to be appropriate for this problem. However, while the seminal scheduling algorithm from PluTo [START_REF] Bondhugula | A Practical Automatic Polyhedral Program Optimization System[END_REF] is able to expose tilable loops and generate tiled code, it can only handle fixed tile sizes. This is because parametric tiling is not an affine transformation: a tiled affine code with parametric tile sizes is no longer affine. A consequence is that existing polyhedral tools cannot be used to evaluate the I/O cost of such a code. These tools include PolyFeat [START_REF] Bao | Static and Dynamic Frequency Scaling on Multicore CPUs[END_REF], which computes an approximation of the number of capacity misses for arbitrary affine programs, as well as other algorithms [START_REF] Bao | Analytical modeling of cache behavior for affine programs[END_REF][START_REF] Chatterjee | Exact Analysis of the Cache Behavior of Nested Loops[END_REF][START_REF] Ghosh | Cache Miss Equations: a Compiler Framework for Analyzing and Tuning Memory Behavior[END_REF] which focus on precisely conflict misses. Some of these algorithms are restricted to a small class of programs, and most of them can only model a one-level cache. They all need to consider fixed parameters and fixed tile sizes. Other works [START_REF] Li | Analytical Cache Modeling and Tilesize Optimization for Tensor Contractions[END_REF][START_REF] Shirako | Analytical Bounds for Optimal Tile Size Selection[END_REF][START_REF] Yuki | Automatic Creation of Tile Size Selection Models[END_REF] implemented ad hoc computation of this cost function for very restricted sub-classes of programs, but our algorithm is the first to be able to automatically generate a symbolic expression of it for arbitrary parametric tiled affine programs, in a multi-level cache setting.

We showed how we use this cost function to find the best loop permutations and tile sizes for a given architecture, using operations research. Once again, this is out of the scope of polyhedral analysis, as cost functions are not affine. So the usual optimization strategies used by polyhedral compilers [START_REF] Bondhugula | A Practical Automatic Polyhedral Program Optimization System[END_REF][START_REF] Grosser | Polly -Performing Polyhedral Optimizations on a Low-Level Intermediate Representation[END_REF][START_REF] Verdoolaege | Polyhedral Parallel Code Generation for CUDA[END_REF], mostly based on parametric integer programming [START_REF] Feautrier | Parametric Integer Programming[END_REF], cannot be used. Some literature exists on tile size selection, but it mostly deals with performance model design and is either handcrafted for specific kernels and access function patterns [START_REF] Liu | Tile Size Selection for Optimized Memory Reuse in High-level Synthesis[END_REF][START_REF] Prajapati | Simple, Accurate, Analytical Time Modeling and Optimal Tile Size Selection for GPGPU Stencils[END_REF] or uses machine learning [START_REF] Yuki | Automatic Creation of Tile Size Selection Models[END_REF] and even cache simulators [START_REF] Mehta | Tile Size Selection Revisited[END_REF]. Polyhedral analysis is only used for the tiling transformation. Renganarayana and Rajopadhye [START_REF] Renganarayanan | Positivity, Posynomials and Tile Size Selection[END_REF] showed that most performance models for tile size selection used in the literature are polynomials that are operations research-friendly. Our automatically generated cost function, which matches the distinct-access model promoted by Ferrante et al. [START_REF] Ferrante | On Estimating and Enhancing Cache Effectiveness[END_REF], fits into this category.

Conclusion

We have presented an algorithm which computes an upper bound on the data movement complexity of a polyhedral program. It outputs a symbolic formula on the problem sizes and the cache size, and advises a loop permutation and tiling which minimizes it. Additionally, we have improved an algorithm which computes the lower bound of the data movement complexity, by exploiting the fact that some parameters are much smaller than the cache size. Both algorithms were implemented as a fully automatic tool and evaluated on multiple convolution and tensor contraction kernels. We have evaluated the numerical tightness of these bounds by evaluating them for the problem sizes of tensor contraction and convolution benchmarks.

Figure 1 .

 1 Figure 1. Organisation of the different components of IOOpt

  and access function 𝑓 𝐴 (𝑖, 𝑗, 𝑘) = (𝑖, 𝑘); • for array 𝐵, with domain M 𝐵 = {𝑥, 𝑦 | 0 ≤ 𝑥 < 𝑁 𝑘 ∧ 0 ≤ 𝑦 < 𝑁 𝑗 } and access function 𝑓 𝐵 (𝑖, 𝑗, 𝑘) = (𝑘, 𝑗); • and for array 𝐶, with domain M 𝐶 = {𝑥, 𝑦 | 0 ≤ 𝑥 < 𝑁 𝑖 ∧ 0 ≤ 𝑦 < 𝑁 𝑗 } and access function 𝑓 𝐶 (𝑖, 𝑗, 𝑘) = (𝑖, 𝑗).

  Tf × Nw.Nc.Nf + 1 Nw.Tf × Nw.Nc.Nf.(Nx -1) = Nc.Nf.(Nx + Nw -1) Tf The same computation for arrays Out and Filter yields (assuming the outermost reuse dimensions are respectively 𝑥 and 𝑓 ): 𝐼𝑂 Out = Nc.Nf.Nx/Tc 𝐼𝑂 Filter = Nc.Nf.Nw

4 if 6 P := ∅; 7 forall 9 𝑅 10 P 11 PFigure 2 .

 467910112 Figure 2. Application of Algorithm 1 to the 1D convolution kernel

  for (b = 0; b < B; b++) for (c = 0; c < C; c++) for (f = 0; f < F; f++) for (x = 0; x < X; x++) for (y = 0; y < Y; y++) for (h = 0; h < H; h++) for (w = 0; w < W; w++) { Out[f,x,y,b] += Image[x+h,y+w,c,b] * Filter[f,h,w,c]; } (a) Convolution kernel. 𝜙 1 (𝑏, 𝑐, 𝑓 , 𝑥, 𝑦, ℎ, 𝑤) = (𝑏, 0, 𝑓 , 𝑥, 𝑦, 0, 0) (from reduction) 𝜙 2 (𝑏, 𝑐, 𝑓 , 𝑥, 𝑦, ℎ, 𝑤) = (𝑥 + ℎ, 𝑦 + 𝑤, 𝑐, 𝑏) (from Image) 𝜙 3 (𝑏, 𝑐, 𝑓 , 𝑥, 𝑦, ℎ, 𝑤) = (𝑓 , ℎ, 𝑤, 𝑐) (from Filter) 𝜙 𝑠𝑑 (𝑏, 𝑐, 𝑓 , 𝑥, 𝑦, ℎ, 𝑤) = (ℎ, 𝑤) (small dimensions) (b) Homomorphisms used in Brascamp-Lieb inequality.

Figure 5 .

 5 Figure 5. Classes of Tensor Contraction kernels from the TCCG benchmarks. The classes are determined from the number of dimensions of the arrays (Out / In1 / In2), and the number of shared dimensions (s. d.) between arrays (Out+In1 / Out+In2 / In1+In2)

1 + 1 +√ 𝑆+1- 1 +

 111 KernelI/O bounds: lower (LB) and upper (UB)TC abcde-efbad-cf UB = 2𝐴𝐵𝐶𝐷𝐸𝐹 √ 𝑆+1-𝐶𝐹 LB = max 𝐴𝐵𝐶𝐷𝐸 + 𝐸𝐹 𝐵𝐴𝐷 + 𝐶𝐹, -2 + 𝐹 -𝑆 + 2𝐶 + 2𝐴𝐵𝐷𝐸 + 2𝐴𝐵𝐶𝐷𝐸 (𝐹 -1) √ 𝑆 TC abcd-dbea-ec UB = 2𝐴𝐵𝐶𝐷𝐸 √ 𝑆+1-𝐶𝐸 LB = max 𝐴𝐵𝐶𝐷 + 𝐷𝐵𝐸𝐴 + 𝐸𝐶, -2 + 𝐸 -𝑆 + 2𝐶 + 2𝐴𝐵𝐷 + 2𝐴𝐵𝐶𝐷 (𝐸-1) √ 𝑆 TC abc-bda-dc UB = 2𝐴𝐵𝐶𝐷 𝐶𝐷 LB = max 𝐴𝐵𝐶 + 𝐵𝐷𝐴 + 𝐷𝐶, -2 + 𝐷 -𝑆 + 2𝐶 + 2𝐴𝐵 + 2𝐴𝐵𝐶 (𝐷-1) √ 𝑆 TC abcdef-dega-gfbc UB = 2𝐴𝐵𝐶𝐷𝐸𝐹𝐺 √ 𝑆+1-1

Figure 7 .

 7 Figure 7. I/O bounds (lower is blue and upper is orange) of tensor contractions (TCCG) and convolutions (Yolo) for different cache sizes

Grenoble Alpes (ANR-19-P3IA-0003), and by the Bpifrance Programme d'Investissements d'Avenir (PIA) as part of the ES3CAP project.


Dalil Ichalal
email: dalil.ichalal@univ-evry.fr

Saïd Mammar

State and fault estimation of a class of nonlinear systems with slow internal dynamics: Asymptotic Decoupling approach

Keywords: State estimation, Unknown Input Observers (UIO), Asymptotic Decoupling concept, Internal Dynamics, Nonlinear Lipschitz Systems, LMIs

This paper presents an Unknown Input Observer (UIO) for a class of Lipschitz nonlinear systems affected by disturbances/faults and having slow internal dynamics after the decoupling of the Unknown Input (UI). Thanks to the Asymptotic Decoupling notion, the constraint related to the internal dynamics is bypassed. The nonlinear UIO presents time varying dynamics that allows to obtain the Asymptotic Decoupling notion. It will be shown that the proposed approach allows to enhance the convergence rate of the UIO. The convergence enhancement of the state estimation error is demonstrated via a Lyapunov analysis and the established conditions are expressed in terms of Linear Matrix Inequalities. Finally, simulations will be given in order to illustrate the proposed approach with some comparisons to existing techniques.

INTRODUCTION

Fault diagnosis and state estimation are important tasks in modern automatic control design. It allows to estimate precisely the state and detect a fault earlier and exploit these informations for control, supervision and maintenance, and consequently, enhance the performances of the system and the security of human operator and the system.

Model-based approaches have reached a certain maturity these last years and several diagnosis problems have been solved thanks to such approaches [START_REF] Ding | Model-Based Fault Diagnosis Techniques Design Schemes, Algorithms, And Tools[END_REF]; [START_REF] Gertler | Fault detection and diagnosis in engineering systems[END_REF]; [START_REF] Patton | Fuzzy observers for non-linear dynamic systems fault diagnosis[END_REF]. The main commonly used approach is the one exploiting the observers as residual generators through mild modifications. For example, the design purpose differs from the observer design in the sense that the residual generators are designed in such a way to maximize the influence of faults on the residual and to minimize the effect of eventual disturbances. This concept can be found in many approaches using the H ∞ and H -techniques [START_REF] Mazars | Computation of a reference model for robust fault detection and isolation residual generation[END_REF][START_REF] Mazars | A QMI approach to the robust fault detection and isolation problem[END_REF]. In [START_REF] Jaimoukha | A matrix factorization solution to the H -/H ∞ fault detection problem[END_REF], a factorization approach is provided to solve the H -/H ∞ fault detection problem. The same problem is explored in [START_REF] Liu | An {LMI} approach to minimum sensitivity analysis with application to fault detection[END_REF], [START_REF] Wang | An {LMI} approach to Hindex and mixed H -/H ∞ f ault detection observer design[END_REF] by using the sensitivity analysis. In [START_REF] Zhong | An LMI approach to design robust fault detection filter for uncertain LTI systems[END_REF], a LMI approach is provided for fault detection by introducing a reference model in order to shape the residual signal. This approach is developed for uncertain LTI systems affected by faults on both the state and the output equations. The problem of actuator fault diagnosis is generally more complex to study compared to sensor ones due the fact that the actuator faults do not affect directly the outputs. These techniques can be enhanced by a simple analysis of the system instead of applying the H -approach directly. In another hand, more interesting approaches based on Proportional Integral Observers can be exploited in order to estimate both the state of the system and the fault affecting it. This technique is suitable when a knowledge on the dynamic of the system is available. This is the case for incipient or abrupt fault which can be approximated by second order polynomials of the form f (t) = a 1 t + a 0 which satisfies the property f (t) = 0. Of course, it is possible to derive similar results for more general models of the faults (see [START_REF] Gao | Unknown input observer-based robust fault estimation for systems corrupted by partially decoupled disturbances[END_REF]; [START_REF] Ichalal | Simultaneous state and unknown inputs estimation with PI and PMI observers for Takagi-Sugeno model with unmeasurable premise variables[END_REF] for more details). There exists also other techniques based on Unknown Input Observers which are suitable for system where no information is known on the faults. In such a situation, the decoupling approach is more suitable instead of simultaneous estimation of the state and the faults by PI observers. However, the use of UIO requires two main conditions [START_REF] Darouach | Full-order observers for linear systems with unknown inputs[END_REF]: The decoupling condition, and the stability of the internal dynamics. Since the internal dynamics cannot be moved or controlled in order to enhance the performances of the observer the dynamics of the UIO is constrained by such dynamics. Recently, in [START_REF] Ichalal | On unknown input observers of linear systems: Asymptotic unknown input decoupling approach[END_REF], a solution is proposed in order to handle such constraint. It consists on the use of the notion "Asymptotic Decoupling" instead of the classical "Algebraic Decoupling". This notion allows to enhance the convergence rate of the UIO in the transient phase.

In this paper, the Asymptotic Decoupling approach is used to construct an UIO in order to estimate both the state and the faults, while guaranteeing a good convergence rate even in the presence of stable but slow internal dynamics. It will be shown that this fact allows to estimate correctly the UI compared to the classical approaches.

PROBLEM STATEMENT

Let us consider the nonlinear systems of the form

ẋ(t) = Ax(t) + Ed(t) + f (x(t), u(t)) y(t) = Cx(t) (1)
where

x ∈ R n , u ∈ R n u , d ∈ R n d
f (x, u) -f (y, u) ≤ ν x -y , ∀x ∈ D ⊆ R n
(2) The UI d(t) may model disturbances and / or faults affecting the state equation. In order to derive the proposed UIO, the following assumptions are also assumed Assumption 1. (1) The internal dynamics are stable

(2) The unknown input

d(t) is bounded d(t) ∞ ≤ σ , σ > 0 (3) The system without UI is observable (4) rank(CE) = rank(E)
The last hypothesis is commonly used in LTI UIO design in order to ensure that there exists a constant matrix that is able to decouple the Unknown Inputs (UI) (See for example [START_REF] Darouach | Full-order observers for linear systems with unknown inputs[END_REF]). This paper focuses on this class of systems but it is possible to generalize the approach to systems that do not satisfy this rank condition (see for example [START_REF] Ichalal | On unknown input observers for lpv systems[END_REF], Ichalal et al. (2015), [START_REF] Barbot | An algebraic framework for the design of nonlinear observers with unknown inputs[END_REF], [START_REF] Floquet | On sliding mode observers for systems with unknown inputs[END_REF] for linear, LPV or nonlinear systems). In the rest of the paper it is assumed that the effect of the UI on x is taken into account in the subset D ensuring (2).

UNKNOWN INPUT OBSERVER DESIGN

Structure of the UIO

The proposed UIO takes the following form

             ż(t) = N β z(t) + G β f ( x(t), u(t)) + L β y(t) x(t) = z(t) -M β y(t) M β = (1 -β (t))H H = -E(CE) - β (t) = κ exp(-αt) η(t) = ν(1 + β (t)) (3)
where (.) -denotes the pseudo-inverse of a matrix and κ, α and ν are positive parameters to design and η(t) is a positive bounded function satisfying the following definitions [START_REF] Ichalal | On unknown input observers of linear systems: Asymptotic unknown input decoupling approach[END_REF]. The notation X β means that the matrix X depend on the time-varying parameter β (t). The structure of the matrices will be defined later.

The following definition presents a class of functions which will be used a a time-varying decay rate [START_REF] Ichalal | On unknown input observers of linear systems: Asymptotic unknown input decoupling approach[END_REF]

(2) η 1 ≤ η(t) ≤ η 2 , with 0 < η 1 < η 2 Notice that lim t→+∞ t 0 η(τ)dτ = +∞
The objective is to design the matrices N β , G β , L β and H in order to ensure asymptotic state estimation and to enhance the convergence rate even in the presence of stable but slow internal dynamics.

The state estimation error is given by e

(t) = x(t) -x(t) = (I + M β C)x(t) -z(t) (4
) By defining the matrix P β = I + M β C and computing the time derivative of e(t), it follows

ė(t) = P β A + Ṗβ -L β C -N β P β x(t) + P β Ed(t) + N β e(t) + P β f (x(t), u(t)) -G β f ( x(t), u(t))
(5)

Under the conditions

P β A + Ṗβ -L β C -N β P β = 0 (6) P β = G β (7)
The state estimation error dynamics becomes

ė(t) = N β e(t) + P β f (x(t), x, u(t)) + S β d(t) (8) where f (x(t), x, u(t)) = f (x(t), u(t)) -f ( x, u(t))
and S β = P β E.

In one hand, as discussed in [START_REF] Zemouche | Observers for a class of Lipschitz systems with extension to h ∞ performance analysis[END_REF]; Phanomchoeng et al. ( 2011); [START_REF] Ichalal | Decoupling unknown input observer for nonlinear quasi-lpv systems[END_REF], instead of using the classical Lipschitz condition, it is much more interesting to use the polytopic form of the Lipschitz condition with the aid of the Differential Mean Value Theorem. Therefore, the term f (x(t), x(t), u(t)) can be expressed in the following polytopic form

f (x(t), x(t), u(t)) = r ∑ i=1 h i (z(t)) Āi e(t) = Φ h e(t) (9)
with z ∈] min(x, x), max(x, x)[. From this equivalent expression, the state estimation error dynamics is written equivalently as follows

ė(t) = N h,β e(t) + S β d(t) (10)
where

N h,β = N β + P β Φ h .
On the other hand, due to the Asymptotic decoupling concept, the matrix S β converges to zero asymptotically i.e. lim t→∞ S β = 0.

Consequently, if the matrix N h,β is stable, the state estimation error e(t) converges asymptotically to zero.

The matrix S β can be expanded as follows

S β = P β E = E + (1 -β (t))HCE
And by construction, the matrix H = -E(CE) -satisfies (I + HC)E = 0, it comes

S β = -β (t)HCE
And under the assumptions that E 2 ≤ λ ) and β (t) is positive, one obtains the following norm bound

S β 2 ≤ λ β (t) (11)
Consequently, due to the definition of the function β (t), the norm S β 2 converges to zero when t → ∞. Then, if the state estimation error dynamics is stable, it will be asymptotically converging to zero since S β goes to zero too.

Stability analysis

Let us consider the Lyapunov function candidate

V (e(t)) = e T (t)Xe(t), X = X T > 0 (12) Its time derivative along the trajectory of ( 8) is expressed by

V (e(t)) = e T (t) N T h,β X + XN h,β e(t) + e T (t)XS β d(t) + d T (t)S T β Xe(t) (13)
By adding and subtracting the term -η(t)e T (t)Xe(t) + γd T (t)S T β S β d(t) where η(t) is a positive bounded and monotonic decreasing scalar function (0 < η 1 ≤ η(t) ≤ η 2 ) and γ is a positive scalar (for more details of the function η(t) please refer to [START_REF] Ichalal | On unknown input observers of linear systems: Asymptotic unknown input decoupling approach[END_REF]). Therefore, the inequality ( 13) is equivalent to

V (e(t)) = ζ T (t)Ξ h,β ζ (t) -η(t)V (e(t)) + γd T (t)S T β S β d(t) (14)
where

Ξ h,β = N T h,β X + XN h,β + η(t)X XS β S T β X -γS T β S β , ζ (t) = e(t) d(t) (15)
Finally, if Ξ h,β ≤ 0 and taking into account (11) the time derivative of the Lyapunov function can be bounded as follows V (e(t)) ≤ -η(t)V (e(t)) + γλ 2 β 2 (t) d(t) 2 2 (16) The solution of the inequality ( 16) is given by V (e(t)) ≤ V (e(0)) + γλ 2 σ 2 ϖ e This proves that V (e(t)) converges asymptotically to zero with the time varying rate η(t).

LMI design procedure

In the previous section, the stability analysis is given and it is established under the condition that Ξ h,β ≤ 0. Straightforwardly, this condition can be expressed in a polytopic form. Firstly, let us recall the matrices N h,β = P β (A + Φ h ) + Ṗβ -K β C (18) which comes from the equation ( 6) where

K β = L β + N β M β and P β = I + M β C.
The matrix Ṗβ is expressed as follows

P β = I + M β C = I + (1 -κ exp(-αt
))HC Consequently, its time derivative is given by Ṗβ = κα exp(-αt))HC Then, the matrix N h,β can be expressed as follows

N h,β = (A + Φ h ) + (1 -κ exp(-αt))HC(A + Φ h ) + κα exp(-αt)HC -K β C (19) or in compact form as follows N h,β = A h,β -K β C where A h,β = (I + (1 -β (t))HC)(A + Φ h ) + αβ (t)HC Knowing that the function 0 < β (t) ≤ κ, a first polytopic form is obtained as follows          A h,β = 2 ∑ j=1 µ j (β (t)) Ā h, j S β = 2 ∑ j=1 µ j (β (t))S j (20)
where

Ā h,1 = (I + (1 -κ)HC)(A + Φ h ) + ακHC (21) Ā h,2 = (I + HC)(A + Φ h ) (22) S 1 = -κHCE (23) S 2 = 0 (24)
From this polytopic form, we construct the gain matrix K β and the positive function η(t) in the same polytopic form as follows

K β = 2 ∑ j=1 µ j (β (t))K j , η(t) = 2 ∑ j=1 µ j (β (t))η j
which allows to get the matrix N h,β in the form

N h,β = 2 ∑ j=1 µ j (β (t)) Ā h, j -K j C = 2 ∑ j=1 µ j (β (t)) N h, j
Consequently, the matrix Ξ h,β can be expressed as follows

Ξ h,β = 2 ∑ j=1 2 ∑ k=1 µ j (β )µ k (β )Θ h, j,k
where :

Θ h, j,k = N T h, j X + X N h, j + η j X XS j S T j X
-γS T j S k Now, by using the polytopic form of the matrix Φ h , it allows to obtain

Ξ h,β = r ∑ i=1 2 ∑ j=1 2 ∑ k=1 h i (z(t))µ j (β )µ k (β )Θ i, j,k where Θ i, j,k = N T i, j X + X N i, j + η j X XS j S T j X -γS T j S k (25)
Then, sufficient conditions for the negativity of Ξ h,β are given by Θ i, j,k ≤ 0, i = 1, ..., r, j, k = 1, 2

Finally, change of variables K j = XK j is performed in order to obtain LMI conditions (See Ichalal and Mammar (2020) for more details).

Algorithm

The UIO design procedure is summarized as follows:

(1) Define the parameters κ and α (2) Compute the matrix H = -E(CE) -and compute the matrix M β (3) Compute the polytopic form of the nonlinear term (9) (4) Construct the polytopic matrices (20) and then the matrices (25) (5) Solve the LMIs

X = X T > 0 (26) N T i, j X + X N i, j + η j X XS j S T j X -γS T j S k ≤ 0 (27) i = 1, ..., r j, k = 1, 2
and obtain the gain matrix K β (6) Compute the remaining matrices:

N h,β = P β (A + Φ h ) + Ṗβ -K β C L β = K β -N h,β H

Unknown Input Estimation

After estimating the state x(t) of the system, it is often important to estimate the Unknown Inputs. Indeed, in the fault diagnosis problems or fault tolerant controllers, the availability of a precise estimation of the faults is highly required for better performances in diagnosis or fault tolerant control tasks. In order to obtain such an estimation, the following expression is used

d(t) = (CE) -( ẏ -CA x(t) -C f ( x(t), u(t)))
which is nothing than the inversion of the system, but since the state estimation x(t) is precise, the estimation of the unknown input is also precise (see the simulation example). This estimation can be used for both fault diagnosis or in fault tolerant controllers to compensate the effect of the faults. Notice that due to the inherent presence of noises in the measurements y, the use of recent signal differentiation algorithms are suitable, instead of classical derivation, in order to have a precise estimation of the time derivative of y. For instance, it is possible to use the High Order Sliding Mode Differentiator [START_REF] Levant | Quasi-continuous high-order sliding-mode controllers[END_REF] or ALIEN differnetiator [START_REF] Barbot | An algebraic framework for the design of nonlinear observers with unknown inputs[END_REF] or Algebraic asymptotic differentiator [START_REF] Ibrir | Online exact differentiation and notion of asymptotic algebraic observers[END_REF].

Another application of the proposed observer is in the classical bank of observers for actuator fault detection and isolation in the case where the number of inputs is greater than the number of outputs.

SIMULATION RESULTS AND COMPARISONS

Let us consider the system (1) defined by the following matrices

A = 0 1 0 -3 -2 1 -5 -6 -0.03 , E = 0 1 0 , C = [ 1 1 0 ] , f (x, u) = [ 0 u(t) x 2 sin(x 2 ) ]
The input is u(t) = sin(t) and the fault signal is defined by

d(t) =      0 0 ≤ t < 5 4 5 ≤ t < 15 t -11 15 ≤ t < 20 t -11 + 2 sin(1.5t) cos(0.5t) t ≥ 20
In order to construct the UIO, the time varying function β (t) = 10 exp(-0.1t) (κ = 10 and α = 0.1) is considered and the matrix H is computed from H = -E(CE) -. By following the procedure of constructing the matrices of the UIO, the proposed LMIs provide a solution and the obtained results are depicted in the figures 1 and 2. Notice that, for comparison, the classical UIO is constructed according to the procedure given in [START_REF] Chen | Unknown input observer design for a class of nonlinear systems: an LMI approach[END_REF] with a Lipschitz constant ν = 10 in (2) computed for a subset D = {x ∈ R 3 : |x 2 | ≤ 10}. This same subset is used also, in the proposed approach, to express the polytopic form of the nonlinear term f (x, u) as in (9). From the figure 1, it can be seen that the slow internal dynamics affect the convergence rate of the observer, especially, for the third state which is only detectable with the classical Exact Decoupling approach. However, with the proposed Asymptotic Decoupling notion, the state estimation is better. Additionally, since the state estimations are used to compute an estimation of the fault, these lasts will be considerably affected if the estimated states are obtained from a classical UIO which is not the case with the proposed UIO as illustrated in figure 2. This paper presents a novel approach to design an observerbased fault diagnosis by using a new UIO. The particularity of the proposed UIO is its ability to handle stable but slow internal dynamics of a system and then enhance the convergence rate of the observer. The proposed UIO uses a time varying function aiming to ensure the asymptotic decoupling property. This property allows to preserve the observability (lost due exact decoupling) at least in a short time interval in the transient phase. Consequently, this allows to enhance the convergence rate of the state estimation error dynamics, and then provide an acceptable estimation of the system states and faults.

  s)ds β 2 (τ)dτ < ϖ, ∀t

Fig. 1 .Fig. 2 .

 12 Fig. 1. Real states (blue) and estimations (red). Classical UIO (Left), Proposed UIO (Right)

  and y ∈ R n y represent, respectively, state, known input, unknown input and output vectors with corresponding dimensions. A, E and C are real matrices with appropriate dimensions. f (x, u) is a Lipschitz nonlinear vector function depending on x and u which satisfies the condition