Dmitrij Celov

Remigijus Leipus

Anne Philippe

Asymptotic normality of the mixture density estimator in a disaggregation scheme *

Keywords: random coefficient AR(1), long memory, aggregation, disaggregation, mixture density

. The moving average representation of the aggregated process is investigated. A small simulation study illustrates the result.

Introduction

Aggregated time series data appears in different fields of studies including applied problems in hydrology, sociology, statistics, economics. Considering aggregation as a time series object, a number of important questions arise. These comprise the properties of macro level data obtained by small and large-scale aggregation in time, space or both, assumptions of when and how the inverse (disaggregation) problem can be solved, finally, how to apply theoretical results in practice.

Aggregated time series, in fact, can be viewed as a transformation of the underlying time series by some (either linear or non-linear) specific function defined at (in)finite set of individual processes. In this paper we consider a linear aggregation scheme, which is natural in applications. In practice it is found convenient to approximate individual data by simple time series models, such as AR(1), GARCH(1, 1) for instance (see [START_REF] Lewbel | Aggregation and simple dynamics[END_REF], [START_REF] Chong | The polynomial aggregated AR(1) model[END_REF], [START_REF] Zaffaroni | Contemporaneous aggregation of linear dynamic models in large economies[END_REF][START_REF] Zaffaroni | Contemporaneous aggregation of GARCH processes[END_REF]), whereas more complex individual data models do not provide an advantage in accuracy and efficiency of estimates, and usually are very difficult to study from the theoretical point of view.

Aggregation by appropriately averaging the micro level time series models can give intriguing results. It was shown in Granger (1980) that the large-scale aggregation of infinitely many short memory AR(1) models with random coefficients can lead to a long memory fractionally integrated process. It means that the properties of an aggregate time series may in general differ from those of individual data.

It is clear however that the weakest point of the aggregation is a considerable loss of information about individual characteristics of the underlying data. Roughly speaking, an aggregated time series can not be as informative about the attributes of individual data as the micro level processes are. On the other hand, using some special aggregation schemes, which involve, for instance, independent identically distributed "elementary" processes with known structure (such as AR(1)), enables to solve an inverse problem: to recover the properties of individual series with the aggregated data at hand. This problem is called a disaggregation problem.

Different aspects of this problem were investigated in Dacunha-Castelle, Oppenheim (2001), [START_REF] Leipus | Orthogonal series density estimation in a disaggregation scheme[END_REF], [START_REF] Celov | Time series aggregation, disaggregation and long memory[END_REF]. The last two papers deal with the asymptotic statistical theory in the disaggregation problem: they present the construction of the mixture density estimate of the individual AR(1) models, the consistency of an estimate, and some theoretical tools needed here. Resuming the previous research, the major objective of the present paper is to obtain the asymptotic normality property of the mixture density estimate, that enlarges the range of applications, solving the accuracy of simulation studies, statistical inference, forecasting and other problems.

Section 2 describes the disaggregation scheme, including the construction of mixture density estimate proposed by [START_REF] Leipus | Orthogonal series density estimation in a disaggregation scheme[END_REF], and formulates the main result of the paper. Important issues about the moving average representation of the aggregated process are discussed in Section 3. The proof of the main theorem and auxiliary results are given respectively in Section 4 and Section 7. Some simulation results are presented in Section 5.

Preliminaries and the main result

Consider a sequence of independent processes Y (j) = {Y (j) t , t ∈ Z}, j ≥ 1 defined by the random coefficient AR(1) dynamics

Y (j) t = a (j) Y (j) t-1 + ε (j) t , (2.1)
where ε (j) t , t ∈ Z, j = 1, 2, . . . are independent identically distributed (i.i.d.) random variables with Eε j) , j = 1, 2, . . . are i.i.d. random variables with |a| ≤ 1 and satisfying

(j) t = 0 and 0 < σ 2 ε = E(ε (j) t ) 2 < ∞; a, a (
E 1 1 -a 2 < ∞. (2.2)
It is assumed that the sequences {ε

(j)
t , t ∈ Z}, j = 1, 2, . . . and {a, a (j) , j = 1, 2, . . . } are independent.

Under these conditions, (2.1) admits a stationary solution Y (j) and, according to [START_REF] Oppenheim | Aggregation of random parameters Ornstein-Uhlenbeck or AR processes: some convergence results[END_REF], the finite dimensional distributions of the process

X (N ) t = 1 √ N N j=1 Y (j) t , t ∈ Z,
weakly converge as N → ∞ to those of a zero mean stationary Gaussian process X = {X t , t ∈ Z}, called the aggregated process. Suppose that random coefficient a admits a density ϕ(x), absolutely continuous with respect to the Lebesgue measure, which by (2.2) satisfies

1 -1 ϕ(x) 1 -x 2 dx < ∞. (2.3)
Any density function satisfying (2.3) will be called a mixture density.

Note that the covariance function and the spectral density of aggregated process X coincides with those of Y (j) and are given, respectively, by

σ(h) := Cov(X h , X 0 ) = σ 2 ε 1 -1 x |h| 1 -x 2 ϕ(x)dx (2.4) and f (λ) = σ 2 ε 2π 1 -1 ϕ(x) |1 -xe iλ | 2 dx.
(2.5)

The disaggregation problem deals with finding the individual processes (if they exist) of form (2.1), which produce the aggregated process X with given spectral density f (λ) (or covariance σ(h)). This is equivalent to finding ϕ(x) such that (2.5) (or (2.4)) and (2.3) hold. In this case, we say that the mixture density ϕ(x) is associated with the spectral density f (λ).

In order to estimate the mixture density ϕ(x) using aggregated observations X 1 , . . . , X n , [START_REF] Leipus | Orthogonal series density estimation in a disaggregation scheme[END_REF] proposed the estimate based on a decomposition of function

ζ(x) = ϕ(x)(1 -x 2 ) -α in the orthonormal L 2 (w (α) )-basis of Gegenbauer polynomials {G (α) k (x), k = 0, 1, . . . }, where w (α) (x) = (1 -x 2 ) α , α > -1. This decomposition is valid (i.e. ζ belongs to L 2 (w (α) )) if 1 -1 ϕ 2 (x) (1 -x 2 ) α dx < ∞, α > -1. (2.6) Let G (α) n (x) = k j=0 g (α)
n,j x j . The resulting estimate has the form k,j (σ n (j) -σn (j + 2)), (2.8) σ2 n,ε = σn (0) -σn (2) is the consistent estimator of variance σ 2 ε and σn (j) = n -1 n-j i=1 X i X i+j is the sample covariance of the aggregated process. Truncation level K n satisfies

φn (x) = σ-2 n,ε (1 -x 2 ) α Kn k=0 ζn,k G (α) k (x), ( 2
K n = [γ log n], 0 < γ < (2 log(1 + √ 2)) -1 .
(2.9) [START_REF] Leipus | Orthogonal series density estimation in a disaggregation scheme[END_REF] assumed the following semiparametric form of the mixture density:

ϕ(x) = (1 -x) 1-2d 1 (1 + x) 1-2d 2 ψ(x), 0 < d 1 , d 2 < 1/2, ( 2.10)
where ψ(x) is continuous on [-1, 1] and does not vanish at ±1. Then, under conditions above and corresponding relations between α and d 1 , d 2 , they showed the consistency of the estimator φn (x) assuming that the variance of the noise, σ 2 ε , is known and equals 1. In more realistic situation of unknown σ 2 ε , it must be consistently estimated. In order to understand intuitively the construction of estimator σ2

n,ε , it suffices to note that, by (2.4), σ 2 ε = σ(0) -σ(2). Also note that the estimator φn (x) in (2.7) possesses property 1 -1 φn (x)dx = 1, which can be easily verified noting that

1 -1 (1 -x 2 ) α G (α) k (x) dx = (g (α) 0,0 ) -1 if k = 0, and = 0 otherwise, implying 1 -1 (1 -x 2 ) α Kn k=0 ζn,k G (α) k (x)dx = ζn,0 /g (α) 0,0 = σn (0) -σn (2) by (2.8).
In this paper, we further study the properties of the proposed mixture density estimator. In order to formulate the theorem about the asymptotic normality of estimator φn (x), we will assume that aggregated process X t , t ∈ Z admits the following linear representation.

Assumption A Assume that X t , t ∈ Z is a linear sequence (2.11) where the Z t are i.i.d. random variables with zero mean, finite fourth moment and the coefficients ψ j satisfy

X t = ∞ j=0 ψ j Z t-j ,
ψ j ∼ cj d-1 , |ψ j -ψ j+1 | = O(j d-2 ), 0 < d < 1/2 (2.12)
with some constant c = 0.

We also introduce the following condition on the mixture density ϕ(x).

Assumption B Assume that mixture density ϕ has a form

ϕ(x) = (1 -x) 1-2d ψ(x), 0 < d < 1/2, (2.13) where ψ(x) is a nonnegative function with supp(ψ) ⊂ [-1, 1], continuous at x = 1, ψ(1) = 0.
Note that, omitting in (2.10) the factor responsible for the seasonal part, we thus obtain the corresponding 'long memory' spectral density with singularity at zero (but not necessary at ±π) and the corresponding behavior of the coefficients ψ j in linear representation (3.2).

Theorem 2.1 Let X t , t ∈ Z be the aggregated process satisfying Assumption A and corresponding to the mixture density given by Assumption B. Assume that (2.6) holds, and d and α satisfy the following condition

-1/2 < α < 5 2 -4d. (2.14)
Let K n be given in (2.9) with γ satisfying

0 < γ < (2 log(1 + √ 2)) -1 1 -max α + 4d - 3 2 , 0 . (2.15)
Then for every fixed x ∈ (-1, 1), such that ϕ(x) = 0, it holds

φn (x) -E φn (x) Var( φn (x)) d -→ N(0, 1). (2.16)
Proof of the theorem is given in Section 4.

Remark 2.1 Suppose that ϕ(x) satisfies Assumption B. Then assumption (2.6) is equivalent to

1 -1 ψ 2 (x)(1 + x) -α dx < ∞ and α < 3 -4d. The last inequality is implied by (2.14). Example 2.1 Assume two mixture densities ϕ(x; d) = C 1 (d)x d-1 (1 -x) 1-2d (1 + x)1 (0,1] (x), 0 < d < 1/2,
(2.17)

where

C 1 (d) = Γ(3-d) 2Γ(d)Γ(2-2d)
, and

ϕ g (x; κ) = C 2 (κ)|x| κ 1 [-a * ,0] (x), κ > 0, (2.18) where 0 < a * < 1, C 2 (κ) = (κ + 1)(a * ) -κ-1 .
According to Dacunha-Castelle and Oppenheim (2001), the spectral density corresponding to ϕ(x; d) is FARIMA(0,d,0) spectral density

f (λ; d) = 1 2π 2 sin |λ| 2 -2d
.

(2.19) Also, since the support of ϕ g lies inside (-1, 1), the spectral density g(λ; κ) corresponding to ϕ g (x; κ) is analytic function (see Proposition 3.3 in [START_REF] Celov | Time series aggregation, disaggregation and long memory[END_REF]).

Consider the spectral density given by

f (λ) = f (λ; d)g(λ; κ), λ ∈ [-π, π].
(2.20)

It can be shown that the mixture density ϕ(x) associated with f (λ) (2.20) is supported on [-a * , 1], satisfies Assumption B with ψ(x) which is continuous function on [-a * , 1] and at the neighborhood of zero satisfies ψ

(x) = O(|x| d ).
This implies the validity of condition (2.6) needed to obtain the corresponding α-Gegenbauer expansion. For the proof of this example and precise asymptotics of ψ(x) at zero see Appendix A. Finally, the aggregated process X, obtained using such mixture density ϕ(x), satisfies Assumption A by Proposition 3.2, which shows that assumptions A and B are satisfied under general 'aggregated' spectral density f (λ) = f (λ; d)g(λ), where g(λ) is analytic function on [-π, π] and the associated mixture density is supported on [-a * , 0] with some 0 < a * < 1.

Remark 2.2 Note that the 'FARIMA mixture density' (2.17), due to factor x d-1 , does not satisfy (2.6) and a "compensating" density such as ϕ g (x; κ) in (2.18) is needed in order to obtain the needed integrability in the neighborhood of zero. Obviously, for the same aim, other mixture densities instead of ϕ g (x; κ) (2.18) can be employed.

Moving average representation of the aggregated process

In order to obtain the asymptotic normality result in Theorem 2.1, an important assumption is that the aggregated process admits a linear representation with coefficients decaying at an appropriate rate (see [START_REF] Bhansali | Approximations and limit theory for quadratic forms of linear processes[END_REF]). The related issues about the moving average representation of the aggregated process are discussed in this section.

From the aggregating scheme follows that any aggregated process admits an absolutely continuous spectral measure. If, in addition, its spectral density, say, f (λ) satisfies

π -π log f (λ)dλ > -∞, (3.1)
then the function

h(z) = exp 1 4π π -π e iλ + z e iλ -z log f (λ)dλ , |z| < 1,
is an outer function from the Hardy space H 2 , does not vanish for |z| < 1 and f (λ) = |h(e iλ )| 2 . Then, by the Wold decomposition theorem, corresponding process X t is purely nondeterministic and has the MA(∞) representation (see Anderson (1971, Ch. 7.6.3))

X t = ∞ j=0 ψ j Z t-j , (3.2)
where the coefficients ψ j are defined from the expansion of normalized outer function h(z)/h(0), ∞ j=0 ψ 2 j < ∞, ψ 0 = 1, and Z t = X t -X t , t = 0, 1, . . . ( X t is the optimal linear predictor of X t ) is the innovation process, which is zero mean, uncorrelated, with variance

σ 2 = 2π exp 1 2π π -π log f (λ)dλ . (3.3)
By construction, the aggregated process is Gaussian, implying that the innovations Z t are i.i.d. N(0, σ 2 ) random variables.

Next we focus on the class of semiparametric mixture densities satisfying Assumption B. As it was mentioned earlier, this form is natural, in particular it covers the mixture densities ϕ 1 (x; d) and ϕ(x) in Example 2.1.

Proposition 3.1 Let the mixture density ϕ(x) satisfies Assumption B. Assume that either

(i) supp(ψ) = [-1, 1] and ψ(x) ≡ ψ(x)(1 + x) 2 d-1 is continuous at -1 and ψ(-1) = 0 with some 0 < d < 1/2, or (ii) supp(ψ) ⊂ [-a * , 1] with some 0 < a * < 1.
Then the aggregated process admits a moving average representation (3.2), where the Z t are Gaussian i.i.d. random variables with zero mean and variance (3.3).

Proof. (i) We have to verify that (3.1) holds. Rewrite ϕ(x) in the form

ϕ(x) = (1 -x) 1-2d (1 + x) 1-2 d ψ(x).
Proposition 4.1 in [START_REF] Celov | Time series aggregation, disaggregation and long memory[END_REF] implies

f (λ) ∼ C 1 |λ| -2d , |λ| → 0, with C 1 > 0. Hence log f (λ) ∼ log C 1 -C 2 log |1 -e iλ |, |λ| → 0, where C 2 > 0. For any ǫ > 0 choose 0 < λ 0 ≤ π/3, such that - log f (λ) -log C 1 C 2 log |1 -e iλ | -1 ≥ -ǫ, 0 < λ ≤ λ 0 .
Since -log |1 -e iλ | ≥ 0 for 0 ≤ λ ≤ π/3, we obtain

λ 0 0 log f (λ)dλ ≥ λ 0 log C 1 -C 2 (1 -ǫ) λ 0 0 log |1 -e iλ |dλ > -∞ (3.4)
using the well known fact that

π 0 log |1 -e iλ |dλ = 0. Similarly, π π-λ 0 log f (λ)dλ > -∞. (3.5) When λ ∈ [λ 0 , π -λ 0 ], there exist 0 < L 1 < L 2 < ∞ such that L 1 ≤ 1 2π|1 -xe iλ | 2 ≤ L 2 uniformly in x ∈ (-1, 1). Thus, by (2.5), L 1 ≤ f (λ) ≤ L 2 for any λ ∈ [λ 0 , π -λ 0 ], and therefore π-λ 0 λ 0 log f (λ)dλ > -∞. (3.6) (3.4)-(3.6) imply inequality (3.1).
The proof in case (ii) is analogous to (i) and, thus, is omitted. 2

Lemma 3.1 If the spectral density g(λ) of the aggregated process X t , t ∈ Z is analytic function on [-π, π], then X t admits representation

X t = ∞ j=0 g j Z t-j ,
where the Z t are i.i.d. Gaussian random variables with zero mean and variance

σ 2 g = 2π exp 1 2π π -π log g(λ)dλ (3.7) and the g j satisfy | ∞ j=0 g j | < ∞, g 0 = 1.
Proof. From Proposition 3.3 in [START_REF] Celov | Time series aggregation, disaggregation and long memory[END_REF] it follows that there exists 0 < a * < 1 such that

g(λ) = σ 2 ε 2π a * -a * ϕ g (x) |1 -xe iλ | 2 dx. (3.8) For all x ∈ [-a * , a * ] and λ ∈ [0, π] we have 1 |1 -xe iλ | 2 ≥ C 3 > 0,
where C 3 = C 3 (a * ). This and (3.8) imply

π 0 log g(λ)dλ > -∞. Finally, | ∞ j=0 g j | < ∞ follows from representation g(λ) = σ 2 g 2π ∞ j=0 g j e ijλ 2
and the assumption of analyticity of g. 2 Proposition 3.2 Let X t , t ∈ Z be an aggregated process with spectral density

f (λ) = f (λ; d)g(λ), (3.9)
where f (λ; d) is FARIMA spectral density (2.19) and g(λ) is analytic spectral density. Then:

(i) if mixture density ϕ g (x) associated with g(λ) satisfies supp(ϕ g ) ⊂ [-a * , 0] with some 0 < a * < 1, then ϕ(x), associated with f (λ), satisfies Assumption B.

(ii) X t admits a linear representation (3.2), where the Z t are Gaussian i.i.d. random variables with zero mean and variance

σ 2 = 2π exp 1 2π π -π log f (λ)dλ = exp 1 2π π -π log g(λ)dλ = σ 2 g 2π
and the coefficients ψ j satisfy

ψ j ∼ ∞ k=0 g k Γ(d) j d-1 , |ψ j -ψ j+1 | = O(j d-2 ), (3.10)
where ψ 0 = 1. (Here, the g k are given in Lemma 3.1.)

Proof. (i) By Corollary 3.1 in [START_REF] Celov | Time series aggregation, disaggregation and long memory[END_REF], the mixture density associated with the "product" spectral density (3.9) exists and has a form where ϕ(x; d) is given in (2.17) and is associated with the spectral density f (λ; d), and ϕ g (x) is associated with the spectral density g(λ). Clearly, this implies that Assumption B is satisfied.

ϕ(x) = C -1 * ϕ(x; d) 0 -a * ϕ g (y)dy (1 -xy)(1 -y/x) +ϕ g (x) 1 0 ϕ(y; d)dy (1 -xy)(1 -y/x) , (
(ii) We have

f (λ; d) = 1 2π ∞ j=0 h j e ijλ 2 with h j = Γ(j + d) Γ(j + 1)Γ(d)
and, recall,

g(λ) = σ 2 g 2π ∞ j=0 g j e ijλ 2 , ∞ j=0 g 2 j < ∞ since, by Lemma 3.1, π -π log g(λ)dλ > -∞. On the other hand, π -π log f (λ)dλ > -∞ implies f (λ) = 1 2π ∞ j=0 ψj e ijλ 2 , ∞ j=0 ψ2 j < ∞
and, by uniqueness of the representation,

ψk = σ g √ 2π k j=0 h k-j g j .
It easy to see that,

k j=0 h k-j g j ∼ h k ∞ j=0 g j ∼ C 4 k d-1 , (3.13)
where

C 4 = Γ -1 (d) ∞ j=0 g j . Indeed, taking into account that h k ∼ Γ -1 (d)k d-1 , we can write k j=0 h k-j g j = Γ -1 (d)k d-1 ∞ j=0 a k,j g j ,
where a k,j = h k-j Γ(d)k 1-d 1 {j≤k} → 1 as k → ∞ for each j. On the other hand, we have |a k,j | ≤ C(1 + j) 1-d uniformly in k and, since the g j decay exponentially fast, the sum ∞ j=0 (1+ j) 1-d |g j | converges and the dominated convergence theorem applies to obtain (3.13).

Hence, we can write

f (λ) = σ 2 g (2π) 2 ∞ j=0 ψ j e ijλ 2 , ψ 0 = 1,
where

ψ j = ψj √ 2π/σ g ∼ C 4 j d-1
. Thus, representation (3.2) and the first relation in (3.10) follows.

Finally, in order to check the second relation in (3.10), it suffices to note that

ψ j -ψ j+1 = j i=0 (h j-i -h j+1-i )g i -g j+1 ,
where h j -h j+1 ∼ C 5 j d-2 and the g j decay exponentially fast. 2

Proof of main result

In order to prove Theorem 2.1, we use the result of [START_REF] Bhansali | Approximations and limit theory for quadratic forms of linear processes[END_REF], who considered the following quadratic form

Q n,X = n t,s=1 d n (t -s)X t X s ,
where the X t are linear sequences satisfying Assumption A and the function d n (k) satisfies the following assumption.

Assumption C Suppose that

d n (k) = π -π η n (λ)e ikλ dλ
with some even real function η n (λ), such that, for some -1 < β < 1 and a sequence of constants m n ≥ 0, it holds

|η n (λ)| ≤ m n |λ| -β , λ ∈ [-π, π]. (4.1)
Denote by E n a matrix (e n (t -s)) t,s=1,...,n , where

e n (t -s) = π -π η n (λ)f (λ)e iλ(t-s) dλ (4.2)
and let E n 2 = n t,s=1 e 2 n (t -s).

Theorem 4.1 [START_REF] Bhansali | Approximations and limit theory for quadratic forms of linear processes[END_REF]] Suppose that assumptions A and C are satisfied. If 2d + β < 1/2 and

r n = o( E n ), (4.3)
where

r n =    m n n max(0,2d+β) if 2d + β = 0, m n log n if 2d + β = 0, (4.4) then, as n → ∞, it holds Var(Q n,X ) ≍ E n 2 and Q n,X -EQ n,X Var(Q n,X ) d -→ N(0, 1). (Here for a n , b n ≥ 0, a n ≍ b n means that C 6 b n ≤ a n ≤ C 7 b n for some C 6 , C 7 > 0.)
Proof of Theorem 2.1. First of all, note that σ2

n,ε P -→ σ 2 ε ,
which easily follows using Theorem 3 in [START_REF] Hosking | Asymptotic distributions of the sample mean, autocovariances and autocorrelations of long-memory time series[END_REF]. Hence, to obtain convergence (2.16), we can replace the factor σ2 n,ε by σ 2 ε in the definition of φn (x). Without loss of generality assume that σ 2 ε = 1. Rewrite the estimate φn (x) in a form

φn (x) = (1 -x 2 ) α Kn k=0 k j=0 g (α) k,j (σ n (j) -σn (j + 2))G (α) k (x) = (1 -x 2 ) α Kn k=0 G (α) k (x) k j=0 g (α) k,j π -π (e iλj -e iλ(j+2) )I n (λ)dλ = π -π η n (λ; x)I n (λ)dλ, (4.5)
where

η n (λ; x) := (1 -x 2 ) α Kn k=0 G (α) k (x) k j=0 g (α)
k,j (e iλj -e iλ(j+2) ) (4.6)

and

I n (λ) = (2πn) -1 | n j=1 X j e ijλ | 2 , λ ∈ [-π, π]
is the periodogram. Now the proof follows from Assumption A and the results obtained in Lemma 4.1 and Lemma 4.2 below, which imply that, under appropriate choice of m n and β, all the assumptions in Theorem 4.1 are satisfied. In particular, by Lemma 4.1, the following bound for the kernel η n (λ; x) holds

|η n (λ; x)| ≤ m n |λ| -β , (4.7)
where

m n = C 8 n γ log(1+ √ 2) , β = α 2 - 3 4 , (4.8)
C 8 is a positive constant, depending on x and α. Clearly, (2.14) implies that -1 < β ≤ 1 2 -2d < 1 2 and 2d + β < 1 2 . Consider the cases 2d+β ≤ 0 or 0 < 2d+β < 1/2. In the case 2d+β ≤ 0, from (4.4), (4.8) we obtain

r n = C 8    n γ log(1+ √ 2) if 2d + α 2 -3 4 < 0, n γ log(1+ √ 2) log n if 2d + α 2 -3 4 = 0.
Hence, by Lemma 4.2, r n E n -1 → 0 because γ log(1 + √ 2) < 1/2. Assume now 2d + β > 0. Then

r n = C 8 n γ log(1+ √ 2)+2d+ α 2 -3 4 and r n E n -1 → 0 by (2.15). 2
The following lemma shows that the kernel η n (λ; x) given in (4.6) satisfies inequality (4.7) with m n and β given in (4.8).

Lemma 4.1 For quantity η n (λ; x) given in (4.6) and for every fixed x ∈ (-1, 1), 0 < |λ| < π it holds

|η n (λ; x)| ≤ C 9 n γ log(1+ √ 2) |λ| (3-2α)/4    (1 -x 2 ) α/2-1/4 if α > -1/2, (1 -x 2 ) α if -1 < α < -1/2,
where C 9 depends on α, and γ is given in (2.9).

Lemma 4.2 Assume that a mixture density ϕ(x) satisfies condition (2.6) and let K n → ∞. Then for every x ∈ (-1, 1), such that ϕ(x) = 0 it holds

E n 2 ≥ C 10 n(1 + o(1)), (4.9)
where C 10 > 0 is positive constant depending on α and x.

Proof of these two lemmas are given in Appendix B.

A simulation study

In order to gain further insight into the asymptotic normality property of the mixture density estimator (2.7), in this section we conduct a Monte-Carlo simulation study. Several examples are considered, which correspond to the mixture densities having different shapes (here we do not pose a question which rigorous aggregating schemes lead to the latter).

The following two families of mixture densities

ϕ(x) = wϕ 1 (x) + (1 -w)ϕ 2 (x), 0 < w < 1,
are considered:

• Beta-type mixture densities defined by

ϕ 1 (x) ∝ x p 1 -1 (1 -x) q 1 -1 1 [0,1] (x), p 1 > 0, q 1 > 0, ϕ 2 (x) ∝ |x| p 2 -1 (a * + x) q 2 -1 1 [-a * ,0] (x), p 2 > 0, q 2 > 0, 0 < a * < 1;
• mixed (Beta and Uniform)-type mixture densities defined by

ϕ 1 (x) ∝ x p 3 -1 (1 -x) q 3 -1 1 [0,1] (x), p 3 > 0, q 3 > 0, ϕ 2 (x) = a -1 * 1 [-a * ,0] (x), 0 < a * < 1.
In order to construct the mixture density estimator, in the first step, the parameters K n and α must be chosen. Preliminary Monte-Carlo simulations showed that the estimator φn (x) has the minimal mean integrated square error (MISE) when the parameter α is chosen to be equal 1 -2d. The justification of this interesting conjecture remains an open problem. This rule also ensures that (2.14) is satisfied. The number of Gegenbauer polynomials K n is chosen according to (2.9). Note that, by construction, the estimator φn (x) is not necessary positive, though it integrates to one.

In Figure 1, we present three graphs and corresponding box plots for the mixture densities of the form above. Cases 1 and 2 correspond to the Betatype mixture densities, Case 3 corresponds to the mixed (Beta and Uniform)type mixture density. The parameter values are presented in Table 1. The box plots are obtained by a Monte-Carlo procedure based on M = 500 independent replications with sample size n = 1500 and bandwidth K n = 3 (we aggregate N = 5000 i.i.d. AR(1) processes). Individual innovations ε (j) t are i.i.d. N(0, 1). Note that the mixture density in Case 2 corresponds to Example 2.1 with the parameters d = 0.2, κ = 0.1 (in the sense of behavior at zero).

w a * (p 1 , q 1 ) (p 2 , q 2 ) (p 3 , q 3 ) d α Case 1 0.8 0.95 (3.0, 1.5) (2.0, 1.0) - 0.25 0.5 Case 2 0.8 0.80 (1.2, 1.6) (1.3, 2.5) - 0.20 0.6 Case 3 0.8 0.90 - - (2.0, 1.2) 0.40 0.2
Table 1: Parameter values in cases 1-3.

Box plots in Figure 1 show that φn approximates the mixture density well when n is sufficiently large. However, when the sample size is relatively small it is difficult to estimate the mixture density of the shape as in cases 2-3. This can be explained by the construction of the estimator which assumes rather smooth form of the mixture density around zero. On the other hand, it is clear that the AR(1) parameter values which are close to zero does not affect the long memory property. For our purposes, an important fact is that the estimator correctly approximates the density at the neighborhood of x = 1. This enables us to estimate the unknown (in real applications) parameter d using a log-log regression on periodogram at the neighborhood of this point (for example Geweke and Porter-Hudak or Whittle-type estimators).

Figure 2 supplements the earlier findings and shows that the distribution of estimator is approximately normal.1 QQ-plots and histograms are given for fixed values x = -0.5 and x = 0.96 correspondingly. We use the same number of replications M = 500 and sample size n = 1500. The last Monte-Carlo experiment aims to show that the decay rate of Var( φn (x)) is n -γ with γ = 1. This ensures that the variance is decreasing fast enough. To do this, we calculate the log-log regression of variance on the length of time series n ∈ {500, 600, . . . , 1400, 1500, 2000, . . . , 5000}. Figure 3 demonstrates the corresponding parameter estimates at different points and shows that γ ≈ 1. where

-3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 Case 1, x = -0.5 Case 1, x = -0.5 -3 -2 -1 0 1 2 3 0.0 0.1 0.2 0.3 0.4 0.5 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 Case 1, x = 0.96 Case 1, x = 0.96 -3 -2 -1 0 1 2 3 0.0 0.1 0.2 0.3 0.4 0.5 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 Case 2, x = -0.5 Case 2, x = -0.5 -3 -2 -1 0 1 2 3 0.0 0.1 0.2 0.3 0.4 0.5 -3 -2 -1 0 1 2 3 -3 -2 -1 0
C = C 1 (d)C 2 (κ)C -1 * is positive constant, ψ 1 (x) := x d-1 (1 + x)1 (0,1] (x) 0 -a * |y| κ (1 -xy)(1 -y/x)
dy, (6.2)

ψ 2 (x) := |x| κ (1 -x) 2d-1 1 [-a * ,0] (x) 1 0 y d-1 (1 -y) 1-2d (1 + y) (1 -xy)(1 -y/x) dy. (6.3)
Denote by F (a, b; c; x) a hypergeometric function

F (a, b; c; x) = Γ(c) Γ(b)Γ(c -b) 1 0 t b-1 (1 -t) c-b-1 (1 -tx) -a dt,
with c > b > 0 if x < 1 and, in addition, c -a -b > 0 if x = 1. Then the corresponding integrals in ψ 1 (x) and ψ 2 (x) can be rewritten as

0 -a * |y| κ (1 -xy)(1 -y/x) dy, = a κ+1 * κ + 1 x(F (1, κ + 1; κ + 2; -a * x) -F (1, κ + 1; κ + 2; -a * /x)) 1 -x 2 ∼ a κ+1 * κ + 1 x, as x → 0+,
and

1 0 y d-1 (1 -y) 1-2d (1 + y) (1 -xy)(1 -y/x) dy = Γ(d)Γ(2 -2d) Γ(2 -2d) F (1, d; 2 -d; 1/x) -xF (1, d; 2 -d; x) 1 -x ∼ Γ(d)Γ(1 -d)|x| d , as x → 0-,
where the last asymptotics follow from the well known properties of the hypergeometric functions (see [START_REF] Abramovitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF]). Thus, from (6.2)-( 6.3) we obtain that

ψ 1 (x) ∼ a κ+1 * κ + 1 x d , as x → 0+, ( 6
-x 2 ) -α η n (λ; x) = Kn k=0 G (α) k (x) k j=0 g (α) k,j (e iλj -e iλ(j+2) ) = (1 -e 2iλ ) Kn k=0 G (α) k (x) k j=0 g (α) k,j e iλj = (1 -e 2iλ ) Kn k=0 G (α) k (x)G (α) k (e iλ ). 1
This and Lemma 7.1 below implies

(1 -x 2 ) -α |η n (λ; x)| ≤ C 11 |λ| -(2α-3)/4 Kn k=0 |G (α) k (x)|(1 + √ 2) k . (7.1)
Now, using the fact that for all -1 < x < 1

|G (α) k (x)| ≤    C 12 (1 -x 2 ) -α 2 -1 4 if α > -1/2 C 12 if α < -1/2, α = -3/2, -5/2, . . .
(see inequality (7.33.6) in [START_REF] Szegö | Orthogonal Polynomials[END_REF] and (3.9) in [START_REF] Leipus | Orthogonal series density estimation in a disaggregation scheme[END_REF]) and (2.9), we get from (7.1)

(1 -x 2 ) -α |η n (λ; x)| ≤ C 13 |λ| -(2α-3)/4 (1 + √ 2) Kn = C 13 |λ| -(2α-3)/4 e Kn log(1+ √ 2) ≤ C 9 |λ| -(2α-3)/4 n γ log(1+ √ 2) . 2 Lemma 7.1 For all k ≥ 0, α > -1, (α = -1/2) and 0 < |λ| < π it holds |(1 -e 2iλ )G (α) k (e iλ )| ≤ C 11 (1 + √ 2) k |λ| -(2α-3)/4 ,
where constant C 11 depends on α.

Proof. Theorem 8.21.10 of [START_REF] Szegö | Orthogonal Polynomials[END_REF] implies that for the usual (nonnormalized) Gegenbauer polynomials with α > -1, α = -1/2 it holds

C (α+1/2) k (e iλ ) = Γ(k + α + 1 2 ) Γ(k + 1)Γ(α + 1 2 ) z k (1-z -2 ) -α-1/2 +O(k α-3/2 |z| k ), (7.2)
where the complex numbers w = e iλ and z are connected by the elementary conformal mapping

w = 1 2 (z + z -1 ), z = w + (w 2 -1) 1/2 , (7.3)
and z satisfies |z| > 1 (thus, λ = 0, ±π).

Recall that the normalized Gegenbauer polynomials G

(α) k (z) are linked to C (α+1/2) k (z) by equality G (α) k (z) = γ -1/2 k C (α+1/2) k (z), where γ k = π 2 2α Γ(k + 2α + 1) (k + α + 1 2 )Γ 2 (α + 1 2 )Γ(k + 1)
.

Therefore, in terms of the normalized Gegenbauer polynomials, (7.2) reads as follows

G (α) k (e iλ ) = sgn(α + 1/2)2 α π 1/2 b k z k (1 -z -2 ) -α-1/2 + O(k -1 |z| k ), (7.4)
where

b k = (k + α + 1/2) 1/2 Γ(k + α + 1/2) Γ 1/2 (k + 1)Γ 1/2 (k + 2α + 1) → 1 as k → ∞. From (7.3) we obtain for w = e iλ w 2 -1 = 1 4 z 2 (1 -z -2 ) 2 ,
which together with (7.4) yields

(1 -e 2iλ )G (α)
k (e iλ ) = -

sgn(α + 1/2)2 α 4π 1/2 b k z k+2 (1 -z -2 ) -α+3/2 + O(k -1 |z| k ).
Since |z| > 1 and z 2 -1 = 2(e 2iλ -1) + 2e 3iλ/2 (e iλ -e -iλ ) 1/2 , we have

|1 -z -2 | ≤ |z 2 -1| ≤ 2|e 2iλ -1| + 2|e iλ -e -iλ | 1/2 = 4| sin λ| + 2 √ 2| sin λ| 1/2 ≤ (4 + 2 √ 2)|λ| 1/2 . (7.5)
So that, by (7.4)-(7.5),

|(1 -e 2iλ )G (α) k (e iλ )| ≤ C 14 b k |z| k |λ| -(2α-3)/4 , (7.6)
where C 14 = C 14 (α). Finally, the straightforward verification shows that sup λ∈[-π,π]

|e iλ + (e 2iλ -1) 1/2 | = 1 + √ 2.
This completes the proof of lemma.

2

Proof of Lemma 4.2. Using (4.2), (4.6) rewrite the coefficients of

E n e n (t-s) = (1-x 2 ) α Kn k=0 G (α) k (x) k j=0 g (α) k,j π -π
f (λ)(e iλ(t-s+j) -e iλ(t-s+j+2) )dλ.

Using the expression of the covariance function of an aggregated process, we have for t -s + j ≥ 0 (7.9)

Since the last term A 3,n is nonnegative by construction, this will prove (4.9). At points x where ϕ(x) = 0 we have

A 1,n = ϕ 2 (x) (1 -x 2 ) 2α |m|<n (n -|m|)x 2|m|
∼ n ϕ 2 (x)(1 + x 2 ) (1 -x 2 ) 2α+1 , as n → ∞, which gives (7.8). Consider term A 2,n . By (7.7), using the similar argument as in the case of term B 1,n . This completes the proof of (7.9) and of the lemma.

A 2,n = |m|<n (n -|m|)ψ |m| (x) ∞ k=Kn+1 G (α) k (x)ψ |m|,k = ϕ(x) (1 -x 2 ) α ∞ k=Kn+1 G (α) j (x)
2

  .7) where the ζn,k are estimates of the coefficients ζ k in the α-Gegenbauer expansion of the function ζ(x) = ∞ k=0 ζ k G

Figure 1 :

 1 Figure 1: True mixture densities (solid line) and the box plots of the estimates. Number of replications M = 500, sample size n = 1500.

Figure 2 :

 2 Figure 2: QQ plots and histograms of the estimates at points x = -0.5 and x = 0.96. Number of replications M = 500, sample size n = 1500.

Figure 3

 3 Figure 3: log-log scale regression of the variance of φn (x) as a function of n. The variance is estimated using M = 500 independent replications.

  .4)ψ 2 (x) ∼ Γ(d)Γ(1 -d)|x| κ+d , as x → 0 -.

  (e iλ(t-s+j) -e iλ(t-s+j+2) )dλ = σ(t -s + j) -σ(t -s + j + 2) s+j ϕ(y)dy.Thus, assuming σ2 ε = 1, for t -s ≥ 0 we have e n (t -s) = (1 -x 2) dy (m is a nonnegative integer), appearing in the last expression is nothing else but the kth coefficient, ψ m,k , in the α-Gegenbauer expansion of the functionψ m (x) = x m ϕ(x) (1 -x 2 ) α ,(7.7)which obviously satisfies ψ m ∈ L 2 (w (α) ). Therefore,e n (t -s) = (1 -x 2 ) x 2 ) α ψ |t-s| (x) -ψ |m|,k , |m| < n, we have (1 -x 2 ) -2α E n 2 = |m|<n (n -|m|) ψ |m| (x) -|m|)ψ |m| (x)R n (m) + |m|<n (n -|m|)R 2 n (m) =: A 1,n -2A 2,n + A 3,n .Now, we prove that, as n → ∞,A 1,n ∼ C 15 n,(7.8)where C 15 = C 15 (x) > 0 is some positive constant, andA 2,n = o(n).

  x 2 ) α (B 1,n -B 2,n -B 3,n )y 2 ) α 1 + xy 1 -xy satisfies φx ∈ L 2 (w (α) ) and K n → ∞, the sum ∞ k=Kn+1 in B 1,n vanishes (as the tail of the convergent series). So that, B 1,n = o(n) and, similarly, B 3,n = o(n). xy) 2 dy = o(1)

The Shapiro-Wilk test confirms that in most cases normality hypothesis is consistent with the data.

Appendix A. Proof of Example 2.1 By Corollary 3.1 in Celov et al. (2007), the mixture density ϕ(x), x ∈ [-a * , 1] associated with f (λ) (2.20) is given by equality (3.11), where ϕ g (x) ≡

The research was supported by bilateral France-Lithuanian scientific project Gilibert and Lithuanian State Science and Studies foundation (V-07058).


Guillaume Aupy

Olivier Beaumont

Lionel Eyraud-Dubois

Sizing and Partitioning Strategies for Burst-Buffers to Reduce IO Contention

Keywords: scheduling, IO, PFS, Burst-Buffers ordonnancement, IO, PFS, Burst-Buffers

Burst-Buffers are high throughput and small size storage which are being used as an intermediate storage between the Parallel File System (Parallel File System) and the computational nodes of modern HPC systems. They can allow to hinder to contention to the Parallel File System, a shared resource whose read and write performance increase slower than processing power in HPC systems. A second usage is to accelerate data transfers and to hide the latency to the Parallel File System. In this paper, we concentrate on the first usage. We propose a model for Burst-Buffers and application transfers. We consider the problem of dimensioning and sharing the Burst-Buffers between several applications. This dimensioning can be done either dynamically or statically. The dynamic allocation considers that any application can use any available portion of the Burst-Buffers. The static allocation considers that when a new application enters the system, it is assigned some portion of the Burst-Buffers which cannot be used by the other applications until that application leaves the system and its data is purged from it. We show that the general sharing problem to guarantee fair performance for all applications is an NP-Complete problem. We give a polynomial time algorithms for the special case of finding the optimal buffer size such that no application is slowed down due to Parallel File System contention, both in the static and dynamic cases. Finally, we provide evaluations of our algorithms in realistic settings. We use those to discuss how to minimize the overhead of the static allocation of buffers compared to the dynamic allocation.

Introduction

The I/O bottleneck is becoming a major issue in current HPC systems. This is especially striking when considering their recent evolution. For instance, when Los Alamos National Laboratory moved from Cielo to Trinity, the peak performance moved from 1.4 Petaflops to 40 Petaflops (×28) while the I/O bandwidth moved to 160 GB/s to 1.45TB/s (only ×9) [START_REF]The trinity project[END_REF]. The same kind of results can be observed at Argonne National Laboratory when moving from Intrepid (0.6 PF, 88 GB/s) and to Mira (10PF, 240 GB/s) and to Aurora (expected 180PF and 1TB/s) [1]. The main storage technology is still based on hard drives, that have shown a better capacity to scale up in terms of storage capacity than speed. On the other hand, the major issue is still on handling I/O peaks, and the overall average usage is still relatively low. For instance, taking the former supercomputer Intrepid at Argonne National Laboratory as an example, the aggregate I/O throughput is lower than one-third of the peak I/O bandwidth for 99% of the time [START_REF] Carns | Storage access characteristics of computational science applications[END_REF][START_REF] Carns | Understanding and improving computational science storage access through continuous characterization[END_REF]. Nevertheless, applications running regularly perform I/Os, to read their initial data, to write intermediate and final results on reliable storage and to enforce reliability using checkpoint restart techniques. In typical HPC systems, very few applications (usually one) are enough to saturate the bandwidth to the Parallel File System, so that delays are experienced if the transfers are not coordinated. Mechanisms such as Clarisse [START_REF] Isaila | Clarisse: A middleware for data-staging coordination and control on large-scale hpc platforms[END_REF] have thus been introduced to reorganize transfers to the Parallel File System at system level.

On the other hand, the usage of HPC systems makes I/O more and more crucial and complex. First, in the framework on the convergence between HPC and BigData [START_REF] Reed | Exascale computing and big data[END_REF], HPC systems are now also used to run BigData applications. One main characteristic of BigData workload is that they are dominated by read operations. Second, the MTBF (Mean Time Between Failures) of HPC systems is decreasing [START_REF] Cappello | Toward exascale resilience: 2014 update[END_REF][START_REF] Bosilca | Unified model for assessing checkpointing protocols at extreme-scale[END_REF] and Checkpoint/Restart (C/R) strategies are necessary to enforce reliable computations in a failure prone system. C/R strategies mostly induce write operations. Third, HPC applications themselves consume a lot of I/O bandwidth (see Section 3.2) and they alternate read, compute and write phase to the Parallel File System, that cannot be overlapped. We consider in the present paper both read and write accesses to the Parallel File System.

When running several such applications, even if the overall bandwidth is enough to cope in the long term with the required data transfers, the bursty nature of both read and write operations and the lack of synchronization between applications induces I/O peaks, that in turn degrade the aggregated bandwidth, as noted in [START_REF] Tang | Toward managing hpc burst buffers effectively: Draining strategy to regulate bursty i/o behavior[END_REF]. In this context, in order to cope with the limited I/O bandwidth of HPC system, Burst-Buffers have emerged as promising solution [START_REF] Ddn Storage | [END_REF][START_REF] Henseler | Architecture and design of cray datawarp[END_REF][START_REF] Daley | Performance characterization of scientific workflows for the optimal use of burst buffers[END_REF].

Burst-Buffers are first used to accelerate transfers and to hide the latency and the limited bandwidth to the Parallel File System, by acting as a cache between the computational nodes and the Parallel File System. In practice, this use is arguable since on the technological side, the use of NVRAM or SSDs makes it possible to achieve much higher bandwidth than hard disks, but at the price of a limited lifetime [START_REF] Han | Accelerating a burst buffer via user-level i/o isolation[END_REF]. Another potential use of Burst-Buffers is to enable to delay and better schedule the I/O operations to the Parallel File System, by acting as an intermediate storage used to delay write operations and to prefetch read operations, in order to avoid access conflicts and to hide contentions to the user by dealing smoothly with I/O peaks.

There is still no clear consensus on Burst-Buffers architecture (see Section 2.1). In this paper, we consider the simplest model where the Burst-Buffers are not distributed and acts as a potential intermediate centralized layer, with a higher I/O bandwidth but a smaller capacity than the hard disk storage system. It can however be partitioned between the different applications, ensuring that each application has a dedicated allocation.

In the present paper, we therefore consider a set of applications running independently on different computational nodes belonging to the same machine, where the allocation of nodes has been done a priori using a batch scheduler such as SLURM [START_REF] Yoo | Slurm: Simple linux utility for resource management[END_REF]. In order to deal with BigData, Checkpoint/Restart and HPC applications, we consider that the pattern of I/O and processing phases is known in advance, typically through monitoring and historical data [START_REF] Aupy | Periodic i/o scheduling for super-computers[END_REF]. Our goal is both to dimension the Burst-Buffers and to partition it between the different applications in order to limit the application slowdown experienced due to the limited I/O bandwidth to the Parallel File System.

The main contributions of this work are the following:

• A precise and complete model of the platform and applications. We propose several uses of Burst-Buffers (static allocation and dynamic allocation).

• The intractability of the problem in the general case: for a given buffer size, find an allocation that ensures that all applications are treated fairly.

• An optimal algorithm for the problem of minimizing the Burst-Buffers size to ensure that no application is delayed because of inter-application competition to the Parallel File System bandwidth (optimal stretch).

• Finally we provide several evaluations with bounded buffer size to compare the performance of the static and dynamic buffer allocations. We discuss these strategies.

Compared to our previous work [START_REF] Aupy | What size should your buffers to disks be?[END_REF], we consider other types of usage of the Burst-Buffers. On the system side, we consider the possibility of partitioning the Burst-Buffer either statically (give a fixed share of Burst-Buffer to each application) or dynamically (the share of an application can vary with time). We consider a precise application model, where the different phases (in terms of I/O and computation volumes) are known in advance, whereas [START_REF] Aupy | What size should your buffers to disks be?[END_REF] relied on random transfer patterns, where only the probability of being involved in a transfer at any point in time was known. This allows us to provide a more precise estimation of the required Burst-Buffer size, by considering actual events when the buffers are emptied and when the I/O operations between applications interfere.

The rest of this paper is organized as follows. In Section 2, we present the related work on Burst-Buffers architecture, bandwidth allocation and HPC applications models. Then, we detail in Section 3 the application model described above and the architecture of the Burst-Buffers. In Section 4, we prove that the problem of partitioning the Burst-Buffers between several applications is NP-Complete in the strong sense. Despite this result, we prove in Section 6 that a special case can be solved in polynomial time, where the goal is to find both the optimal size of the Burst-Buffers and its optimal partitioning between applications so that the completion time of all the applications is exactly the same as what it would be if each application was running alone with a full access to the Parallel File System. We then compare static and dynamic settings through simulation results in Section 7, and discuss the performance of a classic greedy fair sharing strategy. At last, we discuss the current limitations of our model and possible extensions of the present work and provide concluding remarks in Section 8.

Note that all our results are validated through thorough evaluation. The source code and scripts for those evaluations are available at https://gitlab.inria.fr/ordo-bdx/parted-buffers.

Related Work

Burst-Buffer Architectures and models

There are many implementations of Burst-Buffers. The two most studied characteristics are the location of the buffers and whether the buffers are shared between multiple applications.

Inria

A typical architecture consists in locating the Burst-Buffers between the compute nodes and the Parallel File System (PFS). This is the case of DDN IME [START_REF] Ddn Storage | [END_REF][START_REF] Schenck | Evaluation and performance modeling of a burst buffer solution[END_REF] and Cray DataWarp [START_REF] Henseler | Architecture and design of cray datawarp[END_REF][START_REF] Mubarak | Quantifying i/o and communication traffic interference on dragonfly networks equipped with burst buffers[END_REF][START_REF] Daley | Performance characterization of scientific workflows for the optimal use of burst buffers[END_REF]. In this pseudo-centralized architecture, the Burst-Buffers are often colocated with the I/O nodes. Several management strategies have been proposed. Mubarak et al. [START_REF] Mubarak | Quantifying i/o and communication traffic interference on dragonfly networks equipped with burst buffers[END_REF] study the case where the buffers are shared between the different applications on the platform and used to accelerate transfers and to prevent I/O congestion. On the contrary, in Schenck et al. [START_REF] Schenck | Evaluation and performance modeling of a burst buffer solution[END_REF] and Daley et al. [START_REF] Daley | Performance characterization of scientific workflows for the optimal use of burst buffers[END_REF], applications decide the size of the buffer that should be dedicated to them.

Another solution is a distributed version of Burst-Buffers where the buffers are allocated closer to the compute nodes [START_REF] Kimpe | Integrated in-system storage architecture for high performance computing[END_REF][START_REF] Bautista-Gomez | Fti: high performance fault tolerance interface for hybrid systems[END_REF]. A solution consists in allocating the distributed buffers to the applications using the compute nodes close to buffers [START_REF] Cao | To share or not to share: comparing burst buffer architectures[END_REF]. However, other strategies focus on how to share them between the different applications [START_REF] Kimpe | Integrated in-system storage architecture for high performance computing[END_REF][START_REF] Bautista-Gomez | Fti: high performance fault tolerance interface for hybrid systems[END_REF]. In [START_REF] Mubarak | Quantifying i/o and communication traffic interference on dragonfly networks equipped with burst buffers[END_REF], the interaction between the placement of Burst-Buffers and high-radix interconnect topologies is studied. In the context of fault-tolerance, using a buffer on a different node can allow the implementation of hierarchical checkpointing strategies that provide more resilience than in-node buffer strategies [START_REF] Bautista-Gomez | Fti: high performance fault tolerance interface for hybrid systems[END_REF]. Furthermore, in the case where, because of their costs, the number of buffers in the machine has to be limited, one must choose on which node they should be deployed and between which subset of applications they should be shared.

In [START_REF] Tang | Toward managing hpc burst buffers effectively: Draining strategy to regulate bursty i/o behavior[END_REF], the authors consider the use of Burst-Buffers to serve the I/O bursts of HPC applications. They prove that a basic reactive draining strategy that empties the Burst-Buffer as soon as possible can lead to a severe degradation of the aggregate I/O throughput. On the other hand, they advocate for a proactive draining strategy, where data is divided into draining segments which are dispersed evenly over the I/O interval, and the burst buffer draining throughput is controlled through adjusting the number of I/O requests issued each time. Contrarily to our present work, [START_REF] Tang | Toward managing hpc burst buffers effectively: Draining strategy to regulate bursty i/o behavior[END_REF] does not consider the dimensioning nor the optimal partitioning of the Burst-Buffer, but concentrates on dynamic strategies to actually perform transfers.

Algorithms to deal with Burst-Buffers

When it comes to using Burst-Buffers, several solutions have been proposed. We present and discuss the most common ones.

A natural idea is to use Burst-Buffers as a cache to improve the I/O-performance of applications [START_REF] Schenck | Evaluation and performance modeling of a burst buffer solution[END_REF]. For instance, DDN [START_REF] Ddn Storage | [END_REF] announces bandwidth performance 10-fold that of PFS using their Burst-Buffers. The idea is to move the I/O to the Burst-Buffer as a temporary stage between compute nodes and the Parallel File System (whether the data is incoming or outgoing). Thanks to the higher bandwidth of the Burst-Buffers, this has the advantage of improving the I/O transfer time while pipelining the (slowest) phase of sending/receiving data from the PFS with the compute phase of the application. However, as was noted by Han et al. [START_REF] Han | Accelerating a burst buffer via user-level i/o isolation[END_REF], this idea is not viable, (i) Burst-Buffers are based on technologies that are extremely expensive with respect to hard drives, (ii) they are currently based on SSD technology, that is known to have a limited rewrite lifespan [START_REF] Han | Accelerating a burst buffer via user-level i/o isolation[END_REF]. Thus, the large number of I/O operations in HPC applications would decrease their lifespan too fast. In the solution we propose, not all data transfers go through the Burst-Buffer but only those necessary to avoid I/O congestion.

The second natural idea proposed in the literature is indeed to use Burst-Buffers to prevent I/O congestion [START_REF] Liu | On the role of burst buffers in leadership-class storage systems[END_REF][START_REF] Herbein | Scalable i/o-aware job scheduling for burst buffer enabled hpc clusters[END_REF] while maintaining their lifespan. To achieve this goal, the applications use the direct link to the PFS (see Figure 1) when its bandwidth B is not exceeded. When the bandwidth is exceeded by the set of transfers, then the higher bandwidth of the Burst-Buffer is used to complement the bandwidth of the PFS. This is one of the solutions advocated by DDN in [START_REF] Ddn Storage | [END_REF]. The intuition behind this strategy is that the average use of PFS bandwidth is usually small enough, but that Burst-Buffers are crucial to deal with applications' (simultaneous) bursts. This corresponds to the model depicted in Figure 1 that will be used throughout this paper. Finally, a large part of the literature on Burst-Buffers shows how to use them with a specific application workflow [START_REF] Daley | Performance characterization of scientific workflows for the optimal use of burst buffers[END_REF][START_REF] Schenck | Evaluation and performance modeling of a burst buffer solution[END_REF]. Specifically, they consider systems where applications have dedicated and pre-allocated Burst-Buffers, and where the application can explicitly control its data transfers and the use of the Burst-Buffer. This must be done for each application and is very platform dependent. In practice, only few applications have the human-power to implement this. By opposition, our work is only architecture dependent and does not require any additional work from application developers. However, we believe our results can also be used by applications developers if they want to estimate the size of Burst-Buffers that they would need based on their application characteristics. In [START_REF] Daley | Performance characterization of scientific workflows for the optimal use of burst buffers[END_REF] and [START_REF] Da Silva | On the use of burst buffers for accelerating data-intensive scientific workflows[END_REF], the authors analyze workflows used in HPC system (CAMP and SWarp in [START_REF] Daley | Performance characterization of scientific workflows for the optimal use of burst buffers[END_REF], CyberShake in [START_REF] Da Silva | On the use of burst buffers for accelerating data-intensive scientific workflows[END_REF]) in order to model their accesses to the storage system and to identify opportunities to leverage the capabilities of the Burst Buffer of NERSC's Cori system based on Cray DataWarp [START_REF] Henseler | Architecture and design of cray datawarp[END_REF].

Models

In this section we present the models used in this work.

Machine model

We assume that we have a parallel platform made up of N identical unit-speed nodes, composed of the same number of identical processors. We model the long-term storage system (Parallel File System or PFS) as a single file server with input bandwidth B (in a system with several file servers, B would represent their aggregate bandwidth). In addition to this file server, the platform is equipped with Burst-Buffers of input bandwidth B BB and size S. Figure 1 provides a schematic view of this model.

In this work, we assume that the Burst-Buffers can be partitioned either statically or dynamically between applications. In the case of static partitioning, once a new application enters the system, the scheduler decides on the share of the buffer allocated to this application. It cannot be modified as long as the application has not left the system. In the dynamic model, the share of each application can change at runtime. Formally, we write these two strategies:

1. Static: when an application A k enters the system, it is allocated a volume S k of the Burst-Buffers. S k remains constant throughout its execution, and must respect the following constraint: at any time, if {A k } k≤n are running on the system, then Note that part of this model has been verified experimentally to be consistent with the behavior of Intrepid and Mira, super-computers at Argonne [START_REF] Gainaru | Scheduling the i/o of hpc applications under congestion[END_REF] and Jupiter, a machine at Mellanox [START_REF] Aupy | Periodic i/o scheduling for super-computers[END_REF]. In order to be compliant with the strategies described in [START_REF]Trinity / NERSC-8 Use Case Scenarios[END_REF], we introduce (see Section 6 for details) the possibility of progressively providig Burst-Buffer resources to an application before it starts (to prefetch its input data) and to progressively remove Burst-Buffer resources from an application after it ends (to release processing resources as soon as possible).

Application Model

We consider scientific applications running simultaneously on a parallel platform. In the present study, there is no interaction with the batch scheduler and we assume that the set of processing resources provided to the application is given. With respect to I/Os, applications consist in a sequence of three consecutive (and possibly nil) actions:

1. Data fetching from disks (read); 2. computations (compute); and

data uploading on disks (write).

Formally, the application A k is released at time r k and consists of n k iterations. Iteration i ≤ n k of A k consists of three consecutive non-overlapping phases: a read phase, where R k,i denotes the volume of data read, at read bandwidth b r k ; a compute phase, where l k,i denotes the compute time; and a write phase, where W k,i denotes the volume of data to be written at write bandwidth b w k . We assume that the phases can not be overlapped for a given application: reading must be finished before the computation can start, and similarly the computation must be finished before starting to write. This constraint is representative of many applications, whose memory requirements prevent to fetch data for the next phase in advance if the data for the previous phase still occupies the memory. We however assume that the input data of the reading phases can be prefetched in a burst buffer if its size allows it: this data does not depend on the results of the previous computations. A more generic model taking data dependencies into account is out of the scope of this paper.

In practice, b r k and b w k depend on the resources allocated by the batch scheduler and are given for A k . Hence, an application can be written as:

A k = (r k , b r k , b w k , Π n k i=1 (R k,i , l k,i , W k,i )) . (1)
We also denote by

C min k = r k + n k i=1 R k,i b r k + l k,i + W k,i b w k (2)
the earliest an application could finish given its parameters and assuming that the system is not slowing it down. In practice this bound is hard to achieve on a machine: while the computations are done independently because each application uses its own nodes, the applications compete for sending and receiving data during their I/O phase on a dedicated I/O network, what results in congestions and delays between I/O nodes of the platform and the Parallel File System. We discuss the conditions necessary to reach this lower bound in Section 5.2. Execution Model In the execution of a schedule, there are notable events. Specifically, for each phase (R k,i , l k,i , W k,i ) of application A k we can define three events:

A k BBs PFS (f w k + o w k ) ≤ b w k f w k o w k o r k f w k + o r k A k BBs PFS (f r k + i r k ) ≤ b r k f r k i r k i w k f r k + i w k
1. The beginning of the read to the application that coincides with the end of the previous write from the application denoted t r k,i (with t r k,1 = r k ).

2. The beginning of the compute of the application that coincides with the end of the read to the application denoted t c k,i .

3. The beginning of the write from the application that coincides with the end of the computation phase denoted t w k,i .

We consider that above phases coincide without loss of generality, since the amount of I/O resources can vary with time and can be zero at the beginning and the end to encompass for delays. Finally, we denote by C k the end of the last write phase (coinciding with the end of the execution of A k ). For each application A k , the following set of functions (defined at each instant t) describe data movements (see Figure 2):

• Part of its output (write) data is sent to the PFS at rate f w k , and the rest to the Burstbuffers at rate o w k . This is a part of the W k,i phase.

• Part of its input (read) data is collected from the PFS at rate f r k , and the rest from the Burst-buffers at rate i r k . This is a part of the R k,i phase.

By definition, the following properties on

i r k + f r k (resp. o w k + f w k ) hold: 1. It is only non zero on the intervals [t r k,i , t c k,i ] (resp. [t w k,i , t r k,i+1 ]); 2. It is bounded by b r k (resp. b w k ); and 3. t c k,i t r k,i f r k (t)dt = R k,i (resp. t r k,i+1 t w k,i f w k (t)dt = W k,i ).
Independently of the current phase of the application, the buffer itself can prefetch or write data from/to the PFS. We denote by i w k and o r k the function of time expressing the rate at which this is done.

Inria

Optimization problem

We are now interested by the performance model of our applications. Let us consider application

A k = (r k , b r k , b w k , Π n k i=1 (R k,i , l k,i , W k,i
)) and let us first assume that it is running alone on the machine. In order to perform the I/O operations (R k,i ,W k,i ), several strategies can be used:

• Without burst-buffers, the I/O operations take a time of

R k,i min(B, b r k ) and W k,i min(B, b w k )
• With Burst-Buffers (and no size constraints), to execute R k,i , the buffer can prefetch the data at rate B while A k is performing its previous iterations. When A k is done with W k,i-1 , the data can be obtained in R k,i /b r k units of time. Similarly, to execute W k,i , the application sends the data on the buffer at rate b w k , then the buffer sends it on the PFS at rate B. With capacity constraints on Burst-Buffers, A k follows a mix of above two behaviors, depending on how much data can be stored into the buffers.

Finally, let us define the stretch of A k (s(A k )) in a schedule. Given C k , the end of the execution of A k in the schedule, and given C min k , the earliest date when A k may finish (as defined by Eq.( 2)), the stretch of A k is given by

s(A k ) = C k C min k . (3)
For both strategies X ∈ {Static, Dynamic}, we consider two different problems:

Definition 1 (X-Buffer-Size(ρ)). Find a schedule that minimizes the total size S of the Burst-Buffers with strategy X, for a maximum stretch of ρ (∀k, s(A k ) ≤ ρ).

Definition 2 (X-Stretch(S)). Find a schedule that minimizes the maximum stretch (max k s(A k )) with strategy X, where the total buffer size is bounded by S.

Dominant Schedules

A schedule is defined by the list of functions (

f w k , f r k , i w k , i r k , o w k , o r k )
, ∀k which describe the rates of data transfers. In general, the description of these functions over time is not polynomial in the size of the input problem. In this section, we prove that we can focus on strategies where each of these functions is constant between the different events of the schedule (as defined above). We call such schedules: Dominant Schedules. Hence, a schedule can be fully described by the set of events (t r k,i , t c k,i , t w k,i ) k,i , and the values of functions (

f w k , f r k , i w k , i r k , o w k , o r k
), ∀k at these events, what provides a polynomial size description whose correctness (with respect to resource limitations) can be checked in polynomial time.

Theorem 1. Given a schedule S = (f w k , f r k , i w k , i r k , o w k , o r k ) k , there exists a (dominant) schedule S = ( f w k , f r k , ĩ w k , ĩr k , õw k , õr k ) k such that:
• all applications have the same stretch as S;

• the total buffer size used is the same as S;

• between any two events(t r k,i , t c k,i t w k,i ) i,k of S, all functions f ∈ sched, are constant.

Proof. Let us denote by e 0 < • • • < e n the list of events of S. For f ∈ S, let us define f ∈ S by:

∀i, f : x ∈ [e i , e i+1 ] →
e i+1 e i f (t)dt ei+1-ei . This transformation satisfies the following constraints: • Application-specific constraints:

-All the data transfer necessary for a read/write phase is performed during those phases; -The maximum application bandwidths of S are never larger than b r k and b w k ; • PFS-specific constraint:

-The total PFS bandwidth is never larger than B;

• Buffer-specific constraints:

-The total buffer peak is never larger than S (X = Dynamic); -The sum of individual buffer peaks is never larger than S (X = Static).

-The data leaving the buffer is not larger than the data entering the buffer.

Let us first note that ∀f ∈ S, ∀j:

max [ej ,ej+1] f ≤ max [ej ,ej+1] f. ( 4
)
since f is the average of f on those intervals.

Application-specific constraints To check that the solution is valid, we need to check that for all R k,i , all necessary data for the read phase is actually read:

t c k,i t r k,i f r k (t) + ĩr k (t) dt = R k,i .
Denote j 0 and j 1 s.t. e j0 = t r k,i and e j1 = t c k,i .

ej 1 ej 0 f r k (t) + ĩr k (t) dt = j1-1 i=j0 ei+1 ei f r k (t) + ĩr k (t)dt = j1-1 i=j0 ei+1 ei ei+1 ei f r k (x) + i r k (x)dx e i+1 -e i dt = j1-1 i=j0 ei+1 ei f r k (x) + i r k (x)dx = R k,i
The last equality comes from the fact that S is a valid solution. Similarly we obtain that for W k,i we need

t r i+1 t w k,i f w k (t) + õw k (t)dt = W k,i .
In addition, from Eq. ( 4) and because S is a valid schedule, one can check that for all j:

max [ej ,ej+1] f r k + ĩr k ≤ max [ej ,ej+1] f r k + i r k ≤ b r k max [ej ,ej+1] f w k + õw k ≤ max [ej ,ej+1] f w k + o w k ≤ b w k Inria
PFS-specific constraint From Eq. ( 4), and because S is a valid schedule, one can verify that for all j:

max [ej ,ej+1] k f r k + ĩ w k + f w k + õr k ≤ max [ej ,ej+1] k f r k +i w k +f w k +o r k ≤ B
Buffer-specific constraint The occupation of the buffer at time x ≥ e 0 due to data from A k is given by

x e0 õw k (t) -õr k (t) + ĩ w k (t) -ĩr k (t) dt
For any j, let us define by S

k the volume of data in the buffer of A k at time e j :

S (j) k = ej e0 (o w k (t) -o r k (t)) + (i w k (t) -i r k (t)) dt.
Let x ∈ [e j0 , e j0+1 ]:

x e0 õw k (t) -õr k (t) + ĩ w k (t) -ĩr k (t) dt = j0-1 j=0 ej+1 ej õw k (t) -õr k (t) + ĩ w k (t) -ĩr k (t) dt + x ej 0 õw k (t) -õr k (t) + ĩ w k (t) -ĩr k (t) dt = j0-1 j=0 ej+1 ej ((o w k (t) -o r k (t)) + (i w k (t) -i r k (t))) dt + (x -e j0 ) ej 0 +1 ej 0 ((o w k (t) -o r k (t)) + (i w k (t) -i r k (t))) dt e j0+1 -e j0 = S (j0) k +
x -e j0 e j0+1 -e j0 S (j0+1) k -S (j0) k .

In the Dynamic case, one can check that ∀x ∈ [e j , e j+1 ],

k x e0 õw k (t) -õr k (t) + ĩ w k (t) -ĩr k (t) dt = k S (j) k + x -e j e j+1 -e j S (j+1) k -S (j) k ≤ max k S (j) k , k S (j+1) k ≤ B.
Similarly, in the Static case, for any application k, ∀x ∈ [e j , e j+1 ]:

x e0 õw k (t) -õr k (t) + ĩ w k (t) -ĩr k (t) dt ≤ max S (j) k , S (j+1) k ≤ S k .
With the same reasoning, one can also check that ∀x ≥ e 0 : ensuring that in the new schedule, the data does not leave the buffer before it enters it. This concludes the proof that the solution S is valid, and that its performance is identical to S (same time events, same buffer sizes).

In the future, we therefore restrict the search to Dominant Schedules.

Complexity Results

In this Section, we provide a NP-hardness proof for X-Stretch and X-Buffer-Size, using a reduction from the well-known 3-Partition problem (3-Part).

Theorem 2. X-Stretch(0) and X-Buffer-Size(ρ) for any fixed ρ such that 1 < ρ ≤ 2 are NP-complete.

The reduction we consider is adapted from [START_REF] Gainaru | Scheduling the i/o of hpc applications under congestion[END_REF] where it is applied to a the problem of scheduling periodic applications.

Proof. We consider the associated decision problem: given a set of K applications A k and a platform, does there exist a schedule of stretch at most ρ without using any buffer?

We have shown that by considering Dominant Schedules, the problem belongs to NP. We use a reduction from 3-Part. Consider an arbitrary instance I 1 of 3-Part: given an integer B and 3n integers a 1 , . . . , a 3n , s.t. 3n i=1 a i = nB, can we partition the 3n integers into n triplets

I 1 , • • • , I n , each of sum B?
We build the following instance I 2 of X-Stretch(0) and X-Buffer-Size(ρ): the maximum bandwidth of the I/O system is B, there are 3n applications released simultaneously (r k = 0) with one phase each

(n k = 1), R k,1 = W k,1 = a k , l k,1 = d, where d = n+1-2ρ ρ-1 . The maximum bandwidth of application A k is b r k = b w k = a k ,
so that for each application, its communication phase takes time at least 1. Note that C min k = d + 2. We study whether there exists a solution of buffer size S = 0 and with a stretch not greater than ρ (by definition of d, we have ρ = d+n+1 d+2 ), or equivalently, is there a schedule such that all applications finish before time d + n + 1. Note that the definition of d also ensures that d ≥ n -1 as long as ρ ≤ 2n n+1 (which holds if n is large enough).

We now prove that I 1 has a solution if and only if I 2 does.

Assume first that I 1 has a solution. Let us call I 1 , • • • , I n the n triplets of I 1 . By definition we have i∈It a i = B. We construct the followint solution for instance I 2 : if k ∈ I t , then R k,1 is scheduled from time t -1 to time t at maximum rate a k . Then W k,1 is scheduled from time t + d to time t + d + 1 at maximum rate a k (see Figure 3).

Inria

It is a valid solution for A k with respect to the read, compute and write constraints. The stretch of A k is t + d + 1/d + 2 ≤ ρ. Furthermore, since for all t i∈It a i = B, it is also a valid solution with respect to the IO bandwidth constraint.

Assume now that I 2 has a solution. By definition of the stretch, the latest an application can terminate is at time d + n + 1.

There cannot be any I/O movement between time n and d + 1:

• Write data are not ready yet: the minimum time needed is d + 1 units of time for any application;

• Read data should be over: the minimum time needed once data is read is d + 1.

Since the total I/O volume is k R k,1 + W k,1 = 2nB, then the I/O bandwidth must be used at full capacity from time 0 to time n, and from time d + 1 to time n + d + 1.

Lemma 1. In a solution to instance I 2 of length d + n + 1, if an application reads some data between time t and time t + 1, then its read phase finishes at time t + 1.

Proof. All read phases finish at time t ≤ n since it takes a minimum time of d + 1 to do the rest of the computation.

We show the result by contradiction. Let us assume that the claim is not true. Let t be the first time such that an I/O transfer occured between time t and t + 1 but did not complete by time t + 1.

Let V > 0 be the volume of this transfer from time t to t + 1. Since the total volume of I/O transferred from time 0 to t + 1 is at most B(t + 1), the amount of finished read phases is at most B(t + 1) -V .

Because of the structure of I 2 , the amount of write data available before time d + t + 2 is at most B(t + 1) -V . This contradicts the fact that the I/O bandwidth has been used at full capacity from time d + 1 to time d + t + 2, hence showing the result.

Let us consider any time interval [i, i + 1] for i ≤ n -1, then the communication link must be fully used, and communications must take time 1. This implies that the read phase of A k is performed at maximum rate a k . Therefore, partitioning the applications according to the interval in which their read phase happens provides a valid solution to the 3-Part problem. Theorem 3. For any S ≥ 0, Static-Stretch(S) is NP-complete.

Proof. We can extend the previous reduction from 3-Part with an additional application. With the same notations as above, we introduce another application A n+1 , with release time n, a single write phase of size S + Bρ, and a bandwidth b r k = S + Bρ. Furthermore, we ensure that d ≥ n + 1, which holds as long as ρ ≥ 2n+2 n+3 . The minimum execution time of A n+1 is C min n+1 = 1. To achieve stretch at most ρ, this application must therefore complete within time n + ρ. Since d ≥ n + 1, this happens between the read and the write phases of the other applications. However, during this time, only an amount Bρ of data can be sent to the PFS ; the rest of the write phase of A n+1 needs to be sent to the burst buffer. Thus in a solution of stretch ρ and burst buffer size S, all the buffer size needs to be allocated to A n+1 . In the Static model, this implies that no buffer remains available for the other applications, and the proof of the previous theorem allows to conclude.

Burst-Buffer Lower Bounds for the Execution of a Single Application

In Section 3, we have provided a lower bound on the execution time of a single application given its characteristics:

C min k = r k + n k i=1 R k,i b r k + l k,i + W k,i b w k .
In general this lower bound is not reachable. Indeed, to be able to reach it, A k needs to read and write at maximum bandwidth during its read and write phases. This is not typically doable for example if b r k > B r or b w k > B w . As noted in [START_REF]Trinity / NERSC-8 Use Case Scenarios[END_REF], Burst-Buffers are expected to accelerate applications by (i) accelerating the transfers to and from the Parallel File System by using the Burst-Buffers as a cache for writing (buffering) and reading (pre-fetching) data and (ii) enabling to reorganize the communications to the Parallel File System in order to avoid contention when accessing it. In this section, we focus on a single application running on the platform, and show how to optimally dimension the Burst-Buffers to minimize its runtime.

Description of a solution achieving optimal makespan

Let us consider application A k . Let C min k be the lower bound of C k , the makespan of A k as defined in Eq [START_REF]The trinity project[END_REF]. We study a solution that achieves the optimal makespan C min

k = r k + n k i=1 R k,i b r k + l k,i + W k,i b w k . time

Read

Figure 4: Data read by A k (when running in isolation) from the Parallel File System when assuming B = +∞ (solid line) and when using a Burst-Buffer (dashed line)

Let us consider the case where there are no constraints on the bandwidth to the Parallel File System, i.e. B = B r = B w = +∞. The minimal time for the read, compute and write phases are respectively R k,i /b r k , l k,i and W k,i /b w k . Since these phases cannot overlap, then for each iteration i:

1. t r k,i = r k + i-1 j=1 R k,j b r k + l k,j + W k,j b w k
. The read phase must take exactly R k,i /b r k units of time, hence it is performed at bandwidth b r k . This situation is depicted on the solid line in Fig. 4 and in what follows, we denote by R ∞ k (t) the value of the solid line at time t. The slope of R ∞ k (t) is therefore either 0 (during compute and write phases) or b r k during read phases.

Inria

2. t c k,i = r k + i-1 j=1 R k,j b r k + l k,j + W k,j b w k + R k,i b r k
. The compute phase lasts l k,i units of time.

3. t w k,i = r k + i-1 j=1 R k,j b r k + l k,j + W k,j b w k + R k,i b r k + l k,i
. This situation is similar to the read phase. We denote by W ∞ k (t) the value of the solid line at instant t (Fig. 5). The slope of W ∞ k (t) is therefore either 0 (during read and compute phases) or b w k during write phases.

Burst-Buffer size necessary to

A k to achieve C min k in isolation.
In the case where both b w k ≤ B w and b r k ≤ B r and if A k is running alone on the platform, then the solution that performs transfers as soon as possible and at maximal rate trivially achieves C min k and Burst-Buffer are only needed if several applications are running simultaneously, what will be considered in Section 6). This is not the case if the application can write (resp. read) at speed b w k > B w (resp. b r k > B r ) on the Burst-Buffer. In this case, in order to achieve C min k , A k must use Burst-Buffer resources, and our goal is this section is to find the minimal Burst-Buffer size, denoted as P k , so as to achieve such an exectution time. As it is generally the case in practice, we assume that A k read operations do not depend on previous write operations and can be stored from the Parallel File System to the Burst-Buffer in advance.

Let us first concentrate on the write operations onto the Parallel File System of an isolated application, in presence of a limited bandwidth B w to the PFS. The dashed line in Figure 5 depicts the volume of data written to the Parallel File System. Writing the data of the first writing phase to the PFS cannot start before time step r k + l k,1 + R k,1 /b r k and must last at least W k,1 /B w since B w is an upper bound on the achievable bandwidth to the PFS. In order to achieve makespan C min k , the second read phase must start at time

r k + l k,1 + R k,1 /b r k + W k,1 /b w k .
At that date, given the limited bandwidth B w to write data onto the Parallel File System, the Burst-Buffer is used to store the data that could be written to the Burst-Buffer (at rate b w k ) but not on the Parallel File System (at rate B w ). At any time step, the solid lie represents the overall volume of data sent by A k (either to the Burst-Buffer or the Parallel File System) and the difference between the solid and dashed plots represent the minimal volume of data that must be stored onto the Burst-Buffer. After time

r k + l k,1 + R k,1 /b r k + W k,1 /b w k ,
the Burst-Buffer is being emptied onto the Parallel File System, at rate B w following the model described in Section 3. This transfer stops either when the Burst-Buffer is empty (i.e. solid and dashed plots cross) or at time

r k + l k,1 + R k,1 /b r k + W k,1 /b w k + R k,2 + l k,2
(when a new write operation must start).

From then, and until time

r k + l k,1 + R k,1 /b r k + W k,1 /b w k + R k,2 + l k,2 + W k,2 /b w k
, the solid and dashed plots diverge again, meaning that the minimal amount of Burst-Buffer storage increases during this time interval at rate b w k -B w . Therefore, following this algorithm, it is possible to determine the minimal volume of Burst-Buffer for write operations that enables to process A k within the deadline C min k , and that corresponds to the maximal difference between the solid and dashed plots.

The situation for read operations is symmetric and is depicted in Figure 4, where the solid plot depicts the volume of data that must be read by the nodes running A k and the dashed plot depicts the volume of data that must be read from the Parallel File System (and sent either to A k nodes or prefetched on the Burst-Buffer). The algorithm to build the dashed plot is very similar to the algorithm described above in the case of Figure 5, and the necessary Burst-Buffer size for read operations to run A k in isolation with minimal makespan is given by the maximal difference between the plots on Figure 4.

In these derivations, we have considered that it is possible to fully use bandwidth B w (resp. B r ) when writing on (resp. reading from) the Parallel File System. Let us now consider the case where we add an additional constraint stating that the overall bandwidth (the aggregation of incoming and outgoing bandwidth) is bounded by B. Although this problem is more complicated, it can be solved in polynomial time by solving a Linear Program in rational numbers. Let us denote by e k l the ordered set of events (see Section 3 for a formal definition), i.e. instants where a read I/O phase, a processing phase or a write I/O phase starts (plus time e k 0 = r k ). According to Theorem 1, we only need to specify the amount of data transferred to and from the PFS at each of these events. Let us thus denote by

w PFS k (e k l ) (resp. r PFS k (e k l )
) the overall volume of data written by the application to the Parallel File System (resp. read from the Parallel File System either to the Burst-Buffer or the application) before time e k l . Let us also denote by Peak k the size of the Burst-Buffer necessary to handle both write and read operations to the Burst-Buffer.

Then, our goal is to minimize Peak k under the following constraints:

                   ∀l w PFS k (e k l ) ≤ W k (e k l ) ∀l r PFS k (e k l ) ≥ R k (e k l ) ∀l r PFS k (e k l+1 ) -r PFS k (e k l ) ≤ B r • (e k l+1 -e k l ) ∀l w PFS k (e k l+1 ) -w PFS k (e k l ) ≤ B w • (e k l+1 -e k l ) ∀l w PFS k (e k l+1 ) -w PFS k (e k l ) + r PFS k (e k l+1 ) -r PFS k (e k l ) ≤ B • (e k l+1 -e k l ) ∀l W k (e k l ) -w PFS k (e k l ) + r PFS k (e k l ) -R k (e k l ) ≤ Peak k
The first two constraints ensure that the amount of data written (resp. read) from the PFS is smaller (resp. larger) than the maximal (resp minimal) volume to achieve C min k . The following three constraints enforce that neither the read, nor the write and the overall bandwidths are exceeded. At last, the last constraint states that the Burst-Buffer size Peak k is enough to store both the amount of data written to the Burst-Buffer but not yet to the Parallel File System and the amount of data read from from the Parallel File System but not yet transmitted to the computation nodes.

Burst-Buffer Size to compensate for contention between multiple applications

We have seen in Section 4 that X-Buffer-Size(ρ) is NP-complete for 1 < ρ ≤ 2. In this section, we provide a polynomial-time algorithm to solve X-Buffer-Size(1), both in the Static Inria and Dynamic cases. It is computed via a Linear Program whose constraints are detailled in Section 6.2.

6.1 Data transfers from and to the Parallel File System in presence of a Burst-Buffer

Let us now compute the minimal Burst-Buffer size S * necessary to achieve completion time C min k , ∀k, even if all applications compete for the bandwidth to the Parallel File System. In this case, the Burst-Buffer will also be used to avoid contentions to the Parallel File System, by prefetching (before later reading) or storing (before later sending) data eventually to be read from or written to the Parallel File System.

The solid plot in Figure 6 is analogous to the solid plot in Figure 5 and depicts the evolution with time of the data volume that must be written by A k , either to the Parallel File System or to the Burst-Buffer, and is denoted by W ∞ k (t). On the other hand, the dashed plot in Figure 5 depicts the evolution of the data actually written to the Parallel File System in presence of a Burst-Buffer of size S w k , and is denoted by W S k k (t). Indeed, let us now consider the evolution with time of the minimal volume of data that can be sent to the Parallel File System in presence of a Burst-Buffer of size S w k , when achieving the same makespan C min k . At time C min k , all data produced by A k to be written on the Parallel File System must have been transferred either to the Burst-Buffer or to the Parallel File System, since the nodes allocated to A k must be released. Therefore, at least W k,i -S w k must have been written to the PFS, since the Burst-Buffer can hold at most S w k volume. The same holds at all the instants between

C min k -W k l k /b w k and C min
-W k l k /b w k -R k l k /b r k is at least l k -1 W k,i -S w k .
Using a trivial induction, we can prove that before instant C min k , the minimal volume sent to the Parallel File System is given by the dashed plot in Figure 6 and in what follows, we will denote by W S w k k (t) the value of the dashed plot at instant t. After time C min k , the data to be written to the Parallel File System that still reside in the Burst-Buffer must eventually be transferred to the Parallel File System to release space on the Burst-Buffer. The corresponding amount of Burst-Buffer storage is progressively released by A k during this transfer, and can be used by other applications. As previously stated, we search for a solution such that all tasks complete within time Ckmin but note that in our model, the Burst-Buffer is considered as a safe storage and may therefore hold data that should be eventually written to the Parallel File System. This emptying strategy is very similar to the ultimate draining strategy described in [START_REF]Trinity / NERSC-8 Use Case Scenarios[END_REF].

On the other hand, let us now consider any increasing function W k (t) whose plot remains between W S w k k (t) and W ∞ k (t) (i.e. the solid and dashed plots on Figure 6) and whose slope is always at most min(b w k , B w ). Then, W k (t) represents the volume of data written to the Parallel File System in a valid strategy that makes use of a Burst-Buffer of size at most S w k . Indeed, let us denote by W BB k (t) the volume written on the Burst-Buffer at time t, i.e.

W BB k (t) = W ∞ k (t) -W k (t). Then, since W k (t) is increasing and the slope of W ∞ k (t) is at most b w k , then the slope of W BB k (t)
is no more than b w k , that corresponds to the maximal bandwidth to the Burst-Buffer. Then, at any any instant the sum of the slopes of W BB k (that may be either positive, when the Burst-Buffer is filled, or negative, when the Burst-Buffer is emptied to the Parallel File System) and W k (t) is equal to the slope of W ∞ k (t), so that at any instant, the strategy is valid and transfers exactly the same volume of data from A k as for W ∞ k (t). Let us now consider the read phases of application A k . The situation without Burst-Buffer is depicted in the solid plot in Figure 7 (identical to the solid plot in Figure 4) and we will denote by R ∞ k (t) the corresponding value, that represents the minimal amount of data that must be read from the Parallel File System in order to achieve optimal makespan when there is no constraint on the bandwidth to the Parallel File System (B r = +∞). Then, the dashed plot, denoted as R S r k k (t) in what follows, represents the maximal amount of data that can be read from the PFS if A k benefits from a Burst-Buffer of size S r k . Then, as for write operations, any increasing function As already noted, we can consider two different settings for the partitioning of the Burst-Buffer.

R BB k (t) such that R ∞ k (t) ≤ R BB k (t) ≤ R
In the Static case, a buffer of size S k is allocated to A k during all its lifetime, whereas in the Dynamic case, the amount of Burst-Buffer allocated to A k can vary over time. In both cases, we assume that at any given time step, the volume S k of the Burst-Buffer allocated to A k can be shared between prefetch (size S r k ) and intermediate storage (size S w k ) usages. In what follows, Inria we will therefore rather consider a Burst-Buffer for A k of size S k that can be arbitrarily split between read and write operations. In the static case, we assume that the buffer allocated to A k is progressively released from time 0 (size 0) to time r k (size S k ) and is progressively removed from time C k (size S k ) to the end of the schedule, to be compliant with the model of [START_REF]Trinity / NERSC-8 Use Case Scenarios[END_REF].

Let us denote (see Section 3 for a formal definition) the set of points in Figures 6 and7 at which the slope of any function R ∞ k (t) or W ∞ k (t) can change. Let us first remark that all these events, and thus their relative ordering, do do not depend on the values of S r k and S w k , ∀k. As previously, we will denote by e l the time of the l-th event, i.e. e l ≤ e l+1 . In what follows, we define a set of linear constraints that must be satisfied by functions R ∞ k (t) and W ∞ k (t) when t ∈ {e l } and we rely on Theorem 1 to rebuild a valid strategy for R ∞ k (t) and W ∞ k (t), ∀t by using linear interpolation for all functions between the instants t ∈ {e l }.

Therefore, the following set of linear constraints S Static describes the set of valid solutions:

                                           ∀k, l R ∞ k (e l ) ≤ R BB k (e l ) ≤ R S r k k (e l ) ∀k, l W ∞ k (e l ) ≥ W BB k (e l ) ≥ W S r k k (e l ) ∀k, l 0 ≤ R BB k (e l ) ≤ R BB k (e l+1 ) ∀k, l 0 ≤ W BB k (e l ) ≤ W BB k (e l+1 ) ∀l, k (R BB k (e l+1 ) -R BB k (e l )) ≤ B r (e l+1 -e l ) ∀l, k (W BB k (e l+1 ) -W BB k (e l )) ≤ B w (e l+1 -e l ) ∀k, l, W BB k (e l ) ≥ W ∞ k (e l ) -S w k (e l ), ∀k, l, R BB k (e l ) ≤ R ∞ k (e l ) + S r k (e l ), ∀k, l, S w k (e l ) + S r k (e l ) = S l k , ∀k, l ∈ Z k , S l k ≤ S k , ∀l,
k S l k ≤ S In the last constraints of S Static , Z k denotes all the events in the interval [r k , C min k ], i.e. when processing nodes are actually allocated to A k . During this interval the amount of Burst-Buffer allocated to A k is exactly S k . On the other hand, as previously stated, the Burst-Buffer is progressively allocated to A k before r k and progressively released from A k after C min k . We can use S Static to optimize FindOptimalSize Static : Minimize S under the constraints S Static .

Linear Program to Compute the Optimal Burst-Buffer Size (Dynamic Case)

The linear program to compute the optimal Burst-Buffer size in the Dynamic case is very similar. We obtain S Dynamic by removing from S Static the constraint ∀k, l ∈ Z k , S l k ≤ S k , that states that the size of the Burst-Buffer allocated to A k cannot change and remains equal to S k when nodes are actually allocated to A k (i.e. ∀l ∈ Z k ). This leads to optimization problems FindOptimalSize Dynamic : Minimize S under the constraints S Dynamic .

Simulation Results

In this Section, we report extensive simulations to evaluate our linear programming formulations, and compare them with a classic fair sharing approach. In order to do so we instantiate our evaluation based on characteristics of the Intrepid platform, and based on a set of applications as described in our previous work [START_REF] Aupy | What size should your buffers to disks be?[END_REF]

Setup

We consider a set of 4 applications described in the APEX report [START_REF] Lanl | [END_REF] which represent the majority of the load at the LANL. The characteristics of these applications are provided in Table 1. We simulate the execution of these applications on a platform similar to the Intrepid Blue Gene/P supercomputer, used by the Argonne National Laboratory between 2008 and 2014, which was ranked number 3 on the June 2008 Top 500 list. This platform has 96,000 cores, the bandwidth to the file system is B = 160GB/s, and the bandwidth per core is b = 0.02GB/s. We assume that most of the I/Os of these applications come from periodic checkpoints. We estimate the checkpointing period using the checkpoint optimal period given by P = 2 * C * µ # cores , following [START_REF] Aupy | Checkpointing algorithms and fault prediction[END_REF]. In this formula, C denotes the checkpointing duration and µ denotes the MTBF (Mean Time Between Failure) of the individual nodes of the platform. In the simulations, we consider different possible values for the MTBF, ranging from 5 years to 50 years.

To build the actual workload trace, we select a set of 30 applications, where each application is picked from the four application models described in Table 1, with a probability proportional to its usage ratio as reported in [START_REF] Lanl | [END_REF] (Frequency in Table 1). These applications are scheduled in FIFO order on the cores, what provides starting and ending times for each application. In order to compare results, we consider several target values for the IO load of the applications (namely 20%, 50% and 80%), defined in the following way. Let us first remark that an application with checkpoint size s and period P induces an average bandwidth load of s P over the course of its execution. Then, when an application starts or ends, the total required bandwidth is updated, and the maximum value over time (normalized by B) provides the IO need induced by running applications. Once this IO need is computed, the checkpoint sizes of all the applications are multiplied by a constant factor (and the checkpointing periods are adequately recomputed) so as to obtain the targeted IO load. To number of cores needed by each application can be read in Table 1 and it can be used to determine the maximal bandwidth b k at which a given application can communicate with the Burst-Buffer and the Parallel File System (see Section 3). To model the processing time between two checkpoints, we add to the checkpointing period an additional 15% variability.

Therefore, using the values from the APEX data set [START_REF] Lanl | [END_REF] and summarized in Table 1 and the description of the Intrepid platform, we are able to instantiate all platform and application characteristics needed by the linear programs of Section 6. By changing the MTBF (from 5 years to 50 years) and by scaling checkpoint sizes, we are moreover able to study different hardware characteristics and different system loads, while keeping realistic application characteristics. In turn, the linear programs of Section 6 compute an optimal buffer size S OPT (both in static and dynamic cases) and its partitioning among applications (in the static case) to process all the applications with a maximum stretch of 1, i.e. as if the bandwidth to the Parallel File System was enough to cope with all transfers at any time.

In order to evaluate the influence of the size of the Burst- several Burst-Buffer sizes, ranging from 0 × S OPT to 3 × S OPT . To compare the results of the optimal solutions computed in Section 6 to what could be achieved using a dynamic system level strategy, we introduce a greedy strategy, , that shares the bandwidth to the Parallel File System in the following way: an application is write-active if it is in a write phase or it has output data in its buffer, and the write bandwidth is shared equally between all write-active applications (in the limit of their own output rate). A similar policy is used for sharing the read bandwidth. Additionally, we assume that the Burst-Buffer (of size 0 ≤ S ≤ S OPT ) is partitioned between applications proportionally to their respective share in the optimal solution computed by the LP. We then compute the maximum and average stretch for the applications, to evaluate the slowdown incurred by the non-optimal bandwidth sharing and sub-optimal Burst-Buffer size.

Results

The performance of the above greedy system level approach, that shares equally the bandwidth to the Parallel File System in case of conflicts to the resource access but partitions the Burst-Buffer as in the solution of the LP, is depicted in Figure 8 for both the Max Stretch and Average Stretch metrics, for different MTBF values (as noticed before, a larger MTBF induces rarer and larger checkpoints) and different system loads (either 20%, 50% or 80%). All points in Figure 8 corresponds to 10 executions with different sets of applications, where the length of the processing phases between checkpoints vary from an execution to another (with a 15% maximal variation), and the darker area shows the typical variability of the results. As expected, the stretch decreases when the Burst-Buffer increases, but only up to a certain size. If a size of 1 × S OPT is enough to achieve a solution with stretch 1 for all applications when using the solution of the LP proposed in Section 6, a size of 1.5 × S OPT is in general required to achieve the lowest stretch using the greedy strategy. Moreover, in particular for large MTBF (thus larger and rarer checkpoints) and load values, the limit value for the stretch is larger than 1, up to 1.15 for the maximum stretch when the load is 80% and the MTBF is 50 years. Nevertheless, even for high loads, the greedy strategy (with the optimal partitioning computed by the LP) is able to achieve close to optimal results when the MTBF is smaller than 10 years, what is in general considered as a reasonable assumption. At last, since the stretch is much higher than 1 for any buffer size below 1 × S OPT , these plots show that the value S OPT returned by the LP is crucial in order to set the size of the Burst-Buffer to a value that does not induce large stretch values while limiting the hardware cost.

The cost induced by opting for a static partitioning is analyzed in Table 2, that displays the ratio between the Burst-Buffer size required when using a static partitioning and the Burst-Buffer size required when using a dynamic partitioning, for different values of the MTBF (from 5 to 50 years) and for different load levels (from 20 to 80%). These results show that a static partitioning of the Burst-Buffer induces an overhead in size of 25 to 40% for most settings. This overhead is to be compared to the increased simplicity of deployment, in particular with respect to security and dynamic management issues.

Conclusion

We consider the problem of sizing and partitioning Burst-Buffers in the context of HPC and Data Science applications running on a supercomputer and competing for the access to the Parallel File System. In this context, our goal is to minimize the stretch experienced by the applications. Given the characteristics of the platform and of the application, we first prove a negative result stating that the problem of minimizing the stretch given a Burst-Buffer size is in general NP-Complete. Nevertheless, we also prove that the special case, of clear practical interest, where the goal is to find the minimal Burst-Buffer size and its partitioning between applications that are able to compensate both the limited bandwidth to the Parallel File System and the contention in the access to it, can be solved in polynomial time. At last, we prove that it is possible to derive from this optimal solution a simple dynamic strategy, that can be easily implemented at runtime, and that is able to achieve low stretches for most settings. Our study also enables to precisely assess the cost of partitioning the Burst-Buffer between applications as opposed to using it as non-dedicated resource shared between all applications. This work opens several important perspectives. First, if the behavior is well understood when the Burst-Buffer is large enough, the problem of finding efficient strategies when the Burst-Buffer is too small to achieve an optimal stretch is still open. Along the same direction, it would be of great practical importance to be able to assess the good behavior of the dynamic strategy based on the optimal solution of optimal stretch, both in theory and through a larger set of experiments. Other research directions include extending the model for a more precise data management: taking into account data reuse throughout the execution, and/or considering temporary checkpoint data that could remain on the burst buffer until the next checkpoint if space allows it, instead of being written to the PFS.

Figure 1 :

 1 Figure 1: Modeling of the pseudo-centralized platform.

n

  k=1 S k (t) ≤ S.

Figure 2 :

 2 Figure 2: Schematics of the bandwidth used for output (top) and input (bottom).

Figure 3 :

 3 Figure 3: Communication schedule obtained from a positive 3-Part instance

3 Figure 5 :

 35 Figure 5: Data written by A k (when running in isolation) onto the Parallel File System when assuming B = +∞ (solid line) and when using a Burst-Buffer (dashed line)

k

  and the volume transferred to the Parallel File System between instants C min k

Figure 6 :

 6 Figure 6: Data transfer to the Parallel File System (solid line without contention, dashed line with a burst buffer).

kFigure 7 :

 7 Figure 7: Data transfer from the Parallel File System (solid line without contention, dashed line with a burst buffer).

6. 2

 2 Linear Program to Compute the Optimal Burst-Buffer size (Static Case)

Figure 8 :

 8 Figure 8: Maximum stretch results (top) and average stretch results (bottom) of fair sharing when buffer size varies, for different MTBF and load values.

Inria

  Dynamic: at any time t, if A k is running on the system it can use a volume S k (t) of the Burst-Buffers. The same constraint on the total buffer used holds: at all time, if {A k } k≤n are running on the system, then

n k=1 S k ≤ S.

Inria

2.

Table 1 :

 1 . Characteristics of the applications in APEX data set.

	Workflow	EAP LAP Silverton VPIC
	Frequency	65	21	8	6
	Number of cores (thousands)	16	4	32	30
	Checkpoint size (GB)	3,200 2,000	44,800	3,750
	Typical Walltime (hours)	16	4	32	30

Table 2 :

 2 Buffer, in what follows, we consider Ratio of static to dynamic buffer sizes, for different MTBF and load values.

	Inria

RR n°9213

avenue de la Vieille Tour 33405 Talence Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgments

This work was partially supported in part by the French National Research Agency (ANR) in the frame of DASH (ANR-17-CE25-0004), and the "Investments for the future" Program IdEx Bordeaux -SysNum (ANR-10-IDEX-03-02). We would like to thank Gael Goret and his colleagues from the Data Management team at Bull/ATOS Grenoble.