Alexandre Boritchev

Daniel Eceizabarrena
email: deceizabarrena@bcamath.org

Vilaça Victor

Da

Victor Vilaça

Da A Rocha
email: vrocha3@gatech.edu

Riemann

RIEMANN'S NON-DIFFERENTIABLE FUNCTION IS INTERMITTENT

come

Introduction and motivation

Riemann's non-differentiable function

f (t) = ∞ n=1 sin n 2 t n 2 , t ∈ [0, 2π].
is a celebrated example of a continuous but almost nowhere differentiable function. It was introduced by the eponymous mathematician in the 1860s [START_REF] Weierstrass | Über continuirliche functionen eines reellen arguments, die für keinen werth des letzteren einen bestimmten differentialquotienten besitzen[END_REF], and since then it has been widely studied from an analytic perspective [START_REF] Duistermaat | Selfsimilarity of Riemann's nondifferentiable function[END_REF][START_REF] Gerver | The differentiability of the Riemann function at certain rational multiples of π[END_REF][START_REF] Gerver | More on the differentiability of the Riemann function[END_REF][START_REF] Hardy | Weierstrass' non-differentiable function[END_REF][START_REF] Jaffard | The spectrum of singularities of Riemann's function[END_REF]. Interesting by itself as one of the first examples of an analytically pathological function, there is evidence that it appears naturally in a physical context. Indeed, according to [START_REF] De La Hoz | Vortex filament equation for a regular polygon[END_REF], its complex-valued analogue

R(t) = ∞ n=1 e 2πin 2 t n 2 , t ∈ [0, 1], (1)
seems to be a fair representative of some particular temporal trajectories in the dynamics of vortex filaments. These are described by the Vortex Filament equation, also known as the Binormal Flow equation, a nonlinear PDE strongly linked to the nonlinear Schrödinger equation. These dynamics are associated to turbulent phenomena, so the study of properties of R that are related to turbulence is quite natural. The objective of our work is to give sharp quantitative estimates that measure one of these, the intermittency of R, following the ideas of Frisch [START_REF] Frisch | Turbulence: The Legacy of A.N. Kolmogorov[END_REF]Chapter 8].

In Subsections 1.1 to 1.3, we introduce the physical setting where Riemann's non-differentiable function appears naturally, as well as its relationship with the concept of intermittency in the context of turbulence. After that, in Subsection 2.1 we set notation, and in Subsection 2.2 we define two quantities measuring intermittency. One of them is based on frequency high-pass filtering; the other one makes use of the so-called structure functions, a popular tool in the study of turbulence.

Our main result is the following:

Theorem. Riemann's non-differentiable function [START_REF] Arms | Localized-induction concept on a curved vortex and motion of an elliptic vortex ring[END_REF] is intermittent from the point of view of both high-pass filtering and structure functions.

In Subsection 2.2, a quantitative and sharp version of this result is given in Theorem 1.

X t = X x ∧ X xx , (t, x) ∈ R × R, (VFE)
where X = X(t, x) ∈ R 3 is a curve parametrized by arclength x and time t. Here, ∧ is the usual cross product. Considering the Frenet-Serret frame (T , N , B), where T = X x is the tangent vector, N = T x / T x is the normal vector and B = T ∧ N is the binormal vector, the equation can be equivalently written as

X t = κ B.
This is why (VFE) is also called Binormal Flow equation. Here, κ(t, x) is the curvature which, together with the torsion τ (t, x), satisfies the Frenet-Serret relations

   T N B    x =    0 κ 0 -κ 0 τ 0 -τ 0       T N B    .
The equation (VFE), first obtained by Da Rios in 1906 [START_REF] Da Rios | Sul moto d'un liquido indefinito con un filetto vorticoso di forma qualunque[END_REF], and rederived in 1965 by Arms and Hama [START_REF] Arms | Localized-induction concept on a curved vortex and motion of an elliptic vortex ring[END_REF] from the Euler equation, is an asymptotic model for the dynamics of a space curve: the vortex filament. A remarkable result about this equation was given by Hasimoto in 1972 [START_REF] Hasimoto | A soliton on a vortex filament[END_REF], where he proved that the transformation

Ψ(t, x) = κ(t, x) e i x 0 τ (t,σ) dσ (2)
solves the nonlinear Schrödinger equation

iΨ t + Ψ xx + 1 2 |Ψ| 2 + A(t) Ψ = 0. (NLS)
Here, A is a real, time dependent function that depends on κ, τ and their derivatives (see [START_REF] Banica | Selfsimilar solutions of the binormal flow and their stability[END_REF]). This way, the Hasimoto transformation defines a correspondence between (VFE), the equation of the original problem, and the well-known nonlinear Schrödinger equation. This supplies a method to find solutions of (VFE), since they can be produced from particular solutions to (NLS), as explained in [START_REF] Banica | Selfsimilar solutions of the binormal flow and their stability[END_REF].

Self-similar solutions, studied by Gutiérrez, Rivas and Vega in [START_REF] Gutiérrez | Formation of singularities and self-similar vortex motion under the localized induction approximation[END_REF], and later by Banica and Vega in [START_REF] Banica | Selfsimilar solutions of the binormal flow and their stability[END_REF], are of the form

X(t, x) = √ t G x/ √ t , t > 0, x ∈ R,
where G is a regular function defined on R, obtained from the tangent vector of X in the Frenet-Serret system. These solutions are regular for t > 0, but they develop a singularity in the shape of a corner at time 0, which means that X(0, x) represents the union of two non-parallel half-lines, a V-shaped initial datum. This is the behaviour of filaments of air in a delta wing during a flight [START_REF] De La Hoz | A numerical study of the self-similar solutions of the Schrödinger map[END_REF]. Also, a corner can be created after a reconnection of two different filaments, as can be observed in the rear of a plane or in the study of superfluid helium [START_REF] Schwarz | Three-dimensional vortex dynamics in superfluid 4 He: Line-line and line-boundary interactions[END_REF]. In [START_REF] Banica | The initial value problem for the binormal flow with rough data[END_REF], Banica and Vega generalise the results for the initial V-shaped datum to any curve which is regular except at one point where it has a corner, and they show existence and uniqueness of solutions under suitable conditions. Later, in [START_REF] Banica | Evolution of polygonal lines by the binormal flow[END_REF], they analyse initial data given by polygonal lines, extending the one-corner problem to a many-corner one, where they also prove existence and uniqueness.

We are interested in a similar generalisation, when the initial datum is given by a closed, regular and planar M -sided polygon, studied in [START_REF] De La Hoz | Vortex filament equation for a regular polygon[END_REF]. This means that τ (0, x) = 0, so the initial datum in (NLS) is ψ(0, x) = κ(0, x). Through the Hasimoto transformation (2), a corner-shaped singularity for the Vortex Filament equation (VFE) becomes a Dirac delta in the curvature, so the initial datum is given by a periodic sum of Dirac deltas,

Ψ M (0, x) = k∈Z δ x -2π k M , x ∈ R. (3)
In [START_REF] De La Hoz | Vortex filament equation for a regular polygon[END_REF], the Galilean invariance of (NLS) is used to determine

Ψ M (t, x) = Ψ M (t, 0) k∈Z e -i(M k) 2 t+iM kx .
Moreover, evaluation at scaled rational times t p,q = (2π/M 2 )(p/q), where p, q ∈ Z are coprime, gives

Ψ M (t p,q , x) = 2π M q Ψ M (t p,q , 0) k∈Z q-1 m=0 G(-p, m, q) δ x - 2πk M - 2πm M q , ( 4
)
where G stands for the generalised quadratic Gauss sums

G(a, b, c) = c-1 m=0 e 2πi(am 2 +bm)/c , a, b ∈ Z, c ∈ N.
The solution (4) is a mathematical representation of the Talbot effect, a celebrated optical phenomenon discovered in 1836 [START_REF] Talbot | Facts relating to optical science[END_REF] and rigorously described for the first time in 1881 [START_REF] Lord Rayleigh | On copying diffraction gratings, and on some phenomena connected therewith[END_REF]. This effect is based on the diffraction of waves, and it happens when light crosses a diffraction grating, a plate with uniformly distributed parallel vertical slits. Interference creates narrow copies of the grating in every rational multiple of a fixed distance, where the pattern is reproduced exactly. Applications have been discovered in several fields such as imaging or lithography, and accordingly, further theoretical study has been carried out (see, for instance, [START_REF] Berry | Integer, fractional and fractal Talbot effects[END_REF][START_REF] Matsutani | Wave-particle complementarity and reciprocity of Gauss sums on Talbot effects[END_REF]). It is visually represented by a Talbot carpet [START_REF] Berry | Quantum carpets, carpets of light[END_REF].

The structure of the Talbot effect matches the evolution of the polygonal filament if slits of the grating and corners of the filament are identified. Indeed, since according to the Hasimoto transformation Ψ M shares support with the curvature of X M , from (4) one can deduce that the corresponding filament X M (t p,q , x) is an M q-sided, not necessarily planar polygon. To see this, we must recall that the filament is first defined in (0, 2π) and then extended periodically to R. In (4), a value of k corresponds to the interval (2πk/M, 2π(k + 1)/M ), and when k is fixed, the sum in m goes from 0 to q -1 because the value q corresponds to the interval with k + 1 and m = 0. Since (0, 2π) is covered by the values k = 0, 1, . . . , M -1, and for each of them every m = 0, 1, . . . , q -1 corresponds to a Dirac delta, the filament has M q corners at time t p,q . The resulting polygon is in general not planar.

Also in [START_REF] De La Hoz | Vortex filament equation for a regular polygon[END_REF], the authors consider the temporal trajectories of the corners of the initial filament, represented by X M (t, 2πk/M ) for every fixed k ∈ Z. However, they only study X M (t, 0), since the translation invariance of (NLS) together with the space periodicity of the initial datum (3) imply that all of them have the same structure. Then, they show numerically that X M (t, 0) is extremely similar to the image of the function

φ(t) = k∈Z e -4π 2 ik 2 t -1 -4π 2 k 2 ,
and after a proper rescaling depending on the number of initial sides M , it seems to converge to φ as M → ∞. Therefore, φ is intrinsically related to the evolution of vortex filaments following the binormal flow, and its geometric importance as a representative of a physical trajectory is highlighted. This function is a slight modification of Riemann's non-differentiable function (1), since

φ(t) = - i 2π 2 R(-2πt) + it + 1 12
.

This gives evidence of an intrinsic physical and geometric nature of Riemann's non-differentiable function.

Some geometric results were obtained in [START_REF] Duistermaat | Selfsimilarity of Riemann's nondifferentiable function[END_REF] and also in [START_REF] Chamizo | Differentiability and dimension of some fractal Fourier series[END_REF], where the box-counting dimension of some trajectories related to R and other similar Fourier series were analysed. Further results regarding φ were obtained by the second author in [START_REF] Eceizabarrena | Some geometric properties of Riemann's non-differentiable function[END_REF][START_REF] Eceizabarrena | Asymptotic behaviour and Hausdorff dimension of Riemann's non-differentiable function[END_REF][START_REF] Eceizabarrena | Geometric differentiability of Riemann's non-differentiable function[END_REF]. More precisely, in [START_REF] Eceizabarrena | Asymptotic behaviour and Hausdorff dimension of Riemann's non-differentiable function[END_REF] the Hausdorff dimension of the trajectory represented by its image is analysed, and in [START_REF] Eceizabarrena | Geometric differentiability of Riemann's non-differentiable function[END_REF] it is proved that this trajectory does not have a tangent anywhere. On the physical side, the fact that the trajectories of vortex filaments following the binormal flow are related to turbulent phenomena justifies our interest for the relationship between Riemann's non-differentiable function and turbulence-related properties, such as intermittency and multifractality.

1.2. Turbulence and intermittency. In a series of papers in 1941, Kolmogorov gave a precise description of the properties of turbulent flows, for which the Reynolds number is large enough. These articles, usually referred to as the K41 Theory, are seminal for the theory of turbulence. The idea is that at small scales and away from the boundaries, the behaviour of a turbulent flow is universal: it does not depend on the fluid under consideration, nor on the geometry of the particular mechanisms producing the turbulence. The K41 Theory assumes that the flows are homogeneous (statistically invariant under translations) and isotropic (statistically invariant under rotations). Among the results of these papers, K41 hypotheses imply the self-similarity of both the high-pass filters of the velocity and the velocity increments.

Although this theory gives remarkable results, such as the Kolmogorov four-fifths law [START_REF] Kolmogorov | Dissipation of energy in the locally isotropic turbulence[END_REF], two main issues have been pointed out: the lack of universality, emphasized by Landau in 1944 and reported later in [START_REF] Landau | Fluid mechanics[END_REF] (it was just a footnote in the first version of the book), and the lack of self-similarity for the flows at small scales (in the dissipation range), highlighted by several experiments. The intermittency of a flow is formally a measure of this lack of self-similarity, and several ways to define it can be found in the literature. In [16, Sections 8.2 and 8.3], using the probabilistic description of turbulence, Frisch gives two definitions of intermittency for random stationary functions.

The first idea is to deal with high-pass filters in the Fourier space. Let the velocity of the flow v be a stationary random function, which can be expressed in terms of its Fourier transform

v(t) = R e iωt v(ω) dω.
Deleting the low frequencies, the high-pass filtered signal v >Ω is defined by

v >Ω (t) = |ω|>Ω e iωt v(ω) dω, for Ω > 0.
This filter is used to define the flatness

F (Ω) = (v >Ω (t)) 4 (v >Ω (t)) 2 2 , (5)
a tool to measure the intermittency of a signal. Here • denotes an ensemble average, defined as the mean over many realizations of the flow with different initial conditions and forcing. This could as well be a temporal average if we assume the ergodicity of the flow, or a spatial one if we use its homogeneity. In the study of turbulence, these three definitions are usually considered to be equivalent. We remark that, due to the stationarity of the signal v, the flatness F does not depend on time.

The key point in this definition is that F is constant for self-similar flows. Indeed, if v(t) has self-similar high-pass filters with a scaling exponent h, then v >λΩ law = λ h v >Ω for any λ > 0. As a consequence, F (λΩ) = F (Ω) for every λ > 0. Therefore, the more the flatness F grows with Ω, the less v(t) is self-similar at small scales, and thus by definition the more the flow v(t) is intermittent.

Remark 1.1. Frisch [START_REF] Frisch | Turbulence: The Legacy of A.N. Kolmogorov[END_REF]Section 8.2] suggests that, in the setting above, the inverse of the flatness measures the time the studied signal is "on". This remark leads to a more visual definition of intermittency: a signal is said to be intermittent if it displays activity during only a fraction of time, which decreases with the scale under consideration.

The second idea for a definition of intermittency is to use the structure functions, which play a key role in the literature of turbulence (see [START_REF] Frisch | Turbulence: The Legacy of A.N. Kolmogorov[END_REF][START_REF] Jaffard | The spectrum of singularities of Riemann's function[END_REF] and also [START_REF] Boritchev | Sharp estimates for turbulence in white-forced generalised Burgers Equation[END_REF][START_REF] Boritchev | Decaying turbulence in the generalised Burgers Equation[END_REF][START_REF] Boritchev | Multidimensional potential Burgers turbulence[END_REF] for rigorous results for Burgers turbulence). Let v be the velocity of a turbulent flow. For r, l ∈ R 3 , we define the longitudinal velocity increment δv(r, l) between two points separated by l by

δv(r, l) = (v(r + l) -v(r)) • l l ,
For p ≥ 1 and > 0, we can now define the p th structure function S p by

S p ( ) = (δv( )) p , ( 6
)
where there is no dependence on r by homogeneity of the turbulence. Moreover, the above only depends on = l by isotropy. Then, the flatness with respect to structure functions, which we call G, is

G( ) = S 4 ( ) S 2 ( ) 2 = (δv( )) 4 (δv( )) 2 2 . ( 7
)
Once again, the key point is that G is constant for self-similar flows. Indeed, if δv( ) is self-similar with scaling exponent h, then δv(r, λ ) law = λ h δv(r, ) for any λ > 0. We obtain thus G(λ ) = G( ) for all λ > 0. Therefore, we have the same characterization as for the flatness F : the more G grows while goes to 0, the further v is from being self-similar at small scales, and thus by definition the more v is intermittent.

Remark 1.2. Through the hypothesis of the self-similarity of the velocity increments, the K41

Theory gives the exponent h = 1/3 and thus

S p ( ) ∼ →0 p 3 .
Several experiments, some of which are listed in [START_REF] Frisch | Turbulence: The Legacy of A.N. Kolmogorov[END_REF], show that the above does not hold when p = 3 (in particular when p is large). The Kolmogorov four-fifths law is concerned with the special case p = 3, treated without the hypothesis of self-similarity: assuming homogeneity, isotropy and finite dissipation, then

S 3 ( ) = - 4 5 ε ,
where ε is the mean energy dissipation per unit mass.

We mentioned that Riemann's non-differentiable function represents a trajectory of the filament in several experiments related to turbulent flows. The goal of this paper is thus to analyse the intermittency of Riemann's function adapting the two definitions of flatness given above.

1.3. Turbulence, multifractality and Riemann's non-differentiable function. Multifractality is a concept deeply related to fully developed turbulence in fluids. Experiments have shown that the velocity of a turbulent fluid follows a very erratic pattern, so a classification of its analytic regularity might offer relevant information about its evolution. Assuming the latter is described by some function f , such a classification is offered by the so-called spectrum of singularities, which we define in the following lines.

For α > 0, a function f is said to be locally α-Hölder regular in a point t 0 , which is denoted f ∈ C α (t 0 ), if there exists a polynomial P t 0 of degree at most α such that |f (t 0 + h) -P t 0 (h)| ≤ C|h| α for small enough h ∈ R. Let the Hölder exponent of f at t 0 be

H f (t 0 ) = sup {α : f ∈ C α (t 0 )} .
Then, the spectrum of singularities is defined for each α > 0 as the Hausdorff dimension of the set of points t 0 having Hölder exponent α,

d(α) = dim H {t 0 : H f (t 0 ) = α}.
The spectrum of singularities is difficult to measure in experiments, so in an attempt to overcome this, Frisch and Parisi [START_REF] Frisch | On the singularity structure of fully developed turbulence; appendix to Fully developped turbulence and intermittency[END_REF] conjectured a formula relating it with the structure functions [START_REF] Berry | Quantum carpets, carpets of light[END_REF]. In this deterministic context, Jaffard [START_REF] Jaffard | The spectrum of singularities of Riemann's function[END_REF][START_REF] Jaffard | The multifractal formalism for functions Part I: results valid for all functions[END_REF] adapts the definition of structure functions as

S p (f, ) = |f (t + ) -f (t)| p dt. ( 8
)
Unlike the structure functions from the K41 Theory, defined as moments of order p of a stationary random process, this definition uses L p norms, and thus absolute values in the integral. If one assumes that S p (f, ) scales like | | ζ(p) when → 0, the Frisch-Parisi conjecture states that the spectrum is the Legendre transform of ζ,

d(α) = inf p (αp -ζ(p) + 1) . ( 9
)
However, restrictions may arise in the range of validity of [START_REF] Boritchev | Multidimensional potential Burgers turbulence[END_REF]. A more general and complete version of the conjecture was given by Jaffard [START_REF] Jaffard | The multifractal formalism for functions Part I: results valid for all functions[END_REF], who showed that these restrictions can be avoided if

ζ(p) is substituted by η(p) = sup{s | f ∈ B s/p,∞ p
}, where B s/p,∞ p stands for the usual Besov spaces.

Indeed, η(p) generalises ζ(p).

There are examples when the Frisch-Parisi conjecture (9) does not hold, but it has been proved for self-similar functions [START_REF] Jaffard | The multifractal formalism for functions Part II: self-similar functions[END_REF]. In these cases, it is said that f satisfies the multifractal formalism. In [START_REF] Jaffard | The multifractal formalism for functions Part I: results valid for all functions[END_REF], partial results were obtained for generic functions.

Multifractal analysis of Riemann's non-differentiable function was carried out by Jaffard [START_REF] Jaffard | The spectrum of singularities of Riemann's function[END_REF] as a natural continuation of its regularity analysis. Indeed, after the differentiability problem was finally solved by Gerver in [START_REF] Gerver | The differentiability of the Riemann function at certain rational multiples of π[END_REF][START_REF] Gerver | More on the differentiability of the Riemann function[END_REF] following a partial result by Hardy in [START_REF] Hardy | Weierstrass' non-differentiable function[END_REF], the next natural step was to study the local Hölder regularity which, as explained above, is central in the determination of the spectrum of singularities. After an important initial contribution by Duistermaat [START_REF] Duistermaat | Selfsimilarity of Riemann's nondifferentiable function[END_REF], Jaffard proved that

d(α) =      4α -2, α ∈ [1/2, 3/4], 0, α = 3/2, -∞, otherwise,
showing that R is a multifractal function, since there is a continuum of Hölder exponents with non-trivial Hausdorff dimension. His proof is based on wavelet theory. He also checked the validity of the Frisch-Parisi conjecture ( 9), therefore proving that R satisfies the multifractal formalism. This suggested a turbulent nature of R, as was later observed in the context of the evolution of vortex filaments.

Statement of the result

2.1. Setting and notation. In this paper, we work with the circle T = R/Z, and unless otherwise stated we consider complex-valued functions defined on T. For p ≥ 1, we denote by L p the Lebesgue space on the circle, L p (T), and by p the Lebesgue sequence space p (Z). We work with functions such that the corresponding Fourier series

f (t) = n∈Z a n e n (t)
are absolutely convergent, where (e n ) n∈Z is the orthonormal basis of L 2 (T) defined by e n (t) = e 2πint and a n ∈ C are the Fourier coefficients of f . In particular, this implies that f is continuous, and therefore f ∈ L p (T) for every p ∈ [1, +∞].

In the case of Riemann's non-differentiable function, we use the notation

R(t) = ∞ n=1 e 2πin 2 t n 2 = ∞ k=1 σ k k e 2πikt ,
where σ k is defined by

σ k = 1, if k is the square of an integer, 0, otherwise.
For two positive functions f and g, we write f g to denote that there exists a constant C > 0 such that f ≤ Cg. We also write f g to denote that f g and g f . If the constants involved depend on some parameter α, we write f α g and f α g. 2.2. Flatness and main result. In the deterministic setting of Riemann's non-differentiable function, definitions ( 5) and ( 7) need to be modified. The standard way to do so is to substitute p-moments by the p-th powers of L p norms, as was done when going from ( 6) to ( 8).

Definition 2.1. For f : T → C and N ∈ N, we define the high-pass filter as the projection of f on Fourier modes above N , defined by

f ≥N (t) = |n|≥N a n e n (t), ( 10
)
where a n are the Fourier coefficients of f . Similarly, the low-pass filter is the projection of f on Fourier modes below N , defined by

f ≤N (t) = |n|≤N a n e n (t), (11)
The flatness of f in the sense of high-pass filtering is given by

F f (N ) = f ≥N 4 L 4 (T) f ≥N 4 L 2 (T)
.

We say that f is intermittent in the sense of high-pass filtering if lim

N →+∞ F f (N ) = +∞.
Remark 2.2. The filters [START_REF] Chamizo | Differentiability and dimension of some fractal Fourier series[END_REF] and [START_REF] Da Rios | Sul moto d'un liquido indefinito con un filetto vorticoso di forma qualunque[END_REF] can also be defined with strict inequalities analogously:

f >N (t) = |n|>N a n e n (t), f <N (t) = |n|<N a n e n (t).
S f,p ( ) = T |f (t + ) -f (t)| p dt. ( 12
)
The flatness of f in the sense of structure functions is

G f ( ) = S f,4 ( ) S f,2 ( ) 2 .
We say that f is intermittent in the sense of structure functions if

lim →0 G f ( ) = +∞.
Remark 2.4. If there is no risk of confusion regarding f , we write S p ( ) instead of S f,p ( ).

Remark 2.5. In [START_REF] Chamizo | Differentiability and dimension of some fractal Fourier series[END_REF], N measures small perturbations in high-frequency Fourier oscillations, while in [START_REF] Duistermaat | Selfsimilarity of Riemann's nondifferentiable function[END_REF], measures small variations in space. Both capture the small-scale behaviour of R. The result in Theorem 1 is consistent with this fact. Indeed, if we identify N with -1 , both F R (N ) and G R ( -1 ) measure the same phenomenon.

The asymptotics of both definitions of flatness are given by the following theorem:

Theorem 1. Let R be Riemann's non-differentiable function [START_REF] Arms | Localized-induction concept on a curved vortex and motion of an elliptic vortex ring[END_REF]. There exist N 0 ∈ N and 0 < 0 < 1 such that for N > N 0 and < 0 , we have

F R (N ) log N, G R ( ) log( -1 ).
Consequently, R is intermittent in the sense of both high-pass filtering and structure functions.

Remark 2.6. The set R([0, 1]) is not self-similar, but the asymptotic behaviour of R [START_REF] Duistermaat | Selfsimilarity of Riemann's nondifferentiable function[END_REF] and also Figure 1 reveal at least the presence of some approximate self-similar structure. Therefore, if intermittency is a measure of the lack of self-similarity, R should have weak intermittent properties, and its flatness should show this. The logarithmic growth of both F R and G R in Theorem 1 agrees with this interpretation.

In Section 3, we prove the high-pass filter part of Theorem 1 and in Section 4, we prove the structure function part.

Intermittency in the sense of high-pass filters

To prove the part of Theorem 1 concerning high-pass filters, we will use the Littlewood-Paley decomposition of R, as well as a result of Zalcwasser [START_REF] Zalcwasser | Sur les polynomes associés aux fonctions modulaires θ[END_REF] on the L 4 norm of the sum of square-phased exponentials. Both results are stated in Appendices A.1 and A.2.

We estimate the L 2 norm of the high-pass filter first.

Lemma 3.1. For every N ≥ 2, R ≥N L 2 (T) N -3/4 .
Proof. By Plancherel's theorem, we get

R ≥N 2 L 2 (T) = 1 0 ∞ n=N σ n n e 2πint 2 dt = ∞ n=N σ n n 2 = ∞ n= √ N 1 n 4 ∞ √ N dt t 4 = 1 N 3/2 .
To compute the L 4 -norm of the high-pass filter, one may try to use Plancherel's theorem again, since f 4 4 = f 2 2 2 holds. However,

R 2 ≥N (t) = ∞ k=2N k-N n=N σ n σ k-n n (k -n) e 2πikt , ( 13
)
whose Fourier coefficients are related to k-N n=N σ n σ k-n , the number of ways in which k can be written as a sum of two squares both of which are greater than N . The study of such sums is a classical problem in number theory and can be very technical. Instead, the Fourier series of R can be decomposed in frequency pieces that act almost independently by the Littlewood-Paley decomposition. We use this technique to prove the following lemma.

Lemma 3.2. There exists

N 0 ∈ N such that R ≥N L 4 (T) N -3/4 (log N ) 1/4 , ∀N ≥ N 0 .
Proof. According to Appendix A.1, the Littlewood-Paley decomposition of R ≥N is

R ≥N (t) = ∞ j=1 ∆ j R ≥N (t), (14)
where the Littlewood-Paley pieces are

∆ 0 R ≥N (t) = 1≤n<A σ n n e 2πint and ∆ j R ≥N (t) = A j ≤n<A j+1 σ n n e 2πint , ∀j ∈ N.
The value of A will be chosen later (see [START_REF] Hasimoto | A soliton on a vortex filament[END_REF]). We define j(N ) as the index corresponding to the piece containing the N -th Fourier coefficient, the only one satisfying A j(N ) ≤ N < A j(N )+1 . Then, ∆ j R ≥N = 0 for every j < j(N ). By the Littlewood-Paley theorem, we may write the inequality

R ≥N L 4 (T)   j≥j(N ) |∆ j R ≥N | 2   1/2 L 4 (T) ≥ ∆ i(N ) R L 4 (T) , (15)
where i(N ) = j(N ) + 1, using that ∆ i(N ) R = ∆ i(N ) R ≥N . The choice of i(N ) comes from the fact that it is the first complete Littewood-Paley piece after j(N ), which is truncated as a consequence of the high-pass filter.

Let us estimate ∆ i(N ) R L 4 (T) . As in [START_REF] Eceizabarrena | Some geometric properties of Riemann's non-differentiable function[END_REF], we use Plancherel's theorem to write

∆ i(N ) R 4 L 4 = (∆ i(N ) R) 2 2 L 2 =   A i(N )+1 n=A i(N ) σ n n e 2πint   2 2 L 2 2A i(N )+1 k=2A i(N ) n σ n σ k-n n (k -n) 2 , ( 16
)
where the index n must satisfy

A i(N ) ≤ n ≤ A i(N )+1 and A i(N ) ≤ k -n ≤ A i(N )+1
. In both cases, n A i(N ) and k -n A i(N ) . Hence, we can take the denominators outside the sum:

∆ i(N ) R 4 L 4 1 A 4i(N ) 2A i(N )+1 k=2A i(N ) n σ n σ k-n 2 = 1 A 4i(N ) A i(N )+1 n=A i(N ) σ n e 2πint 4 L 4 . ( 17
)
This is a sum of exponentials with squared phases, whose L p norms were computed by Zalcwasser [START_REF] Zalcwasser | Sur les polynomes associés aux fonctions modulaires θ[END_REF]. First, the triangle inequality gives

A (i(N )+1)/2 n=A i(N )/2 e 2πin 2 t L 4 A (i(N )+1)/2 n=1 e 2πin 2 t L 4 + A i(N )/2 n=1 e 2πin 2 t L 4 , ( 18
)
so by Zalcwasser's theorem in Appendix A.2 we get, for large enough N ,

∆ i(N ) R L 4 1 A i(N ) A i(N )+1 log A (i(N )+1)/2 1/4 + A i(N ) log A i(N )/2 1/4 1 A i(N ) A i(N ) log A i(N )/2 1/4 A -3i(N )/4 log A i(N ) 1/4 . ( 19
)
On the other hand, using the reverse triangle inequality in (18), we get

A (i(N )+1)/2 n=A i(N )/2 e 2πin 2 t L 4 A (i(N )+1)/2 n=1 e 2πin 2 t L 4 - A i(N )/2 n=1 e 2πin 2 t L 4 .
Let us denote the constants in Zalcwasser's theorem for p = 4 by 0 < c < C. We get

∆ i(N ) R L 4 1 A i(N ) c A i(N )+1 log A (i(N )+1)/2 1/4 -C A i(N ) log A i(N )/2 1/4 = 1 A i(N ) (A i(N ) log A i(N )/2 ) 1/4 (c A(1 + i(N ) -1 )) 1/4 -C 1/4 A -3i(N )/4 log A i(N ) 1/4 c 1/4 A 1/4 -C 1/4 . ( 20
)
Finally, choose A so that

(cA) 1/4 -C 1/4 = 1. (21)
Observe that in the proof of [START_REF] Gerver | The differentiability of the Riemann function at certain rational multiples of π[END_REF] and [START_REF] Gerver | More on the differentiability of the Riemann function[END_REF] we may replace i(N ) by any j ≥ i(N ), so we have proved that

∆ j R L 4 A -3j/4 log A j 1/4 , ∀j ≥ i(N ). (22)
Coming back to j = i(N ), from [START_REF] Eceizabarrena | Geometric differentiability of Riemann's non-differentiable function[END_REF], since A i(N )-1 ≤ N < A i(N ) , we get

R ≥N L 4 (T) N -3/4 (log N ) 1/4 .
Lemma 3.2 suffices to prove that the flatness of Riemann's non-differentiable function tends to infinity. However, we can be more precise and show that the lower bound in the lemma is sharp.

Lemma 3.3. There exists

N 0 ∈ N such that R ≥N L 4 (T) N -3/4 (log N ) 1/4 , ∀N ≥ N 0 .
Proof. Applying the triangle inequality in the Littlewood-Paley decomposition (14), we write

R ≥N L 4 ≤ ∆ j(N ) R ≥N L 4 + j≥i(N ) ∆ j R L 4 . ( 23
)
Using [START_REF] Hardy | Weierstrass' non-differentiable function[END_REF] we can estimate ∆ j R L 4 for j ≥ i(N ). To deal with the index j(N ), following the arguments in ( 16), ( 17), ( 18) and ( 19) and using A j(N ) ≤ N < A j(N )+1 , we write

∆ j(N ) R ≥N L 4 1 N A j(N )+1 n=N σ n e 2πint L 4 1 N   A j(N )+1 n=1 σ n e 2πint L 4 + N n=1 σ n e 2πint L 4   1 N (N log N ) 1/4 = N -3/4 (log N ) 1/4 . ( 24
)
On the other hand, using [START_REF] Hardy | Weierstrass' non-differentiable function[END_REF], we bound

j≥i(N ) ∆ j R L 4 j≥i(N ) A -3j/4 log A j 1/4 j≥i(N ) j 1/4 A -3j/4 . ( 25
)
By Hölder's inequality one can write j≥i(N )

j 1/4 A -3j/4 ≤   j≥i(N ) j A -3j/4   1/4   j≥i(N ) A -3j/4   3/4 . ( 26
)
The second sum is geometric and equals A -3i(N )/4 1 -A -3/4 A -3i(N )/4 . The first one can be computed differentiating power series. Indeed, for |r| < 1 we write j≥i(N ) j r j ≤ j≥i(N ) (j + 1) r j = d dr j≥i(N ) r j+1 = d dr

r i(N )+1 1 -r = i(N ) r i(N ) 1 -r 1 + 1 i(N ) 1 + r 1 -r ≤ 2 1 -r i(N ) r i(N ) .
The last inequality is satisfied when i(N ) > 1 + r/(1 -r), which holds when N > A 1+r/(1-r) . Choosing r = A -3/4 , from ( 25) and ( 26), we get j≥i(N ) ∆ j R L 4 i(N ) 1/4 A -3i(N )/4 N -3/4 (log N ) 1/4 .

Combining this with ( 23) and ( 24), we finally obtain

R ≥N L 4 N -3/4 (log N ) 1/4 ,
for large enough N .

From Lemmas 3.1, 3.2 and 3.3, the proof of the high-pass filter part of Theorem 1 is immediate. Consequently, with the same restrictions on m as above, Bounding the sine trivially and using [START_REF] Eceizabarrena | Some geometric properties of Riemann's non-differentiable function[END_REF] with N = -1 /2, we get ( -1 ) -3 log -1 3 log( -1 ).

I ≤ -1 /2 ( ) 2 2
I > -1 /2 ( ) 2 2 2 ≤ k>l -1   m, k-m> -1 /2 σ m σ k-m m(k -m)   2 = (R ≥ -1 /2 ) 2 2 2 = R ≥ -1 /2
Remark 4.3. The method above can be generalised to compute the asymptotic behaviour of every S 2k , k ∈ N.

Proof of Theorem 1 (Part 2)

. By Lemmas 4.1 and 4.2, for small enough we get G R ( ) 3 log( -1 ) ( 3/2 ) 2 = log( -1 ).

Figure 1 . 1 . 1 .

 111 Figure 1. The image of Riemann's non-differentiable function R.

For

  any a, b ∈ R such that a < b, we write b n=a = n∈[a,b]∩Z and b n>a = n∈(a,b]∩Z .

Definition 2 . 3 .

 23 Let p ≥ 1, f : T → C a bounded and measurable function and ∈ [0, 1]. The structure functions of f are defined by

Proof of Theorem 1 ( 4 L 2 k 2 m

 1422 Part 1). For large enough N determined by Lemmas 3.2 and 3.3, we may writeF R (N ) = R ≥N 4 L 4 (T) R ≥N (T) N -3 log N N -3 = log N,so we look for both upper and lower estimates for I ≤ -1 /2 ( ) 4 , but an upper bound suffices forI > -1 /2 ( ) 4 .By Parseval's theorem, we haveI ≤ -1 /2 ( ) ) sin(π(k -m) ) m(k -m) σ m σ k-m .where the index m satisfies 1≤ m ≤ -1 /2 and 1 ≤ k-m ≤ -1 /2.As above, since 1/2 ≤ sin(t)/t ≤ 1 for |t| ≤ π/2, we get b σ m σ k-m .

4 √- 1 / 2 k=1 e k 2 4 4

 4124 By Zalcwasser's theorem in Appendix A.2, for small enough we getI ≤ -1 /2 ( ) 4 4 = I ≤ -1 /2 () frequency piece, using again Parseval's theorem, we can writeI > -1 /2 ( ) -m> -1 /2 sin(πm ) sin(π(k -m) ) m(k -m) σ m σ k-m .

4 4 ,

 44 so by Lemma 3.3 we obtain thatI > -1 /2 ( ) 4 4

Acknowledgments

The authors would like to thank Valeria Banica and Luis Vega, who made this collaboration possible through funding and advice, and Francesco Fanelli, Evelyne Miot and Dario Vincenzi, for several useful discussions. AB and VVDR have been supported by the IUF grant of Valeria Banica. Moreover, they would like to acknowledge the support of ANR ISDEEC. DE has been supported by the Ministry of Education, Culture and Sport (Spain) under grant FPU15/03078. DE and VVDR have been supported by the ERCEA under the Advanced Grant 2014 669689 -HADE and also by the Basque Government through the BERC 2018-2021 program and by the Ministry of Science, Innovation and Universities: BCAM Severo Ochoa accreditation SEV-2017-0718. VVDR is also supported by NSF grant DMS 1800241.

and therefore lim N →∞ F R (N ) = +∞.

Intermittency in the sense of structure functions

We begin by observing that for t, ∈ [0, 1], making the elementary change of variables t → t -/2, the structure functions S p ( ) can be described in terms of the increment function I( , t),

so that if we write I( ) = I( , •), we have

Proof. By Parseval's theorem,

On the one hand, since sin(t

On the other hand, we have the upper bound

Lemma 4.2. There exists

Proof. We decompose I( ) in low and high frequencies so that by the triangle inequality we get

On the other hand, the Fourier coefficients of I( ) are positive for every 1 ≤ k ≤ -1 . This, together with Parseval's theorem, implies that

In short, we have Theorem 2. Let p > 1, A > 1 and f (t) = n∈Z a n e 2πint a function in L p (0, 1). Consider the decomposition

Then, there exist constants B 1 , B 2 > 0 depending on p such that

A.2. A theorem of Zalcwasser. The following result, contained in [START_REF] Zalcwasser | Sur les polynomes associés aux fonctions modulaires θ[END_REF], is crucial to our proof.


Ana Gainaru
email: ana.gainaru@vanderbilt.edu

Guillaume Pallez
email: guillaume.pallez@inria.fr

Making Speculative Scheduling Robust to Incomplete Data

Keywords: HPC scheduling, stochastic applications, performance modeling, discrete and continuous estimators, sampling

In this work, we study the robustness of Speculative Scheduling to data incompleteness. Speculative scheduling has allowed to incorporate future types of applications into the design of HPC schedulers, specifically applications whose runtime is not perfectly known but can be modeled with probability distributions. Preliminary studies show the importance of speculative scheduling in dealing with stochastic applications when the application runtime model is completely known. In this work we show how one can extract enough information even from incomplete behavioral data for a given HPC applications so that speculative scheduling still performs well. Specifically, we show that for synthetic runtimes who follow usual probability distributions such as truncated normal or exponential, we can extract enough data from as little as 10 previous runs, to be within 5% of the solution which has exact information. For real traces of applications, the performance with 10 data points varies with the applications (within 20% of the full-knowledge solution), but converges fast (5% with 100 previous samples).

Finally a side effect of this study is to show the importance of the theoretical results obtained on continuous probability distributions for speculative scheduling. Indeed, we observe that the solutions for such distributions are more robust to incomplete data than the solutions for discrete distributions.

I. MOTIVATION

With the drive to incorporate big data and machine learning applications to HPC, new types of applications need to be considered when submitting jobs using HPC schedulers. These applications have heterogeneous, dynamic and data-intensive requirements and their performance is widely dependent on the input data. The convergent nature of the code or the continuous updates to the used algorithms makes the prediction of their execution time extremely hard. Currently, the best solution to understanding the behavioral patterns of such applications is to use a model of their execution time (for instance a probability distribution) constructed based on historical data.

In our previous work [START_REF] Gainaru | Speculative scheduling for stochastic hpc applications[END_REF], we have shown that using the current HPC schedulers with applications whose execution time is not precisely known can lead to poor response time and low utilization. Indeed, many of such software systems are reservation-based (for instance Slurm [START_REF] Yoo | Slurm: Simple linux utility for resource management[END_REF], Torque [START_REF] Staples | Torque resource manager[END_REF] and Moab [START_REF] Capit | A batch scheduler with high level components[END_REF]), meaning that the schedule is computed thanks to an estimate (often an overestimate) given by the users. The possible gaps generated by the applications finishing early are then filled using backfilling. However, if one overestimates consistently the execution time of applications, the backfilling work is not enough to cover for the gaps generated. In addition, platform administrators tend to penalize overestimation by decreasing user priority for future submissions.

While one natural direction could be to develop more precise models to try to obtain precise performance estimation [START_REF] Carrington | A performance prediction framework for scientific applications[END_REF], [START_REF] Leo T Yang | Cross-platform performance prediction of parallel applications using partial execution[END_REF], we have chosen a different complementary approach.

In our work [START_REF] Gainaru | Speculative scheduling for stochastic hpc applications[END_REF], we have developed a proof of concept for a novel HPC scheduler that would integrate the uncertainty in the execution time of those new applications. The key idea is the following: in order to deal with uncertainty in the execution time, instead of doing a single reservation, the scheduler should schedule a series of increasing-size reservations (called a strategy): S = (t 1 , t 2 , • • • , t i , • • • ). For a given batch of jobs, it starts by scheduling the first reservations of all the jobs of the local batch. For each job of the batch, either the first reservation was sufficient to finish the job execution (in which case this job is considered executed), or it was not (in this case the execution is failed, the time considered wasted). For all the jobs where the first reservation was not large enough, the scheduler schedules the second reservation, and so on until all the jobs of the batch are executed. We have shown in a previous work [START_REF] Aupy | Reservation strategies for stochastic jobs[END_REF] how to derive the optimal sequence of reservations to minimize the expected cost of a job, given the knowledge of the its execution profile (modeled by a probability distribution). Using a modified version of this algorithm that incorporates backfilling, we also showed [START_REF] Gainaru | Speculative scheduling for stochastic hpc applications[END_REF] that the batch scheduler algorithm derived above was more efficient for both total system throughput, and average job response time! All these previous results rely on the assumption that we know the correct distribution of possible walltimes for a job. An example of such a profile could be: with an input of size 50GB, on 2048 nodes, the job's execution time follows a lognormal distribution of mean 10h and variance 4h. However, in practice, getting the correct application profile (even if it is a probability distribution) is complicated. One way is to use historical data (the last M runs with similar parameters) as a basis for this profile [START_REF] Leo T Yang | Cross-platform performance prediction of parallel applications using partial execution[END_REF]. Figure 1 shows the historic execution times for two stochastic and big data applications from the neuroscience field for the 2017 year together with two possible fits (polynomial or using distributions) of the data. In this work we study the robustness of our solutions under various considerations. We show how one could design a performance estimator in a scheduler to derive efficient reservation strategies, solely based on historical data. We show that the natural strategy of using the data as is (discrete model) and computing the optimal solution is not efficient, but that with a simple transformation of the data to a continuous model, the strategies are often much more efficient.

In this work we make the following contributions:

• We study different applications both synthetic and real to study the impact of incomplete data. • We show that even with minimal application knowledge, the strategies developed in our previous work still perform well. • We show the importance of the theoretical study for continuous distributions, even though most of the data used to compute an application profile is in discrete format.

• Finally, we give a high level intuition on how one would include this into a scheduler. The rest of the paper is organized as follows. Section II presents related works and an overview of the set of applications that show stochastic behavior and would benefit from the presented strategy. In Section III we remind the Speculative Scheduling model, and describe the optimization problem with incomplete knowledge of the application profile. In Section IV, we propose two different algorithms to derive strategies to solve the optimization problem. We then study it more on synthetic data in Section V, before experimenting on actual Fig. 2. The relation between the input data size and total execution time for a cortical surface modeling application traces of applications in Section VI. Finally, we give some pointers on how to best practices for running on HPC systems before giving some concluding remarks.

II. BACKGROUND AND RELATED WORK

Emerging applications normally have a stochastic nature due to two reasons: (i) users not having a full understanding of the runtime behavior and resource requirements prior to execution and (ii) applications changing their runtime behavior and resource requirements during execution [START_REF] Weidner | Rethinking high performance computing platforms: Challenges, opportunities and recommendations[END_REF]. In our experience many big data and emerging applications fit both criteria. They are driven by the dynamic nature of each emerging field, using adaptive algorithms that can have different resource requirements depending on input data characteristics (e.g. the proximity of the input data to the convergence value) or on input parameters. (e.g. a specific area of interest which triggers an in-depth analysis algorithms). Their exploratory codes are in continuous change depending on previous results and thus so is the dependency of their resource requirements to the input data. There are currently a few studies to look at the performance of specific variable applications (e.g. for seismic wave propagation [START_REF] Carrington | High-frequency simulations of global seismic wave propagation using specfem3d globe on 62k processors[END_REF], or for shock hydrodynamics [START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF]). There are also studies that analyze platform workload's performance (e.g. for systems at TACC [START_REF] Evans | Comprehensive resource use monitoring for hpc systems with tacc stats[END_REF] or at Sandia National Laboratory [START_REF] Agelastos | Toward rapid understanding of production hpc applications and systems[END_REF]). The study in [START_REF] Tuncer | Diagnosing performance variations in hpc applications using machine learning[END_REF] characterizes application variation due to resource contention and software-related problems. However, there is nothing in literature on understanding application behavioral patterns for applications like the ones developed by emerging fields that focus on productivity and not performance.

Emerging applications are typically composed of multiple stages that could use multiple programming languages and with different performance characteristics. The total walltime of an application could depend on complex properties inside the input data and do not have a correlation with the input data size (typical for current scientific large-scale applications). Figure II shows the normal pattern in emerging applications, where the execution time can vary several orders of magnitude for the same input data size.

Current HPC schedulers extensively used on today's largescale systems rely on accurate estimates for the requested walltimes and memory requirements. Enforcing static limits for resources leads to two unfavorable scenarios: (i) underestimation when the users request resource requirements lower than the actual needs of their application in which case the runtime system kills the execution the moment the needs of the application exceed the requested amount; or (ii) overestimation in which case the utilization of the system is being wasted and users are being penalized by the cluster policies (by larger wait times in the queue due to either the system having to schedule larger jobs than necessary or by system administrator decreasing the user priority). Datacenter platforms like UC Berkeley's Apache Mesos [START_REF] Hindman | Mesos: A platform for fine-grained resource sharing in the data center[END_REF] and Google's Borg [START_REF] Verma | Large-scale cluster management at google with borg[END_REF] alleviate some of the constraints of current HPC platforms for emerging applications. However they are built for Cloud type workflows, are disruptive in nature and do not allow for an evolution in their design. In addition, emerging applications, often have data which is tied down to a particular site and therefore the only viable option is to run the application on local institutional clusters.

Considerable research has been conducted to improve job runtime estimates. In [START_REF] Gaussier | Improving backfilling by using machine learning to predict running times[END_REF], the authors used a polynomial model to enhance the accuracy of runtime estimates; [START_REF] Mustafa | A machine learning approach for predicting execution time of spark jobs[END_REF] uses regression models to alleviate the problem of underestimation; [START_REF] Tsafrir | Backfilling using systemgenerated predictions rather than user runtime estimates[END_REF] uses the average of the last two runtimes as a prediction. While this last method is simple, it is capable of doubling the accuracy. Based on current research, the Medicalimage Analysis and Statistical Interpretation (MASI) [START_REF] Masi Lab | Medical-image Analysis and Statistical Interpretation (MASI) Lab[END_REF] laboratory at Vanderbilt is using a similar method, averaging the last 10 runs of an application for determining the request time, and doubling it each time it underestimates the walltime. This method adapts to the changes in behavior while decreasing the overestimation caused by asking for the largest estimation. We compare our estimators with this method as well as the HPC classical model.

To conclude this section, we discuss the problematics of learning an unknown probability distribution from its samples. This problematic is one of the most natural and important questions in statistics and has been widely studied [START_REF] Braess | How to achieve minimax expected kullback-leibler distance from an unknown finite distribution[END_REF], [START_REF] Braess | Bernstein polynomials and learning theory[END_REF], [START_REF] Chan | Learning mixtures of structured distributions over discrete domains[END_REF], [START_REF] Kamath | On learning distributions from their samples[END_REF]. The typical strategy is to derive a good estimator (i.e. a function that transforms the samples) and to evalute it on a given loss function (hence the estimators are specific to the loss function). Kamath et al. [START_REF] Kamath | On learning distributions from their samples[END_REF] provide a good overview of the different work of the field. They describe how a theoretical evaluation should be performed. This is orthogonal to our work, and future studies could be dedicated to finding better and proved estimators for our loss function which would improve the performance of speculative scheduling.

III. MODEL AND ESTIMATORS

In this section we formalize the problem at hand. The core idea is the following: when we do not know the theoretical performance model of an application, can we use the sample of data available to find a good reservation strategy? Hence the underlying problem is the following: if we only have samples from the distribution of execution times of an application, can we estimate a reservation strategy that minimizes the expected cost for this distribution? We can also characterize how far are we using this strategy from the best one that has complete information about the application behavior.

A. Stochastic jobs and Reservation strategy

In this work, we use the model from our previous work [START_REF] Aupy | Reservation strategies for stochastic jobs[END_REF], [START_REF] Gainaru | Speculative scheduling for stochastic hpc applications[END_REF]. We consider stochastic jobs, that is jobs whose execution time X is unknown but (i) deterministic: two successive execution on the same input have the same duration; and (ii) randomly and uniformly sampled from a probability distribution law D (PDF: f X ). In this work, we consider that D is not known except for its domain [a, b] (essentially we use a = 0, and consider that the users can give an upper bound on the execution time). In addition, we have a sample of size n: X 1 , . . . , X n also iid (independent and identically distributed) from D.

A reservation strategy S = (t 1 , t 2 , . . . , t i , t i+1 , . . . ) for X is a set of increasing reservation sizes that are submitted to the machine in increasing order following the algorithm:

• A reservation of size t 1 is made. • If X ≤ t 1 ,
then the first reservation was enough for the successful execution of X and we are done. • If X > t 1 , then the reservation failed and the work is considered lost. A second reservation of size t 2 > t 1 is performed. • This process is repeated until the reservation is enough for the successful execution of X.

Finally, for a job of actual length t, we define the cost of a reservation t 1 as:

αt 1 + β min(t, t 1 ) + γ.
where the α part of this cost is the reservation cost, and the β part is the utilization cost. We have shown [START_REF] Aupy | Reservation strategies for stochastic jobs[END_REF] that this cost function could indeed cover many different cost models (from execution time on an HPC platform to cost of usage in the cloud). The studies in this work are done for α = 1 and β = γ = 0. For example, when simulating HPC systems, depending on how α is chosen, the cost will be determined by the total makespan (or response time) of an application by including the wait time in the queue until the scheduler submits the job for execution on the machine and/or the reservation time in the system.

Overall, the total cost for job t with strategy S is:

C(k, t) = k-1 i=1 (αt i + βt i + γ) + αt k + βt + γ
where k is the smallest index in the sequence such that t ≤ t k . The overall expected cost is then given by the equation:

Cost(S, X) = ∞ k=1 t k t k-1 C(k, t)f X (t)dt B.

Model, Optimization problem

In our previous work, we have shown that the general cost model could be written as: Theorem 1]). Given a random variable X and a strategy S = (t 1 , t 2 , . . . , t i , t i+1 , . . . ), we define the Expected Cost of S on X:

Theorem 1 (Expected Cost [2,
Cost(S, X) = β • E[X]+ ∞ i=0 (αt i+1 +βt i +γ) P(X > t i ) (1)
In this work, we use α = 1, β = 0, γ = 0 (RESERVA-TIONONLY model [START_REF] Aupy | Reservation strategies for stochastic jobs[END_REF]). In this case, the HPC system creates rigid reservations which are not adapted depending on how long the application is running (typical for several current large-scale systems). This scenario also assumes that large jobs have higher priority in being scheduled to be executed and thus do not wait longer in the queue.

Given a random variable X, we can compute S opt (X) the strategy that minimizes the Expected Cost on X. In statistics, estimators are often used to estimate an unknown parameter of a distribution. For example, if one knows that the distribution is a gaussian distribution, then given X 1 , • • • , X n the following are unbiased estimators

X = 1 n n i=1 X i for the mean (2)
S 2 n-1 = 1 n -1 n i=1 X i -X 2
for the variance

Their Mean-Square Error is O(1/n).

Using this we can now define a loss function for the estimators that we consider: Definition 1 (Expected Loss (EL)). Given a random variable X. We define the Expected Loss of a strategy S

d EL (S) = Cost(S, X) -Cost(S opt (X), X) Cost(S opt (X), X) (4)
Essentially, this loss function is the information about the proximity to the optimal algorithm. It is a distance to evaluate the quality of the estimators.

We can now define formally the problem that we consider here:

Definition 2 (Sample-based strategy). Given a set of n iid samples X 1 , • • • , X n for a random variable X following an unknown distribution. Find an estimator S n for S opt (X) that minimizes d EL .

IV. ESTIMATORS FOR RESERVATION STRATEGIES

In this section we provide different solutions to the Sample-Based strategy problem introduced in the previous section. We consider that we are given a sample data of size n: a) Discrete-Method: The most natural way to compute an estimate for the reservation strategy is to use the data as is: the n iid samples can be seen as a discrete distribution with n values, each of probability 1/n (assuming there is no redundant value, otherwise k/n if a value appears k times). We know that we can compute in polynomial time the optimal strategy for this distribution [2, Theorem 5]. We denote

X 1 , • • • , X n .
S n disc (X 1 , • • • , X n ) this strategy.
b) Continuous-Method: The second way is to transform the discreet data to a continuous domain and compute the optimal sequence of reservations in this domain. In order to do so, the idea is to transform the discrete cumulative function based on the n sample values into a continuous cumulative function (see the dashed line from Figure 3).

Given this transformation, we now have a cumulative function on a continuous domain. We then interpolate it with a continuous function using two different methods:

• Polynomial (P 10 -S n cont (X 1 , • • • , X n ))
: we find the polynomial of maximum degree 10 1 that minimizes the mean square error to the data. Figure 1 (lower) gives an example where a polynomial interpolation is a good fit.

• Usual-Distribution (UD-S n cont (X 1 , • • • , X n )):
we use the set of distribution available in scipy stats package (alpha, beta, cosine, dgamma, dweibull, exponnorm, exponweib, exponpow, genpareto, gamma, halfnorm, invgauss, invweibull, laplace, loggamma, lognorm, lomax, maxwell, norm, pareto, pearson3, rayleigh, rice, truncexpon, truncnorm, uniform, weibull) to find the best interpolation to the data with respect to the mean square error. Figure 1 (upper) gives an example of an application where a distribution interpolation gives the best fit.

For those two interpolation strategies, we can then compute reservation strategies using the different algorithms provided in [START_REF] Aupy | Reservation strategies for stochastic jobs[END_REF] (in this work, we use the discretization scheme: EQUAL-TIME and divide the support in 500 equal length steps).

V. SIMULATION ON SYNTHETIC APPLICATIONS

In this first set of simulations we evaluate the quality of our estimators on known distributions. This allows to measure exactly the value of d EL since knowing the shape of the data allows us to compute the optimal cost. We do not consider in this section actual traces of applications. We compare our estimators to estimators from the literature.

We are particularly interested whether our estimators are efficient estimators.

In the following we compare different estimators on several probability distributions:

1) Truncated gaussian (mean: 8h, domain: 0 to 20h) with low (2h) and wide (4h) variance; 2) An exponential distribution (mean: 8h); 3) A distribution whose density is the average of the density of two truncated gaussian (domain: 0 to 20h respectively, variance: 2) with different means (4h and 10h), representing an application with two different behaviors. We compare the estimators S n disc , [P 10 /UD ]-S n cont as introduced in Section IV to a semi-clairvoyant estimator which is expected to perform better and which can give us information on a good behavior for our estimators.

Specifically, the Semi-clairvoyant estimator is defined as a version of UD where the set of functions for the interpolation is restricted to the actual distribution of the data.

A. Speed of convergence

In the first set of evaluation, we are interested at the volume of samples needed to obtain a good solution. Hence we vary the number of data elements available (representing the number of previous runs) from 10 to 500 and measure the expected loss. Each evaluation is performed 100 times to evaluate the variation of our results (Figure 4).

To understand the results: an estimated loss of 0.05 with k sample points for estimator E essentially means that with a history of k runs, the strategy given by E gives a cost which is only 5% larger than that of the optimal strategy! Here we observe that if the distribution follows a truncated gaussian distribution with a history of 10 runs the cost of S n cont is only 6% larger than that of the optimal. With 60 runs we are down to 3%. It is interesting to observe the difference depending on the variance of the distribution: with low variance, the number of samples needed to be below 3% is small, but it is hard to get below 1%. With high variance, we need more runs to get below 3%, but we can get very close to optimal when the number of samples increases. The other observation is that S n disc has poorer performance than S n cont . S n disc over-fits the current data and is a inferior predictor of future data points. The results are similar for the exponential distribution (Figures 4 right).

For applications that shift their behavior due to configuration parameters (like in the case of profiles consisting of aggregated distributions), S n cont needs more samples to converge but it is much better than using the discrete data on it's own (Figure 5). These results confirm the intuition that going through a continuous interpolation of data makes the scheduling strategy more robust to variation.

From the results there are two key observations:

• The S n cont strategies are clearly dominating over the S n disc strategy; • The number of samples needed for good performance is relatively low.

B. Precision of polynomial approximation

In the second set of evaluation, we are interested in the impact of the maximum degree k of the polynomial in P k -S n cont on the performance of the fit. When the degree of the for the fitting function. While underfitting is expected to give poor performance, one can expect the same from overfitting the data: the ultimate overfitting of a function is the discrete data, which we showed was not as good as a continuous fit of the data to compute an efficient schedule. We write the result in Table I. Note that in this set of experiments, we have also added results from Multi-Atlas: traces from a real applications where we found that polynomial fit gave specifically poor results. In this case d EL is computed slightly differently (see Section VI).

From these results we make the following observations.

• Under-fitting is a problem for the performance of the algorithm, thus we need to fit polynomials of high enough order to capture the behavior of applications • It seems that over-fitting is less of a problem. However in some cases it shows poor results so it is a good idea to limit the order of the polynomial interpolation. Moreover, we observed that the fitting error of polynomial is not always correlated with the total cost. Overall, while the difference between a fitting with a maximum degree of 2 and 5 are notable, there is never a significant difference between results on a maximum degree of 5 and 15. This is true both for distributions where a polynomial fit gives poor results (Multi-Atlas), and distributions where is gives a good fit.

In future work, it would be interesting to show theoretically the observations here, or to find even better estimators.

VI. EVALUATIONS ON REAL APPLICATIONS

In this section we focus on real applications and study the applicability of our previous observations. For this pur-pose, we analyze the highly stochastic applications used by the The Medical-image Analysis and Statistical Interpretation laboratory at Vanderbilt. The dataset consists of over 30 neuroscience and medical applications for which we have execution logs of their runs on the ACCRE cluster (local to the Vanderbilt university). To compute the strategy, for each application we pick a number of samples (from 50 to 500) either consecutive or randomly chosen from history log and use the two strategies: S n cont and S n disc . We then evaluate their cost on the entire year of data by summing the cost of these strategies over all samples. Hence, the cost for application A is given by the following formula:

Cost(S, A) = t∈A C(k t , t) = N j=0 ( kj i=0 t i ),
where N is the total number of runs of application A and k is chosen so that t kj is the first reservation larger than the walltime of the j-th run (i.e. t kj > X k and t i < X k ∀i < k j ). These evaluations are done 100 time to study the variation of the performance.

Figure 6 shows the cost of each strategy relative to the optimal for 6 applications when varying the number of samples used for computing the sequence. The optimal strategy is computed by considering the discrete algorithm on the entire dataset.

As the number of samples increases, the S n disc strategy converges to S n cont . Note that in some cases, S n disc exceeds the performance of S n cont (e.g. Figure 6a). This is due to the fact that the number of data points available for the evaluation is actually very small, hence S n disc with many samples is almost the same as the computation of the optimal solution. However, in the majority of cases S n cont shows better results after only a few samples being used in the training phase (for example for the cortical surface modeling code, using only 10 runs for training gives a cost only 10% over the optimal). For applications with complex behaviors where the pattern of the application changes dramatically (for example Figure 6f shows execution times between 3 to 78 hours) both S n disc and S n cont have a high performance variation. We believe allowing polynomials of higher degree or more complex interpolations will be beneficial for these applications and continue to improve the results of S n cont . Further we analyzed the impact of how to chose the training samples. For all results we used consecutive runs in the logs in order to simulate a user having access to x previous runs and trying to predict the behavior of future ones. This might create a bias in the results since consecutive runs might share similar patterns. We made the same experiments again, this time randomly selecting the x sample runs (from 50 to 500) from the log. The results show that randomly choosing the samples can decrease the cost by 50% in some cases (a typical example presented in Figure 7). Depending on the number of samples in the execution log, it might be preferable to select a smaller random sample in order to compute the sequence for future runs. The final question that remains, is: in practice, how useful is it to be 20, 10, or even 2% far from the optimal strategy.

In this final evaluation section we try to answer this question by simulating a production system (like ACCRE at Vanderbilt) executing the 6 applications characterized in the previous section. For this purpose, we chose a week of execution from the history logs and extracting the application submission information. We used this data to feed it in an HPC scheduler simulator (ScheduleFlow [START_REF] Gainaru | ScheduleFlow: A simulator for HPC schedulers[END_REF] used in our previous studies) and compute the average job response time, system utilization and average job stretch (the ratio between the job response time and walltime) for different reservation strategies:

• HPC model in which successful completion is guaranteed by using the largest past execution as the estimate.

• The model used by the MASI group, by doing the Average of previous runs • Our Continuous fit model, by using random previous runs to compute the reservations using the S n cont strategy.

In some cases, the past largest execution is smallest than a future submission for the same application. In those cases, we extend all strategies by multiplying the last reservation by 1.5 until the reservation is large enough to complete the job.

Figure 8 presents the results. Using only 10 previous runs, the S n cont strategy decreases the average stretch by over 25%. Due to Multi Atlas runs that require a larger training sets, the average makespan and utilization do not see the same effect. However, including 110 samples into the training, the average job makespan decreases by 33% while the total system utilization increases by almost 50%.

The results show that even with little data one can use linear interpolation to obtain efficient scheduling strategies. While the performance of each individual application is different depending on the fitting functions and the number of samples used, randomly sampling previous runs seems to be a good optimization in all cases. In this conditions, speculative scheduling becomes feasible either integrated into a current HPC scheduler or as an additional layer at submission used to adjust request estimates. Our method can easily be adapted to the characteristics of any system by tunning the α, β and γ parameters. The goal of this work is to show that speculative scheduling can be used, even in the presence of incomplete information. We do not claim to have designed a good predictor for the performance of applications. Additional study on more applications need to be done. For instance, one can expect that the performance of some applications is also time dependent: updates to the software can be made for instance which intrinsically changes the code. Again, more research should be done in this regard.

We believe that we are opening a new direction in performance prediction research, where performance predictors do not need to necessarily predict the exact performance profile of an application, but can describe it with random variable and probability distributions. This may be an easier problem to solve for the field of performance prediction.

In addition, we have shown that even if the performance profile is not correct, we can still find benefits to it in speculative scheduling. Finally, an important takeaway from this work is the importance of the theoretical results for continuous distributions as well as those for discrete. This was not straightforward since historical data is more often given out as discrete samples. However, as we have seen in this work, interpolating the data makes for more robust solutions.

In the future, we plan to use this study as a basis to verify the upcoming theoretical results (e.g. using checkpoints so that the work is not lost; integrating multi-dimension such as memory needs, machine needs and performance models such as malleability).

The classes used to compute the sequence of requests are part of the ScheduleFlow simulator (available for download at: https://github.com/anagainaru/ScheduleFlow). After the successful execution of the script, a csv file will be created (wither dataset file cost.csv or distribution cost.csv. This file holds the results of one experiments per line where each line contains the fields: F u n c t i o n , F i t , P a r a m e t e r s , \ Cost , T r a i n s e t , EML Function is either Optimal, Discrete or Continuous. The fit column contains values only for the Continuous function and can either be a distribution or a order for the polynomial.

The parameters hold the parameters for the best interpolation function, wither the distribution parameters or the k-1 parameters for a k degree polynomial function. The Cost represents the cost for the sequence obtain using the continuous or interpolation fits, EML is the relative cost to the cost of the optimal. The trainset defines the number of samples used for computing the sequence (from 10 to 500).

E. Evaluation and expected result

The experiments ran 100 times each for a training set of 50 to 500 elements (in steps of 50).

f o r i i n r a n g e { 1 . . 1 0 0 } ; do \ p y t h o n c o m p u t e s e q u e n c e c o s t . py t r u n c n o r m ; \ done After running the experiments script, the results will be printed out in the corresponding csv file. This jupyter notebook print sequence cost.ipynb can be used to plot the figures from Section V.

We do not expect the results to be exactly the same as the ones presented in this paper but the observed trend should be similar.

F. Experiment customization

Fig. 1 .

 1 Fig. 1. Performance behavior for two stochastic functions together with two transformations into the continuous domain (top: Functional QA and preprocessing application, bottom: Deep brain structure segmentation application). Runs done during 2017 by the neuroscience department at the Vanderbilt University.

Fig. 3 .

 3 Fig. 3. From a discrete cumulative to a continuous cumulative function, for n = 20 data points, sampled from a truncated normal distribution (mean=8, variance=2, domain=12 to 14).

Fig. 4 .

 4 Fig. 4. Results for the truncated normal distributions (left: low variance, middle: high variance) and the truncated exponential distribution (right)

Fig. 5 .

 5 Fig. 5. Results for the sum of truncated normal distributions. For this one we do not plot the semi-clairvoyant algorithm since the cumulative function is not part of the available regular distributions.

  Fig. 6.

Fig. 7 .

 7 Fig. 7. Results for a structural identification of orbital anatomy application when using consecutive and random sampling

Fig. 8 .

 8 Fig. 8. Results for one week of execution (from top to bottom: system utilization; average job stretch and average job makespan)

Fig. 9 .

 9 Fig. 9. Behavior of the CDF and PDF for the rice and pearson3 distributions, which best fit the datasets presented.

4 )

 4 Datasets: All datasets are generated during runtime by the script (to follow either truncated normal or truncated exponential distributions).

  C. Installationg i t c l o n e \ g i t @ g i t h u b . com : a n a g a i n a r u \ / R e p r o d u c i b i l i t y I n i t i a t i v e . g i t cd R e p r o d u c i b i l i t y I n i t i a t i v e \ / 2 0 1 9 s c a l aD. Experiment workflowDownload or clone from git the source code.To execute the code using a log file of past executions for an application: p y t h o n c o m p u t e s e q u e n c e c o s t . py \ { f o l d e r / d a t a s e t f i l e } Dataset file needs to contain a list of execution times To execute the code using synthetic data for the execution times: p y t h o n c o m p u t e s e q u e n c e c o s t . py \ { d i s t r i b u t i o n } Supported distributions: [truncnorm, expon] The truncated normal distribution uses µ=8, σ=2 with limits [a, b]=[0,20]. The truncated exponential distribution uses µ=1, σ=1.5 with limits [a, b]=[0,9] These values can be changed inside the compute sequence cost.py script.

There are four possible customizations: 1 )

 1 The truncated normal and exponential parameters can be changed from the compute sequence cost.py script by changing the values from the ExponentialRun and TruncNormRun classes 2) The maximum degree of the polynomial interpolation can be changed in the compute sequence cost.py script when initializing W o r k l o a d F i t . P o l y I n t e r p o l a t i o n ( \ m a x o r d e r = 1 cost.py script by adding a class with the distribution parameter similar to ExponentialRun and TruncNormRun. The new class needs to inherit the DistributionRuns class. 4) The distributions used for the distribution interpolation can be changed from the DistInterpolation class in the WorkloadFit file.

TABLE I IMPACT

 I ON d EL OF THE MAX-DEGREE k OF THE POLYNOMIAL IN P k -S n CONT FOR DIFFERENT SAMPLE SIZES.

			Max Degree k
	Sample size = 10	2	5	10	15
	TNormal (µ = 8, σ = 2) 0.253	0.104	0.186	not fitted
	Exponential 0.203	0.091	0.092	0.241
	Sum of TNormal	0.198	0.024	0.021	0.021
	Multi Atlas 2.175	1.405	1.181	1.413
	Sample size = 60	2	5	10	15
	TNormal (µ = 8, σ = 2) 0.108	8.64e -3	5.52e -3	9.62e -3
	Exponential 0.182	2.41e -3	3.12e -3	9.39e -3
	Sum of TNormal	0.163	0.018	0.019	0.017
	Multi Atlas 1.274	1.036	0.811	1.13
	Sample size = 110	2	5	10	15
	TNormal (µ = 8, σ = 2) 0.121	9.61e -3	6.49e -3	9.82e -3
	Exponential 0.179	5.39e -3	8.12e -3	1.53e -2
	Sum of TNormal	0.211	0.031	0.034	0.032
	Multi Atlas 0.829	0.636	0.615	0.721
	Sample size = 510	2	5	10	15
	TNormal (µ = 8, σ = 2) 0.203	8.95e -2	0.103	9.22e -2
	Exponential 0.191	0.015	0.017	0.015
	Sum of TNormal	0.391	0.062	0.052	0.103
	Multi Atlas 1.061	1.26	0.955	0.821
	polynomial increases, more dimensions of freedom open up

Note that the number 10 for the maximum degree of the polynomial has been chosen to give some leniency while not overfitting the data.

https://github.com/anagainaru/ReproducibilityInitiative/tree/master/2019\ scala

ACKNOWLEDGEMENTS

We thank the VUIIS Center for Computational Imaging for sharing de-identified logs without patient or investigator identifiable data.

APPENDIX ARTIFACT DESCRIPTION: MAKING SPECULATIVE SCHEDULING ROBUST TO INCOMPLETE DATA A. Abstract

This artifact contains the code for computing the sequence of requests for an application given its past behavior (together with the corresponding cost as described in Section III), as well as instructions on how to use the code and print the data to reproduce the results from Section V. Section VI uses neuroscience applications that are open source, however the execution logs as well as the input MRI data for each application is not available outside the Vanderbilt University.

B. Description

1) Check-list (artifact meta information):

• Algorithm:

• Program: Python code • Data set: Synthetic application runs created from the truncated normal and truncated exponential distributions • Run-time environment: the code should execute successfully on any environment capable of running python, scipy and numpy • Hardware: any configuration • Execution: the compute sequence cost.py script can receive as input either a scipy distribution (between truncnorm and expon) or a file with a list of execution times (either one per line or separated by space) • Output: the output of the script is a csv file containing the optimal cost as well as (Best Fit, Parameters for the best fit, Cost, EML) for the discrete and continuous fit of the input data (when using 10 to 500 entries for training and all entries when computing the cost) • Experiment workflow: download the code, run compute sequence cost.py script on different distributions; obtain output csv files, run the jupyter notebook to print the data, analyze the output figure • Publicly available?: Yes 2) How software can be obtained (if available): The code used is available on github at https://github.com/anagainaru in the ReproducibilityInitiative repository (link in footnote) 2 .

3) Software dependencies: Python3 is required. The code is using scipy functions for fitting the data to a distribution. The pearson3 and rice have been removed from the list of distributions for which fitting is performed since often the CDF values given by scipy does not correspond with the PDF (Figure 9). Our results do not include these distributions. They can be added in the WorkloadFit.py file in the DistInterpolation class (in which case the outliers need to be removed from the results file).