Denys Dutykh
email: denys.dutykh@univ-savoie.fr

How to overcome the Courant-Friedrichs-Lewy condition of explicit discretizations?

Keywords: Heat equation, Finite differences, Explicit schemes, Implicit schemes, CFL condition MSC: [2010] 65M70, 65N35 (primary), 80M22, 76M22 (secondary) PACS: [2010] 47.11.Kb (primary), 44.35.+c (secondary) Heat equation, Finite differences, Explicit schemes, Implicit schemes, CFL condition

This manuscript contains some thoughts on the discretization of the classical heat equation. Namely, we discuss the advantages and disadvantages of explicit and implicit schemes. Then, we show how to overcome some disadvantages while preserving some advantages. However, since there is no free lunch, there is a price to pay for any improvement in the numerical scheme. This price will be thoroughly discussed below.

In particular, we like explicit discretizations for the ease of their implementation even for nonlinear problems. Unfortunately, when these schemes are applied to parabolic equations, severe stability limits appear for the time step magnitude making the explicit simulations prohibitively expensive. Implicit schemes remove the stability limit, but each time step requires now the solution of linear (at best) or even nonlinear systems of equations. However, there exists a number of tricks to overcome (or at least to relax) severe stability limitations of explicit schemes without going into the trouble of fully implicit ones. The purpose of this manuscript is just to inform the readers about these alternative techniques to extend the stability limits. It was not written for classical scientific publication purposes.

Introduction

In this text we consider the classical linear heat equation:

u t = ν ∇ 2 u , (1.1)
where u(x, t) is a quantity being diffused in some domain Ω ⊆ R d . In physical applications u(x, t) , x ∈ Ω , t > 0 may represent the temperature field, moisture content, vapor concentration, etc. and ν > 0 is the diffusion coefficient. The subscripts denote partial

derivatives, i.e. u t def := ∂u(x, t) ∂t . Finally, ∇ 2 ≡ ∇ • ∇ is the classical d-dimensional
Laplace operator:

∇ 2 def := d i = 1 ∂ 2 ∂x 2 i .
The derivation of this equation for Brownian motion process was given by A. Einstein [START_REF] Einstein | Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen[END_REF] in 1905.

From now on we shall restrict our ambitions on the 1-dimensional case where Ω ≡ [ 0, ℓ ] ⊆ R 1 and equation (1.1) correspondingly becomes:

u t = ν u x x . (1.2)
This equation has to be supplemented by one initial

u t = 0 = u 0 (x) ,
and two boundary conditions:

Φ l t, u(t, 0), u x (t, 0) = 0 , (1.3)
Φ r t, u(t, ℓ), u x (t, ℓ) = 0 . (1.4)
The functions Φ l, r ( • ) have to be specified depending on the practical situation in hands. For example, if we have the Dirichlet-type condition on the left boundary then

Φ l t, u(t, 0), u x (t, 0) ≡ u(t, 0) -u • (t) = 0 ,
where u • (t) is a prescribed function of time. Often, it is assumed that u • (t) ≡ const . The homogeneous Neumann-type condition on the right looks like

Φ r t, u(t, ℓ), u x (t, ℓ) ≡ u x (t, ℓ) = 0 .

Some healthy criticism

Obviously, the heat equation (1.1) is a simplified model obtained after a series of idealizations and simplifications. As a result, equation (1.1) is linear and its Green function * can be computed analytically

† G(x, t) = 1 √ 4 π ν t e - x 2 4 ν t
In particular, one can see that for any sufficiently small t > 0 the function G(x, t) is not of compact support. In other words, the information about a point source initially localized at x = 0 spreads instantly over the whole domain. Of course, the infinite speed of information propagation is physically forbidden. This non-physical feature of heat equation's solutions is a consequence of simplifying assumptions made during the derivation. We just mention here that some nonlinear versions of the heat equation do have fundamental solutions with compact support. However, in practice the solutions to Partial Differential Equations (PDEs) such as (1.1) are computed numerically than constructed analytically. That is why we can hope to correct some non-physical features of the heat equation solution at the discrete level. This is the main topic of the present manuscript.

Manuscript organization. Below we shall review some classical numerical schemes for the heat equation (1.1) in one spatial dimension in Section 2:

• The explicit scheme in Section 2.1 • The implicit scheme in Section 2.2 • The leap-frog scheme in Section 2.3 • The Crank-Nicholson scheme in Section 2.4

In that Section we justify our preference for explicit schemes in time. However, explicit schemes are known to have an important CFL-type stability restrictions on the time step [START_REF] Courant | Über die partiellen Differenzengleichungen der mathematischen Physik[END_REF]. That is why in Section 3 we review some not so widely known schemes which allow to overcome the stability limit while still being explicit in time. In particular, we consider the following alternatives:

• Dufort-Frankel method in Section 3.1 • Saulyev method in Section 3.2 • Hyperbolization method in Section 3.3 * The Green function is also known as the fundamental solution. More precisely, it solves the following problem:

G t = ν G x x , x ∈ R , G t = 0 = δ (x) ,
where δ(x) is Dirac delta function. † This expression is false. To be checked later! Finally, the main conclusions and perspectives of the present study are discussed in Section 4.

Classical numerical schemes

In order to describe numerical schemes in simple terms, consider a uniform discretization of the interval Ω Ω h :

Ω h = N -1 j = 0 [ x j , x j+1 ] , x j+1 -x j ≡ ∆x , ∀j ∈ 0, 1, . . . , N -1 .
The time layers are uniformly spaced as well t n = n ∆t , ∆t = const > 0 , n = 0, 1, 2, . . . The values of function u(x, t) in discrete nodes will be denoted by u n j def := u (x j , t n ) . The space-time grid is schematically depicted in Figure 1.

The Explicit scheme

The standard explicit scheme for the linear heat equation (1.2) can be written as The stencil of this scheme is depicted in Figure 2. In order to complete this discretization we have to find from boundary conditions (1.3), (1.4) the boundary values:

u n+1 j -u n j ∆t = ν u n j-1 -2 u n j + u n j+1 ∆x 2 , j = 1, . . . , N -1 , n 0 . (2.1)
u n+1 0 = ψ l (t n+1 , u n+1 1 , . . . ) , u n+1 N = ψ r (t n+1 , u n+1 N -1 , . . . ) , (2.2)
where functions ψ l, r ( • ) may depend on adjacent values of the solution whose number depends on the approximation order of the scheme (here we use the second order in space). For example, if the temperature is prescribed on the right boundary, then we simply have

u n+1 N = ψ r (t n+1 ) ≡ φ r (t n+1 ) ,
where φ r (t) is a given function of time. On the other hand, if the heat flux is prescribed on the left boundary ν ∂u ∂x = φ l (t) then it can be discretized as

ν -3 u n+1 0 + 4 u n+1 1 -u n+1 2 2 ∆x = φ l (t n+1 ) .
By solving (2.1) with respect to u n+1 j we obtain a discrete dynamical system

u n+1 j = u n j + ν ∆t ∆x 2 u n j-1 -2 u n j + u n j+1
, whose starting value is directly obtained from the initial condition:

u 0 j = u 0 (x j ) , j = 0, 1, . . . , N .
It is well-known that scheme (2.1) approximates the continuous operator to order O(∆t + ∆x 2 ) . The explicit scheme is conditionally stable under the following CFL-type condition:

∆t 1 2 ν ∆x 2 . (2.3)
Unfortunately, this condition is too restrictive for sufficiently fine discretizations.

The Implicit scheme

The implicit scheme for the 1D heat equation (1.2) is given by the following relations: The finite difference stencil of this scheme is depicted in Figure 3. These relations have to be properly initialized and supplemented with numerical boundary conditions (2.2). In the following Sections we shall not return to the question of initial and boundary conditions in order to focus on the discretization. The scheme (2.4) has the same order of accuracy as the explicit scheme (2.1), i.e. O(∆t + ∆x 2 ) . However, the implicit scheme (2.4) is unconditionally stable, which constitutes its major advantage. It could be interesting to have also the second order in time as well. This issue will be addressed in the following Sections.

u n+1 j -u n j ∆t = ν u n+1 j-1 -2 u n+1 j + u n+1 j+1 ∆x 2 , j = 1, . . . , N -1 , n 0 . (2.4)
The most important difference with the explicit scheme (2.1) is that we have to solve a tridiagonal system of linear algebraic equations to determine the numerical solution values u n+1 j N j = 0 on the following time layer t = t n+1 . It determines the algorithm complexity -a tridiagonal system of equations can be solved in O(N) operations (using the simple Thomas algorithm, for example) and it has to be done at every time step.

The Leap-frog scheme

The leap-frog scheme * is obtained by replacing in (2.1) the forward difference in time by the symmetric one, i.e.

u n+1 j -u n-1 j 2 ∆t = ν u n j-1 -2 u n j + u n j+1 ∆x 2 , j = 1, . . . , N -1 , n 0 . (2.5)
This scheme is second order accurate in space and in time, i.e. O(∆t 2 + ∆x 2 ) . Unfortunately, the leap-frog scheme is unconditionally unstable. It makes it un-exploitable in practice. However, we shall use some modifications of this scheme below.

The Crank-Nicholson scheme

We saw above that the first tentative to obtain a scheme with second order accuracy in space and in time was unsuccessful (see Section 2.3). However, a very useful method was proposed by Crank & Nicholson (CN) and it can be successfully applied to the heat equation (1.2) as well:

u n+1 j -u n j ∆t = ν u n j-1 -2 u n j + u n j+1 2 ∆x 2 + ν u n+1 j-1 -2 u n+1 j + u n+1 j+1 2 ∆x 2 , j = 1, . . . , N -1 , n 0 . (2.6)
This scheme is O(∆t 2 + ∆x 2 ) accurate and unconditionally stable (similarly to (2.4)).

That is why numerical results obtained with the CN scheme will be more accurate than implicit scheme (2.4) predictions. The stencil of this scheme is depicted in Figure 5. The CN scheme has all advantages and disadvantages (except for the order of accuracy in time) of the implicit scheme (2.4). At every time step one has to use a tridiagonal solver to invert the linear system of equations to determine solution value at the following time layer t = t n+1 .

Some nonlinear extensions

Since most of real-world heat conduction models used in building physics are nonlinear, it is worth to discuss some nonlinear extensions of the Crank-Nicholson (CN) scheme. For linear problems CN scheme turns out to be the same as the mid-point and trapezoidal rules for Ordinary Differential Equations (ODEs). Indeed, consider a nonlinear ODE:

u = f (u) , u(0) = u 0 . (2.7)
The mid-point and trapezoidal rules consist correspondingly in discretizing (2.7) as follows:

u n+1 -u n ∆t = f u n + u n+1 2 , u n+1 -u n ∆t = f u n + f u n+1 2 .
Now, if we set in formulas above f (u) = ν L • u , where L ≃ ∂ x x is the linear operator which represents the second central finite difference, we recover the CN scheme (2.6). Consider a non-conservative nonlinear heat equation:

u t = k(u) u x x . (2.8)
The straightforward application of the CN scheme to equation (1.1) yields the following scheme:

u n+1 j -u n j ∆t = k(u n j ) u n j-1 -2 u n j + u n j+1 2 ∆x 2 + k(u n+1 j ) u n+1 j-1 -2 u n+1 j + u n+1 j+1 2 ∆x 2 , j = 1, . . . , N -1 , n 0 . (2.9)
However, it is less known that one can apply also the cross-Crank-Nicholson (cCN) scheme:

u n+1 j -u n j ∆t = k(u n+1 j ) u n j-1 -2 u n j + u n j+1 2 ∆x 2 + k(u n j ) u n+1 j-1 -2 u n+1 j + u n+1 j+1 2 ∆x 2 , j = 1, . . . , N -1 , n 0 . (2.10)
We underline that both schemes (2.9) and (2.10) are second order accurate in space and in time, i.e. the consistency error is O(∆t 2 + ∆x 2 ) . However, there is a major advantage of the cCN scheme (2.10) over the classical CN scheme (2.9) in the fact that cCN is linear with respect to quantities evaluated at the upcoming time layer t = t n+1 provided that k(u) is an affine function of u . This fact can be used to simplify the resolution procedure without destroying the accuracy of the CN scheme. Otherwise, for more general diffusion coefficients k(u) the success of operation depends on the easiness to solve nonlinear equation (2.10). It goes without saying that information propagates instantaneously in both CN and cCN schemes.

Information propagation speed

Let us discuss now an important issue of the information propagation speed in the discretized heat equation (1.2). As the initial condition we take the following grid function:

u 0 j = 1 , j = 0 , 0 , j = 0 ,
which corresponds to the discrete Dirac function. In all fully implicit schemes (such as (2.4) and (2.6)) the grid function u 1 j N j = 0 will generally have non-zero values in all nodes (modulo perhaps homogeneous boundary conditions). Thus, we can conclude that information has spreaded instantaneously. On the other hand, as it is illustrated in Figure 1 with grey and red circles, in explicit discretizations the information propagates one cell to the left and one cell to the right in one time step. Thus, its speed c s can be estimated as

c s = ∆x ∆t CFL 2 ν ∆x .
Thus, the value of c s is finite. Of course, in the limit ∆x → 0 we recover the infinite propagation speed, but let us not forget that computations are always run for a finite value of ∆x . We arrived to an interesting conclusion. Even if the continuous heat equation (1.1) possesses an unphysical property, it can be corrected if we use a judicious (in this case explicit) discretization. This is the main reason why we privilege explicit schemes in time. However, these schemes are subject to severe stability restrictions. The rest of the manuscript is devoted to the question how to overcome the stability limit?

There is another computational advantage of explicit schemes over the implicit ones. Namely, explicit methods can be easily parallelized and they allow to achieve almost perfect scaling on HPC systems [START_REF] Chetverushkin | Explicit schemes and numerical simulation using ultrahigh-performance computer systems[END_REF]. Indeed, the computational domain can be split into subdomains, each sub-domain being handled by a separate processor. Since the stencil is local, only direct neighbours are involved in individual computations. The communication among various processes is almost minimal since only boundary nodes have to be shared. This is another good reason to privilege explicit schemes over the implicit ones.

Improved explicit schemes

Below we present several alternative methods which were specifically designed to overcome the stability limitation of the standard explicit scheme (2.1).

Saulyev method

In this Section we describe a not so widely known method proposed by Saulyev in [START_REF] Saulyev | Integration of parabolic equations by the grid method[END_REF] to integrate parabolic equations. For simplicity we focus on the 1-dimensional heat equation (1.2). The first idea of this method consists in rewriting the discrete second (spatial) derivative as

u x x x = x j ≈ u j+1 -2 u j + u j-1 ∆x 2 ≡ u j+1 -u j ∆x - u j -u j-1 ∆x ∆x .
In the finite difference formula above we do not specify intentionally the time layer number.

The next trick consists in writing the following asymmetric finite difference approximation:

u n+1 j -u n j ∆t = ν u n j+1 -u n j ∆x - u n+1 j -u n+1 j-1 ∆x ∆x ,
or after simplifications we obtain

u n+1 j -u n j ∆t = ν u n j+1 -u n j + u n+1 j + u n+1 j-1 ∆x 2 . (3.2)
The last difference relation is slightly different from the Dufort-Frankel method presented above. Moreover, the relation written above is not consistent with the original equation (1.2). That is why we consider the next time layer t = t n+2 and we apply symmetrically the same tricks, i.e.

u n+2 j -u n+1 j ∆t = ν u n+2 j+1 -u n+2 j ∆x - u n+1 j -u n+1 j-1 ∆x ∆x ,
correspondingly). However, this scheme can be recast in an almost explicit form using judicious recurrence relations. Consider the first stage (3.2) of Saulyev's scheme. Similarly to Section 3.1 we introduce for simplicity the parameter

λ def := ν ∆t ∆x 2 . From difference relation (3.2) we find u n+1 j = 1 -λ 1 + λ u n j + λ 1 + λ u n j+1 + λ 1 + λ u n+1 j-1 . (3.4)
The first stage of Saulyev's scheme is computed in rightwards direction (increasing index j ր). From the left boundary condition we find first the value

u n+1 0 = ψ l (t n+1 ) .
This allows us to compute u n+1 1 , u n+1 2 , etc. thanks to the recurrence relation (3.4). At the final step, the value u n+1 N is computed directly from the right boundary condition:

u n+1 N = ψ r (t n+1 ) .
This completes the first stage of computations.

Remark 1. For some types of boundary conditions ( e.g. Robin-type) Saulyev's scheme might require solution of a small dimensional (typically 2 × 2) linear system of algebraic equations to initiate the recurrence (3.4).

Let us make explicit now the second stage (3.3) of Saulyev's method. For this purpose we solve relation (3.3) with respect to u n+2 j :

u n+2 j = 1 -λ 1 + λ u n+1 j + λ 1 + λ u n+1 j-1 + λ 1 + λ u n+2 j+1 . (3.5)
Now it is getting clear that during the second stage of Saulyev's scheme we proceed in the leftwards direction (decreasing j ց). From the right boundary condition we find first

u n+2 N = ψ r (t n+2 ) .
It allows us to compute u n+2 N -1 , u n+2 N -2 , etc. using the recurrence relation (3.5). At the final step, the value u n+2 0 is computed directly from the left boundary condition:

u n+2 0 = ψ l (t n+2 ) .
As a result we obtain a fully explicit resolution scheme without stability related limitations. We notice however that Saulyev's scheme provides consistent results only every second time step or, in other words, after the successive completion of both stages (3.4) and (3.5). The intermediate result is not consistent with the equation (1.2).

Hyperbolization method

We saw above that the Dufort-Frankel scheme is a hidden way to add a small amount of 'hyperbolicity' into the model (1.2). In this Section we shall invert the order of operations: first, we perturb the equation (1.2) in an ad-hoc way and only after we discretize it with a suitable method.

Consider the 1-dimensional heat equation (1.2) that we are going to perturb by adding a small term containing the second derivative in time:

τ u t t + u t -ν u x x = 0 . (3.6)
This is the hyperbolic heat equation already familiar to us since it appeared in the consistency analysis of the Dufort-Frankel scheme. Here we perform a singular perturbation by assuming that τ u t t ≪ u t . The last condition physically means that the new term has only limited influence on the solution of equation (3.6). Here τ is a small ad-hoc parameter whose value is in general related to physical and discretization parameters τ = τ (ν, ∆x, ∆t) . Remark 2. One can notice that equation (3.6) is second order in time, thus, it requires two initial conditions to obtain a well-posed initial value problem. However, the parabolic equation (1.2) is only first order in time and it only requires the knowledge of the initial temperature field distribution. When we solve the hyperbolic equation (3.6), the missing initial condition is simply chosen to be u t t = 0 = 0 .

Dispersion relation analysis

The classical dispersion relation analysis looks at plane wave solutions: u(x, t) = u 0 e i (κ x -ω t) .

(3.7)

By substituting this solution ansatz into equation (1.2) we obtain the following relation between wave frequency ω and wavenumber k :

ω(κ) = -i ν κ 2 . (3.8)
The last relation is called the dispersion relation even if the heat equation (1.2) is not dispersive but dissipative. The real part of ω contains information about wave propagation properties (dispersive if Re ω(κ) κ = const and non-dispersive otherwise) while the imaginary part describes how different modes κ dissipate (if Im ω(κ) < 0) or grow (if Im ω(κ) > 0 ). The dispersion relation (3.8) gives the damping rate of different modes.

The same plane wave ansatz (3.7) can be substituted into the hyperbolic heat equation (3.6) as well to give the following implicit relation for the wave frequency ω :

-τ ω 2 -i ω + ν κ 2 = 0 .
Remark 3. The ad-hoc parameter τ can be chosen in other ways as well. One popular choice consists in taking

τ = ∆x c s ,
where c s is the real physical information speed.

Error estimate

It is legitimate to ask the question how far are solutions u h (x, t) to the hyperbolic equation (3.6) from the solutions u p (x, t) of the parabolic heat equation (1.1) (for the same initial condition). This question for the initial value problem was studied in [7] and we shall provide here only the obtained error estimate. Let us introduce the difference between two solutions: δu (x, t)

def := u h (x, t) -u p (x, t) .
Then, the following estimate holds

| δu (x, t) | τ M 1 + 2 √ π 8 √ 2 τ + 4 √ 2 π 2 2 T ,
where T > 0 is the time horizon and

M def := sup Ω ξ, ζ ∂ 2 u p ∂ t 2 (ξ, ζ) ,
and the domain Ω ξ, ζ is defined as

Ω ξ, ζ def := (ξ, ζ) : 0 ζ t , x - t -ζ √ τ ξ x + t -ζ √ τ .

Discussion

We saw above that only explicit schemes allow to have the finite speed of information propagation in the discretized version of the heat equation (1.1). The ease of parallelization along with excellent scaling properties of the codes obtained with explicit schemes constitute another important advantage to privilege explicit schemes over implicit ones [START_REF] Chetverushkin | Explicit schemes and numerical simulation using ultrahigh-performance computer systems[END_REF]. However, explicit schemes for parabolic equations suffer from very stringent CFL-type conditions ∆t = O(∆x 2 ) on the time step. That is why it is almost impossible to perform long time simulations (required in e.g. building physics applications) using fully explicit schemes such as (2.1). In order to keep explicit discretizations and overcome stringent CFL-type restrictions, a certain number of hybrid schemes were described. These hybrid schemes are based on different ideas. Some schemes rely on the information about the numerical solution on following time layers while keeping the overall scheme explicit using various tricks. Some hybrid schemes (e.g. Dufort-Frankel method described in Section 3.1 and Saulyev's scheme from Section 3.2) can be even unconditionally stable.

Using the local error analysis we noticed that CFL-improved schemes gain the stability by introducing some weak hyperbolicity into the model (1.2). In particular, it is the case of the Dufort-Frankel scheme, while Saulyev method seems to be rather dispersive. This observation suggests that a new method can be derived by introducing this hyperbolicity in a controlled manner and to discretize the perturbed equation later. The method of hyperbolization deforms the equation operator to achieve desired properties * of the numerical solution. However, in all cases the gain in stability results in some loss of accuracy in representing the original continuous operator (1.1). The trade-off between the stability and accuracy has to be made by the end user. The second order accuracy is equivalent to the classical parabolic CFL-type condition. However, the user can choose to degrade intentionally the accuracy to relax the stability restriction up to hyperbolic-type conditions ∆t = O(∆x) (originally found in [START_REF] Courant | Über die partiellen Differenzengleichungen der mathematischen Physik[END_REF]) and even beyond. It is not difficult to see that hybrid schemes presented in this study can be easily generalized to nonlinear cases with source terms in one and more spatial dimensions.

The main goal of this manuscript was to communicate and attract community's attention to these improved discretizations, which can be used in modern building physics simulations where typical time scales are measured rather in months or even years. Our preference goes perhaps to the method of hyperbolization since it has been successfully applied (and validated) even to compressible Navier-Stokes and MHD equations [START_REF] Chetverushkin | Hyperbolic type explicit kinetic schemes of magneto gas dynamic for high performance computing systems[END_REF]. Another advantage of hyperbolization technique is that it can be mathematically derived for gas dynamics from the kinetic theory of Boltzmann.

While preparing this text, the Author discovered also another strategy to relax the the time step limits. This idea can be traced back to the scientific school of R. Temam who proposed to separate the scales for the purposes of numerical simulations. Then, high frequencies, which impose severe stability limits, are treated separately. For a modern introduction to these approaches we refer to [START_REF] Brachet | Stabilized Times Schemes for High Accurate Finite Differences Solutions of Nonlinear Parabolic Equations[END_REF]. There exist also explicit time integration methodsThe so-called Runge-Kutta-Tchebyshev methods. especially designed to have extended stability limits. This research direction was pioneered also by V. K. Saulyev in 1960 [START_REF] Saulyev | Integration of parabolic equations by the grid method[END_REF].

Contents 1

 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Some healthy criticism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Classical numerical schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 The Explicit scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 The Implicit scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 The Leap-frog scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 The Crank-Nicholson scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Some nonlinear extensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Information propagation speed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Improved explicit schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Dufort-Frankel method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Saulyev method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Resolution procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Hyperbolization method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Dispersion relation analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Error estimate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 .

 1 Figure 1. A schematic representation of the uniform discretization in space and time. Red nodes correspond to non-zero values of the discrete solution. Grey nodes correspond to u n j ≡ 0 . The shaded area is an equivalent of the ' light cone' for the initially activated node.

Figure 2 .

 2 Figure 2. Stencil of the explicit finite difference scheme (2.1).

Figure 3 .

 3 Figure 3. Stencil of the implicit finite difference scheme (2.4).

Figure 4 .

 4 Figure 4. Stencil of the leap-frog (2.5) and hyperbolic (3.9) finite difference schemes.

Figure 5 .

 5 Figure 5. Stencil of the Crank-Nicholson (CN) finite difference scheme (2.6).

Figure 6 .

 6 Figure 6. Stencil of the Dufort-Frankel (3.1) finite difference scheme.

* This scheme is called in French as 'le schéma saute-mouton'.

Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Acknowledgments

The author would like to thank Professor B.N. Chetverushkin whose talks available on-line greatly inspired the preparation of this manuscript. I would like also to thank my collaborators J. Berger, S. Gasparin and N. Mendes for bringing my attention to heat and mass conduction problems.

Dufort-Frankel method

Let us take the unconditionally unstable leap-frog scheme (2.5) and slightly modify it to obtain the so-called Dufort-Frankel method:

) where we made a replacement

The scheme (3.1) has the stencil depicted in Figure 6. At the first glance the scheme (3.1) looks like an implicit scheme, however, it is not truly the case. Equation (3.1) can be easily solved for u n+1 j to give the following discrete dynamical system:

where

The standard von Neumann stability analysis shows that the Dufort-Frankel scheme is unconditionally stable.

The consistency error analysis of the scheme (3.1) shows the following interesting result:

So, from the asymptotic expansion for L n j we obtain that the Dufort-Frankel scheme is second order accurate in time and

However, the Dufort-Frankel scheme is unconditionally consistent with the so-called hyperbolic heat conduction equation:

We shall return to this equation below. At this stage we only mention that information propagates with the finite speed in hyperbolic models. and after simplifications we have

Without any surprise the relation (3.3) does not approximate equation (1.2) either. However, both relations (3.2) and (3.3) constitute the so-called Saulyev scheme and are called the first and second stages of Saulyev method correspondingly.

In order to perform the approximation error analysis we take the sum of (3.2), (3.3) and we divide it by two:

After applying local Taylor expansions we obtain

From the last asymptotic expansion we arrive at an important result: Saulyev scheme is second order accurate in space if ∆t = O(∆x 3/2 ) . We underline the fact that this condition coming from the accuracy requirements is weaker than the usual CFL restriction (2.3). For instance, if the user is ready to sacrifice the spatial accuracy to the first order O(∆x) , then it is sufficient to take ∆t = O(∆x) . Without proof we report that Saulyev's scheme is unconditionally stable (see [START_REF] Saulyev | Integration of parabolic equations by the grid method[END_REF] for more details). The stencil of Saulyev's scheme is depicted in Figure 7.

Resolution procedure

At the first glance Saulyev scheme appears as an implicit one since each relation (3.2) and (3.3) contains two terms from the following time layer (t = t n+1 and t = t n+2 By solving this quadratic equation with complex coefficients for ω , we obtain two branches:

This dispersion relation will be analyzed asymptotically with τ ≪ 1 being the small parameter. The branch ω -(κ) is not of much interest to us since it is constantly damped, i.e.

It is much more instructive to look at the positive branch ω + (κ) :

The last asymptotic expansion shows that for small values of parameter τ we obtain a valid asymptotic approximation of the dispersion relation (3.8) for the heat equation (1.2).

Discretization

Equation (3.6) will be discretized on the same stencil as the leap-frog scheme (2.5) (see Figure 4):

The last scheme is consistent with hyperbolic heat equation (3.6) to the second order in space and in time O(∆t 2 + ∆x 2 ) . Indeed, using the standard Taylor expansions we obtain

The stability of the scheme (3.9) was studied in [START_REF] Chetverushkin | Explicit schemes and numerical simulation using ultrahigh-performance computer systems[END_REF] and the following stability condition was obtained: ∆t ∆x τ ν . By taking, for example, τ = ν ∆x we obtain the following stability condition

which is still weaker than the standard parabolic condition (2.3). However, it was reported in [START_REF] Chetverushkin | Hyperbolic type explicit kinetic schemes of magneto gas dynamic for high performance computing systems[END_REF][START_REF] Chetverushkin | An explicit scheme for the solution of the filtration problems[END_REF] that stable computations (even in 3D) can be performed even with ∆t = O(∆x) . The authors explain informally this experimental observation by the fact that usual stability conditions are too 'pessimistic'.


Antoine Toueir
email: toueir@irit.fr

Julien Broisin
email: broisin@irit.fr

Michelle Sibilla
email: sibilla@irit.fr

Goal-Oriented Monitoring Adaptation: Methodology and Patterns

Keywords: Quality requirements, adaptive monitoring, autonomic systems, goal-oriented adaptation

This paper argues that autonomic systems need to make their distributed monitoring adaptive in order to improve their "comprehensive" resulting quality; that means both the Quality of Service (QoS), and the Quality of Information (QoI). Thus, we propose a methodology to design monitoring adaptation based on high level objectives (goals) related to the management of quality requirements. One of the advantages of adopting a methodological approach, is that monitoring reconfiguration will be conducted through a consistent adaptation logic. Starting from a model-guided monitoring framework, we introduce our methodology to assist human administrators in eliciting the appropriate quality goals piloting the monitoring. Moreover, some monitoring adaptation patterns falling into reconfiguration dimensions are suggested and exploited in a cloud provider case-study illustrating the adaptation of Quality-Oriented monitoring.

Introduction

Autonomic systems that are implemented by virtue of their four characteristics self-configuration, self-optimization, self-healing and self-protection, are serving the main principle of making them self-managed to achieve high level objectives [START_REF] Kephart | The vision of autonomic computing[END_REF]. In practice, the four self-* characteristics are realized by implementing the MAPE-K (Monitoring, Analyzing, Planning, Executing -Knowledge) loop modules. This implementation is either embedded within a resource, or distributed over several resources. However, the monitoring module of MAPE-K loop plays a crucial role, since wrong decisions might be taken by the analyzing & planning modules, if they were provided with interrupted or wrong information. Therefore, autonomic systems need to ensure the quality of information (e.g., correctness, freshness, timeliness, accuracy, etc.) exposed by the distributed monitoring modules.

Within autonomic systems, monitoring is usually quality-oriented. In other words, the underlying monitoring instruments metrics and evaluates them against quality specifications expressed via Service Level Agreements (SLAs) or management high level objectives. Since the management system could provide the possibility to renegotiate or modify the QoS specification afterward, and also, various management needs could be distinguished during the management system lifetime, the monitoring system has to adapt its behavior according to these new requirements. To resume, the monitoring of autonomic systems has to be capable of configuring the underlying gathering mechanisms (i.e., polling & notification) carrying the monitoring functions (e.g., measuring, gathering, evaluating, filtering, etc.) starting from quality specification, as well as reconfiguring those mechanisms based on quality requirements.

Most of the time, reconfiguration is held through ad-hoc logic. But this approach isn't suitable for reuse in other scenarios, and also doesn't satisfy high level objectives; which are extended at the autonomic system whole scale. To overcome these issues, first, we adopted the Requirements Engineering methodology to design monitoring adaptation; it starts from high level goals, and ends up with the (re)configuration of monitoring mechanisms. However, reconfiguration questions such as: why to delay launching some monitoring mechanisms? Why to substitute remote agents? What determine the instrumentation of particular set of metrics but not another one? need to be answered. In other words, identifying goals representing the "starting point" for deriving monitoring (re)configuration is a big challenge. Thus, besides adaptation methodology, this paper answers these questions by proposing monitoring adaptation patterns to assist human administrators in designing meaningful adaptations, that increase the overall quality of autonomic systems.

The paper is organized as follows: Section 2 points out the weaknesses of other monitoring adaptation approaches; Section 3 gives an overview of the monitoring framework at the basis of our contributions; the methodology of designing goal-oriented adaptive monitoring is presented in Section 4; the reconfiguration dimensions as well as monitoring adaptation patterns are discussed in Section 5 and then applied on a case-study in Section 6; finally, we conclude by listing our perspectives in Section 7.

Related Work

This section enumerates some existing trends focusing on (i) adapting the QoS monitoring in autonomic systems [START_REF] Grefen | Crossflow: Cross-organizational workflow management in dynamic virtual enterprises[END_REF][START_REF] Roxburgh | Monitoring as an sla-oriented consumable service for saas assurance: A prototype[END_REF][START_REF] Thongtra | An adaptable capability monitoring system[END_REF][START_REF] Munawar | Adaptive monitoring with dynamic differential tracing-based diagnosis[END_REF][START_REF] Nobre | Decentralized detection of sla violations using p2p technology[END_REF], and (ii) designing patterns regarding the distributed deployment as well as the adaptation of MAPE loop modules [START_REF]An architectural blueprint for autonomic computing[END_REF][START_REF] Weyns | On patterns for decentralized control in self-adaptive systems[END_REF][START_REF] Ramirez | Design patterns for developing dynamically adaptive systems[END_REF].

Collecting additional metrics or joining managed resources is addressed in [START_REF] Grefen | Crossflow: Cross-organizational workflow management in dynamic virtual enterprises[END_REF][START_REF] Roxburgh | Monitoring as an sla-oriented consumable service for saas assurance: A prototype[END_REF][START_REF] Munawar | Adaptive monitoring with dynamic differential tracing-based diagnosis[END_REF] either to adapt monitoring to meet SLA modifications, or to deal with the managed scope changes, or even to operate a "minimal" monitoring that is able to be extended in case of SLA violations. Indeed, the capability of scaling up/down the monitored metrics and resources is important as an adaptation action. But, it isn't clear whether this capability could be applied in other scenarios for distinct objectives, if so, how that could be feasible.

Runtime deployment of monitoring resources (i.e., managers, probes) is discussed in [START_REF] Roxburgh | Monitoring as an sla-oriented consumable service for saas assurance: A prototype[END_REF][START_REF] Thongtra | An adaptable capability monitoring system[END_REF][START_REF] Nobre | Decentralized detection of sla violations using p2p technology[END_REF] either to integrate monitoring into the SLA management life-cycle of large scale systems, or to replace failed managers, or even to monitor metrics concerning particular paths or segments. But here also, besides the undeniable gains of deploying monitoring resources during runtime, we don't see how the system administrators can orchestrate the monitoring adaptation of the distributed modules among several collaborating managers. The orchestration isn't a trivial task, because a given quality objective may need to be extended on several managers and treated differently on each of them. Thus, a "simple deployment" of new monitoring resources isn't enough to realize orchestration.

Inspired from the autonomic computing reference architecture proposed in [START_REF]An architectural blueprint for autonomic computing[END_REF], patterns regarding the distribution of the MAPE loop modules were proposed in [START_REF] Weyns | On patterns for decentralized control in self-adaptive systems[END_REF][START_REF] Ramirez | Design patterns for developing dynamically adaptive systems[END_REF]. Those patterns are useful in term of design reuse as well as clarifying the application contexts and benefits, but they target mainly the deployment of the monitoring modules rather than the monitoring behavior itself. That is, once the autonomic system is designed and deployed based on the proposed patterns, the monitoring conserves its behavior. Consequently, the management system knowledge won't exceed a "maximum ceiling" and the management will be limited regarding treating new situations; as it has specific vision reflecting the same aspects all the time.

To resume, most of the studied work focuses either on the auto-configuration of monitoring, or on the reconfiguration of the functional system based on the knowledge produced from the QoS monitoring, or even the reconfiguration of the monitoring itself based on static purposes. Contrary to the studied work, we argue that increasing the QoS of the whole system begins from designing adaptive monitoring, that is derived from and satisfies high level quality requirements, and interpreted by several managers.

The Enriched Adaptive Monitoring Framework

Our approach is based on a 3-layered framework [START_REF] Moui | A cim-based framework to manage monitoring adaptability[END_REF][START_REF] Moui | Information models for managing monitoring adaptation enforcement[END_REF][START_REF] Moui | Managing polling adaptability in a cim/wbem infrastructure[END_REF] illustrated in Figure 1, and defines three capabilities required to control monitoring: being configurable, adaptable and governable. This framework operates monitoring mechanisms without any consideration regarding agents or management protocols.

The configurability layer stands on the Distributed Management Task Force (DMTF) Common Information Model (CIM) standard. This low level layer aims at representing, in addition to the managed resources, the metrics [START_REF] Toueir | Toward configurable performance monitoring: Introduction to mathematical support for metric representation and instrumentation of the cim metric model[END_REF] and their gathering mechanisms as well [START_REF] Moui | A cim-based framework to manage monitoring adaptability[END_REF]. Moreover, those models have been enriched with the concept of Monitoring Mode. The latter encapsulates: metrics (Aspects of Interest to be instrumented), constraints (thresholds to be checked once associated metrics are instrumented), indications (notifications to be raised once metrics are violated), as well as two types of subscriptions to deliver both metrics values and raised indications to the appropriate destinations (see Figure 1). The adaptability layer provides an interface encapsulating operations to be applied on the lower layer models. Those operations constitute a "control interface" to update the attributes of both monitoring modes and metrics, thus the underlying gathering mechanisms will be reconfigured, and consequently the monitoring adapts its behavior. Finally, the governability layer is the top level layer representing the "intelligence" of the monitoring adaptation. To express the quality requirements, it uses Event/Condition/Action (ECA) policies to describe when and how adaptation should take place, that is, they determine the situations during which the adaptability layer operations should be invoked. Thinking that one of the existing software engineering approaches should answer our needs, we looked forward a suitable method for designing monitoring adaptation. The origin of Requirements Engineering (RE) goes back to the need to avoid crucial mistakes committed at the project design phase, and aims at building systems solving real-world problems. This methodology applies iterative activities about "eliciting, evaluating, documenting, consolidating and changing the objectives, functionalities, assumptions qualities and constraints that the system-to-be should meet based on the opportunities and capabilities provided by new technologies" [START_REF] Van Lamsweerde | Requirements Engineering: From System Goals to UML Models to Software Specifications[END_REF]. Among multiple RE approaches, Keep All Objectives Satisfied (KAOS) is adopted as RE goal-oriented method, due to its formal assertion layer that proves correctness and completeness of goals [START_REF] Van Lamsweerde | Requirements engineering in the year 00: a research perspective[END_REF]. In KAOS the system-to-be is divided into various models. The goal models elicit goals representing high level behavioral prescriptions of the system-to-be (the monitoring system in this study). Each goal may serve one or more objectives, and a given goal is realized through the cooperation of several components or actors, the so-called Agents1 . Goals are decomposed into sub-goals via a refinement process (see Section 6), where the most refined goals are called Requirements, or Leaf Goals; in KAOS, each of those leaf goals is assigned to a specific agent in order to be realized. Others goals, depicted within the obstacle models, are deduced from the goal models and prevent the satisfaction of the latter. In our study, since goals converge on the quality of monitoring, their obstacles will be related to the monitoring failures. The operation models are composed of the sets of internal operations carried by agents in order to realize the Leaf Goals. Finally, the object models identify the system-to-be objects.

By iterating refinement process on goals and obstacles, leaf goals will be identified. Once leaf goals are determined, both of ECA policies (to be inserted into the governability layer) and agents (invoking operations of the adaptability layer) will be reconnoitered. Thus, monitoring adaptation is automatically handled, and high level goals remain satisfied. However, human administrators have to refine manually the high level objectives they want to reach, in order to identify the leaf goals. To facilitate this task, we investigated the monitoring aspects that are subject to adaptation. As a result, we have identified various leaf goals belonging to four dimensions (i.e., Spatial, Metric, Temporal, Exchange) [START_REF] Toueir | A goal-oriented approach for adaptive sla monitoring: a cloud provider case study[END_REF]. In Section 5, we pursuit in proposing monitoring adaptation patterns falling into those dimensions, in order to assist human administrators in refining goals.

Dimensions & Patterns

With regard to refinement process, besides the basic AND/OR-decompositions, we rely on some predetermined correct and complete refinement patterns proved mathematically [START_REF] Darimont | Formal refinement patterns for goal-driven requirements elaboration[END_REF]. Those patterns refine Achieve goals of the form P ⇒ ♦Q (see Table 1), and written in Linear Temporal Logic (LTL) classical operators where ♦, and W mean respectively some time in the future, always in the future, and always in the future unless. Starting from a given goal (P ), milestone pattern identifies one or many intermediate goals (R, [...]) that must be reached orderly before reaching the ultimate one (Q). Rather, case pattern identifies the set of different and complete cases (P 1, P 2) for reaching final goals (Q1, Q2) that OR-decompose the ultimate goal (Q). Finally, the guard pattern requires the recognition of a condition (R) before achieving the ultimate goal (Q).

In order to clarify the exploitation contexts, pattern goals and requirements, as well as the various situations in which they may apply to, our pattern structure encompasses: context, pattern refinement, and examples. Notice that we are focusing on adaptation actions taken at the autonomic manager side only. Thus, investigating adaptations at the agent side is out of scope. In addition, the patterns are refined using KAOS graphical language [START_REF] Van Lamsweerde | Requirements Engineering: From System Goals to UML Models to Software Specifications[END_REF].

P ⇒ ♦R R ⇒ ♦Q Case Pattern P ∧ P 1 ⇒ ♦Q1 P ∧ P 2 ⇒ ♦Q2 (P 1 ∨ P 2) Q1 ∨ Q2 ⇒ Q Guard Pattern P ∧ ¬R ⇒ ♦R P ∧ R ⇒ ♦Q P ⇒ P W Q

Exchange Dimension Pattern

Context. Relying on IBM blueprint reference architecture [START_REF]An architectural blueprint for autonomic computing[END_REF], autonomic systems could distribute the MAPE loop over multiple collaborating autonomic managers. Each of which is responsible for managing particular scope of managed resources. Patterns belonging to this dimension are useful to overcome metrics gathering and delivering problems. Those problems could affect either metrics values, communication reliability between the information sources and destinations, or even on trustworthiness of those sources and destinations.

Pattern Refinement. Communications inside autonomic system could be classified according to the entities involved in information exchange (i.e., managers, agents, shared databases). Therefore, we identify three communication classes: Manager-2-Agent, Manager-2-Manager, and Manager-2-Shared Database (see Figure 1). By taking into consideration push and pull modes 2a). Since a manager responds to pull requests in both Manager-2-(Manager/ Shared DB) pull mode communications, it is considered as agent (because it is the source of information); therefore, this case becomes identical to Manager-2-Agent pull mode. Moreover, adaptation actions related to Manager-2-Agent push mode are not treated because they need to be held at the agent side.

Examples. This pattern is suitable for the following cases: (1) Increasing accuracy or precision of pulled/pushed metrics values, by replacing information source. (2) Querying more available agents, or blocking fake agents trying to integrate the distributed management system. (3) Securing the communication between information sources and destinations. (4) Modifying information destination when changing the topology of collaborating autonomic managers.

Metric Dimension Pattern

Context. The main idea behind building autonomic systems is to delegate decisions, that human administrators are used to make, to the autonomic systems themselves. Thus, to be able to make "wise" decisions, monitoring system needs to instrument specific metrics that could be activated/deactivated according to the management needs during runtime. Patterns belonging to this dimension are useful to control the trade-off between constructing richer knowledge and monitoring the information that is necessary for management.

Pattern Refinement. Metric instrumentation must be thought at the whole management system level. In other words, a given autonomic manager could activate/deactivate instrumentation of particular metrics, but when deactivating metrics on that manager, it doesn't mean necessarily that those metrics are "abandoned", because they could be transferred to other collaborating autonomic manager on which they are activated. These two cases are OR-decomposing the first refinement level. Regarding metrics manipulation inside autonomic manager, the second refinement level uses case pattern to cover metric classes. Our research exploits both CIM Metric Model [START_REF]Base metrics profile[END_REF] classifying metrics into Base, Discrete & Aggregation, as well as our mathematical extension [START_REF] Toueir | Toward configurable performance monitoring: Introduction to mathematical support for metric representation and instrumentation of the cim metric model[END_REF] classifying base metrics into Resource, Measurable & Mathematical. Each of these classes is OR-decomposed using Add Aspects and Remove Aspects leaf goals. On the other hand, the transfer of metrics among autonomic managers could be refined through milestone pattern, when metrics are activated on the collaborating manager first (Add Aspects in Figure 2b, as Subgoal 1 in Table 1), and then removed from the delegating one (Remove Aspects, as Subgoal 2).

It is worthy to precise that previously mentioned Aspects are representing "metric definitions", rather than "metric values". The former encompasses attributes related to the nature of metric (e.g., data type, unit), while the latter describes the instrumented values and their relevant contexts.

Examples. This pattern can be applied in the following cases: (1) Troubleshooting, or applying root cause analysis algorithms. ( 2) Modifying the hierarchical topology of the management system by instrumenting aggregated metrics to be exported to other managers or shared DBs. ( 3) "Engineering" the distribution of monitored metrics among autonomic managers.

Spatial Dimension Pattern

Context. As mentioned earlier, in autonomic system, each manager is responsible for managing a set of managed resources. In many cases, the number of users consuming the autonomic system services may oscillate rapidly, or even become quite important in term of size. Thus, managed resources are subject to be joined/withdrawn during runtime. Patterns belonging to this dimension are useful to react to the important changes of the managed resources scope.

Pattern Refinement. As management of autonomic systems is orchestrated by the collaboration of multiple autonomic managers, each of which can act on its own perimeter, as well as the perimeters of its collaborating peers. Thus, the first refinement level uses case pattern to cover these two cases. In fact, acting on its own perimeter is OR-decomposed using Expand and Shrink Monitoring Perimeter leaf goals to join/withdraw resources respectively into its managed scope. Rather, acting on others perimeters is refined using case pattern into deploying a new manager, or soliciting an existing one. First, the case of deploying a new manager is refined using milestone pattern into launching manager (Launch Delegated Manager in Figure 2c, as Subgoal 1 in Table 1), and then, delegating perimeter (Delegation, as Subgoal 2). In its turn, delegation goal is also refined through milestone pattern into joining delegated perimeter on the delegated manager (Expand Perimeter, as Subgoal 1), and then, deleting this perimeter from the delegating manager (Shrink Perimeter, as Subgoal 2). About the second case, where acting is held on existing manager, it is refined twice, using milestone pattern, into delegating the whole perimeter to the delegated manager (Delegation, as Subgoal 1), and then shutting down the delegating one (Shutdown Delegating Manager, as Subgoal 2).

Examples. This pattern is suitable for the following cases: (1) Load balancing of monitoring among autonomic managers. (2) Supporting scalability of the autonomic systems. (3) Decreasing the number of monitoring resources.

Temporal Dimension Pattern

Context. Temporal aspects are decisive factors in adapting monitoring behavior. Notice that the previous patterns are explained without considering time notions, but in fact, they imply some temporal aspects. Patterns belonging to this dimension are useful to overcome, among others, both temporal violations and scheduling problems, as well as to tune the analysis over the instrumented metrics. However, these two cases are far from being exhaustive, and time intervenes in a lot of other cases.

Pattern Refinement. Regarding information exchange, once again, we use case pattern to represent the same cases identified in exchange dimension. Obviously, dealing with information exchange temporal aspects, means that exchange is done iteratively and not once. Thus, Manager-2-Agent case is OR-decomposed into periodic poll, and both Manager-2-(Manager/Shared DB) cases are ORdecomposed into periodic export. We distinguish two levels of temporal granularity: the fine-grained level deals with an individual polling (exporting), whereas the coarse-grained level addresses a collective polling (exporting). Based on this distinction, we identify various leaf goals OR-decomposing periodic poll (export), where: Update Polling (Exporting) Period to update the frequency of a given polling (exporting), Align Polling (Exporting) to launch a set of synchronized parallel pollings (exported metrics) at the same time, and Misalign Polling (Exporting) to launch pollings (exported metrics) according to a "relative offset" delaying their launching moments one another (see Figure 2d).

Regarding metrics calculation, we identify the case of modifying the temporal interval covered by the metric value3 . Therefore, case pattern is used twice to cover all possible metric classes previously mentioned. But, we refine only measurable, mathematical & aggregation metrics, because time has a sense in their calculation, but not the others. Thus, we OR-decompose measurable & mathematical metrics using Update Time Scope Interval, while Update Time Series Interval OR-decomposes aggregation metrics (see Figure 2d).

Examples. This pattern is suitable for the following cases: (1) Controlling (e.g., relaxing, stressing) the monitoring load on autonomic managers, network paths among autonomic managers and shared DBs, as well as remote agents. (2) Tuning temporal parameters of metrics analysis.

Case-Study

Context. Our scenario rolls in a cloud data center providing to the virtual machines (VMs) owners a continuous monitoring of their enforced SLAs metrics. Each VM integrates two agents (primary: MIB-II SNMP & secondary: SBLIM ProviderCmpiBase) providing metrics reflecting the performance level of that VM. In most large scale systems, distributed agents push metrics periodically; in our case, agents push metrics each 10 seconds to specific pre-configured autonomic managers. We assume that our studied SLA template distinguishes two time-slots: metrics must be refreshed at the client side with a freshness falling into the range of [START_REF] Roxburgh | Monitoring as an sla-oriented consumable service for saas assurance: A prototype[END_REF][START_REF] Thongtra | An adaptable capability monitoring system[END_REF][START_REF] Munawar | Adaptive monitoring with dynamic differential tracing-based diagnosis[END_REF][START_REF] Nobre | Decentralized detection of sla violations using p2p technology[END_REF] seconds during the first time-slot, and a range of [30][31][32][33][34][35][36][37][38][39][40] seconds for the second one. The SLAs metrics values are instrumented and delivered automatically through polling and exporting respectively in a manner that, once new SLA is enforced, the autonomic managers pull metrics with the lowest freshness value (3 seconds).

Objectives. Based on the data center management strategies, human administrators identify three high level monitoring goals: Respect Metrics Freshness makes sure that SLAs are monitored appropriately regarding freshness, Minimize Monitoring Cost aims at limiting the resources dedicated to monitoring as much as possible, and React to Gathering Problems operates resilient gathering mechanisms after analyzing the potential reasons of gathering problems.

Patterns. Several patterns could be exploited to refine the first objective. During the first time-slot, we use the temporal pattern to relax polling & exporting by updating their periods (Update Polling & Exporting Period in Figure 3) with respect to the highest freshness range (6 seconds). If delivering freshness violates the highest freshness, that would be a result of overloading manager [START_REF] Toueir | A goal-oriented approach for adaptive sla monitoring: a cloud provider case study[END_REF], thus we apply the spatial pattern as a second alternative, and consequently, a new autonomic manager will be deployed to assist the overloaded one (Launch Delegated Manager, Expand Perimeter & Shrink Perimeter ). As a third alternative, and in case that the overloaded autonomic manager monitors non-SLAs metrics (i.e., physical servers healthiness), the metric pattern could be applied to transfer them to other manager, in order to relax the first one (Add & Remove Aspects). Since the second time-slot freshness ([30-40] seconds) is greater than agents push period (10 seconds), there is no need to poll metrics, nor to export all received metrics. Rather, we apply the temporal pattern to update the exporting period from [START_REF] Roxburgh | Monitoring as an sla-oriented consumable service for saas assurance: A prototype[END_REF][START_REF] Thongtra | An adaptable capability monitoring system[END_REF] to [30-40] seconds (Update Exporting Period ), and consequently, the first time-slot pollings will be stopped (Deactivate Polling).

The second objective is refined using spatial pattern in order to shutdown recently deployed managers, during the first time-slot. Thus, an underloaded manager delegates its whole perimeter to another one, and shutdowns itself (Expand Perimeter, Shrink Perimeter & Shutdown Delegating Manager ). During the second time-slot, autonomic managers already deliver to clients around onethird of the metrics pushed by agents, thus no adaption actions are to be taken in regard with minimizing monitoring resources.

As the third objective deals with gathering problems (e.g., lack of collected information), autonomic managers would act on the exchange pattern. Notably, they substitute either remote agent or communication protocol (Substitute Agent & Protocol in Figure 3). But, if "failures" at the agent side cause information lack, in such case substituting protocol won't solve the problem. Therefore, we act on the metric pattern by launching troubleshooting to acquire more knowledge (Add Aspects) required to determine the appropriate substitution.

For all previous objectives, autonomic managers adapt their monitoring if they recognize adaptation stimuli. Therefore, we exploit guard pattern to apply adaptation actions (Adaptation in Figure 3, as Subgoal 2 in Table 1) as response to specific stimulus (Guard, as Subgoal 1), while maintaining the current monitoring behavior unless adaptation takes place (Unless, as Subgoal 3). resources in their SLA monitoring mode instances. Otherwise the delegated manager will activate the appropriate monitoring mode instances over the transfered resources. However, the Shrink Perimeter leaf goal is realized through withdrawing the determined managed resources from their SLA monitoring modes.

Each metric definition instance belonging to a SLA monitoring mode is systematically associated with a listening instance to update that metric with the pushed values. But, the same metric definition is also associated with a polling instance to pull its values according to the SLA period. However, Deactivate Polling leaf goal is applied on all SLA monitoring mode instances to stop polling and consequently listening continues to instrument associated metrics with the agents push period.

Both Add/Remove Aspects leaf goals are realized respectively through monitoring mode activation/deactivation. Therefore, to transfer the servers healthiness indicators, the corresponding monitoring mode will be activated on the delegated manager and deactivated from the delegating one. Rather, autonomic manager activates the "gathering troubleshooting" mode to determine whether it must substitute agent or protocol. This mode encapsulates metrics collected from remote agents SNMPv2-MIB variables, such as: snmpInBadVersions & snmpIn-BadCommunityNames4 . If the SNMP messages querying agent don't increment the "gathering troubleshooting" mode metrics of that agent, it means that the synthesized SNMP messages are correct, and probably the problem comes from the agent itself; in that case, Substitute Agent is applied.

Conclusion & Perspectives

Based on the Requirements Engineering, we proposed a goal-oriented approach for designing self-managed monitoring in autonomic systems. This approach assists human administrators to adapt the monitoring system behavior regarding quality requirements. We designed four reusable monitoring adaptation patterns falling into reconfiguration dimensions.

About the perspectives, and in order to validate our approach, we need to consider two validation levels. First, the completeness and correctness of the patterns refinement must be validated. Using proof theory, we can avoid missing some necessary requirements, and also, we can discover the available alternatives to refine a given goal [START_REF] Darimont | Formal refinement patterns for goal-driven requirements elaboration[END_REF]. On the other hand, once the monitoring system adaptation is modeled as a transition system (i.e., set of states and transitions) through the Temporal Logic, we need to determine the critical properties to be checked (i.e., invariance, safety, eventuality, fairness, and precedence) [START_REF] Goranko | Temporal logics for specification and verification[END_REF]. Using model checking, the properties satisfaction is checked, either on particular state, or path, or even the whole system model. Apart from validation, and besides enriching patterns, adaptation actions at the agent side need to be investigated, and orchestrated with those applied at the autonomic manager side.

Fig. 2 :

 2 Fig. 2: Dimensional Patterns

Table 1 :

 1 Patterns Refining Achieve Goals (P ⇒ ♦Q)[START_REF] Darimont | Formal refinement patterns for goal-driven requirements elaboration[END_REF]

	Pattern	Subgoal 1	Subgoal 2	Subgoal 3
	Milestone Pattern

Notice that Agents in networks and systems management are entities responding to management requests coming from other management entities called Managers; therefore the term "Agent" in RE has a different meaning.

In pull, entity needing information solicits the one possessing it, that responds with queried information; while in push, entity possessing information reports it to others.

For instance, the throughput is not an instantaneous metric, and the validity of its value equals to the temporal interval through which that value was measured.

These variables describe the number of SNMP messages delivered to a SNMP agent with unsupported version and unknown "community string", respectively .

Applying Leaf Goals. In our case-study, autonomic managers are CIM servers operating model-guided monitoring. Moreover, the healthiness indicators of the physical servers, as well as each SLA template, are converted into a monitoring mode each, whereby metrics and constraints are encapsulated.

The Update Polling/Exporting Period leaf goals are realized through applying the corresponding methods on all SLA monitoring mode instances to stretch the time interval of collecting/delivering metrics.

The Expand Perimeter leaf goal is realized in two manners: first, in case of delegating managed resources to another autonomic manager, whereby the same SLA templates are already enforced (i.e., SLAs monitoring modes are already instantiated), the delegated manager will join the yet transfered managed


Hlib Mykhailenko
email: hlib.mykhailenko@inria.fr

Giovanni Neglia
email: giovanni.neglia@inria.fr

Fabrice Huet
email: fabrice.huet@unice.fr

Which Metrics for Vertex-Cut Partitioning?

In this paper we focus on vertex-cut graph partitioning and we investigate how it is possible to evaluate the quality of a partition before running the computation. To this purpose we scrutinize a set of metrics proposed in literature. We carry experiments with the widely-used framework for graph processing Apache GraphX and we perform an accurate statistical analysis. Our preliminary experimental results show that communication metrics like vertex-cut and communication cost are effective predictors on most of the cases.

Introduction

Many distributed architectures have been proposed to process large graphs, like for example GraphX [START_REF] Xin | Graphx: A resilient distributed graph system on spark[END_REF], Pregel [START_REF] Malewicz | Pregel: a system for large-scale graph processing[END_REF], PowerGraph [START_REF] Gonzalez | Powergraph: Distributed graph-parallel computation on natural graphs[END_REF], and GraphLab [START_REF] Low | Distributed GraphLab: a framework for machine learning and data mining in the cloud[END_REF]. Before starting the actual computation, the graph is partitioned among the executors. Each executor then processes only a subset of the graph, but it needs to periodically share intermediate-computation results with the other executors. Partitioning affects then the total execution time in two different ways. First, the executor computation time is related to the size of the subset it needs to process. Second, the communication time depends on how much overlap there is among the different subsets. As a trivial example, if the graph has different components and each executor gets assigned one of them, executors can in many cases work almost independently. While partitioning can significantly affect the total execution time, finding the best partition for a given graph, cluster setting and application task can be itself a hard computation problem. In practice, existing partitioners rely on heuristics to find good enough trade-offs between balance (subsets should have roughly the same size), and communication (subsets should be as disjoint as possible).

A classic way to distribute the graph is the edge-cut approach where each vertex is assigned to a different partition. An edge is cut, if its vertices belong to two different partitions. In most applications, the executor storing one of the two vertices, say it the source, will need to be updated when the destination value is modified. Edge-cut partioners try then to minimize the number of edges cut. More recently, vertex-cut partitioning has been advocated as a better approach to process graphs with a power-law degree distribution (very common in real-world datasets) [START_REF] Gonzalez | Powergraph: Distributed graph-parallel computation on natural graphs[END_REF]. In this case, edges are mapped to partitions and vertices are cut if their edges happen to be assigned to different partitions. The effect can be qualitatively explained with the presence in a power-law graph of hubs, i.e. nodes with degree much larger than the average. In an edge-cut partition, attributing a hub to a given partition easily leads to i) computation unbalance, if its neighbors are also assigned to the same partition, or ii) to a large number of edges cut and then strong communication requirements. A vertex-cut partitioner may instead achieve a better trade-off, by cutting only a limited number of hubs. Analytical support to these findings is in [START_REF] Bourse | Balanced graph edge partition[END_REF]. For this reason many new graph computation frameworks, like GraphX [START_REF] Xin | Graphx: A resilient distributed graph system on spark[END_REF] and PowerGraph [START_REF] Gonzalez | Powergraph: Distributed graph-parallel computation on natural graphs[END_REF], rely on vertex-cut partitioners.

In this paper we focus on vertex-cut partitioning and we investigate how it is possible to evaluate the quality of a partition before running the computation. Answering this question would be useful both for the data analysts who need to choose the partitioner appropriate for their graphs and for developers who aim at proposing more efficient partitioners. In this paper we make a step towards providing an answer by analysing a set of metrics proposed in literature to evaluate the intrinsic quality of a partition. To this purpose we carry experiments with the widely-used framework for graph processing Apache GraphX and perform an accurate statistical analysis.

The paper is organized as follows. Apache GraphX is presented in Section 2 together with the metrics for graph partitioning which we consider in our study. In Section 3 we discuss the graph partitioning algorithms implemented in GraphX. In Section 4 we discuss our statistical methodology to study the dependency between the execution time of different graph-processing algorithms and partitioning metrics. We discuss the experimental setup and illustrate the results in Section 5. Finally, we conclude in Section 6.

The source code is available here [START_REF]Source code[END_REF].

Background

We start this section by presenting Apache GraphX. Then, we discuss the metrics used to evaluate the quality of graph partitions.

Apache GraphX

Apache Spark [START_REF]Apache Spark[END_REF] is a large-scale distributed framework for data processing. It relies on the concept of Resilient Distributed Dataset [START_REF] Zaharia | Resilient distributed datasets: A fault-tolerant abstraction for in-memory cluster computing[END_REF] (RDD) which is an immutable distributed collection mainly stored in RAM memory. The basic API provides methods for mapping and reducing RDDs.

Apache GraphX [START_REF] Xin | Graphx: A resilient distributed graph system on spark[END_REF] is a widely-used Spark library for graph processing. GraphX programs are executed as Spark jobs. The framework relies on a vertexcut partitioning approach where a vertex can be split among several partitions. There exist many different algorithms for choosing where to map a given edge. We discuss in Section 3 six such algorithms implemented in GraphX.

By default, GraphX considers all input graphs as directed ones, but then an algorithm can work on its input as if it were undirected. This is for example the case of the Connected Components algorithm built-in in GraphX. Other algorithms, like GraphX implementation of PageRank [START_REF] Page | The pagerank citation ranking: bringing order to the web[END_REF], instead assume the input is directed. In this case, if one wants to process an undirected graph, he/she can modify directly the code or pre-processing the input and replace each undirected edge by two edges with opposite directions.

GraphX is executed on a cluster with one master machine and N m machines. A graph is stored in GrahpX as two RDDs: a vertex RDD and an edge RDD. The vertex RDD stores the ids of the vertices and their values. The edge RDD stores source ids of the edges, destination ids of the edges, and the values assigned to the edges. Each of these RDDs is split in N subsets each assigned to a different executor 1 . The vertex RDD is always partitioned by a hash function based on vertex ids, while the edge RDD is partitioned using a user-specified partitioner. Note that the number of partitions N can be different from the number N m of machines (for example a machine can be assigned multiple executors). GraphX distributes the N partitions among the machines in a round robin fashion.

Metrics for partitioning

We introduce here some additional notation to simplify the definitions.

Let us denote the input directed graph as G = (V, E), where V is a set of vertices and E is a set of edges. Vertex-cut partitioning splits the set of edges into N disjoint subsets, E 1 , E 2 , . . . , E N . With some abuse of terminology, we also say that a vertex belongs to the (edge) partition E i , if it is the source or the destination of an edge in E i . Using this convention, an edge can only belong to one partition, but a vertex is at least in one partition and at most in N partitions. Let V (E i ) denote the set of vertices that belong to partition E i . The vertices that appear in more than one partition are called frontier vertices. Each frontier vertex has then been cut at least once. F (E i ) denotes the set of frontier vertices that are inside partition E i . F (E i ) is the set of vertices in partition E i that were not cut, and this set equals to F (

E i ) = V (E i ) -F (E i ).
There exist several metrics to evaluate the quality of a graph partition. In this paper we want to understand how well they fulfill this goal, i.e. to which extent they can be used to predict the final performance. In our study we considered six metrics which were used by the authors of graph partitioning algorithms presented in [START_REF] Rahimian | Distributed Vertex-Cut Partitioning[END_REF] and [START_REF] Guerrieri | Dfep: Distributed funding-based edge partitioning[END_REF]. The first three of them aim to quantify how homogeneous are the sizes of the subsets in the partition and then how balanced are the computation tasks across all the executors. The last three instead target the communication cost, by considering how many vertices are shared across the different subsets.

Balance (denoted as BAL) is the ratio between the maximum number of edges in one partition and the average number of edges in all partitions, BAL =

max i=1,...N (|Ei|) |E|/N . NSD is the normalized standard deviation of partition size, NSD = N i=1 |Ei| |E|/N -1 2 1 N .
Largest partition (denoted as LP) shows the number of vertices in the largest partition of a graph, LP = max i=1,...N |V (E i )|.

Vertex-cut (denoted as VC) shows the number of vertices that were cut

during graph partitioning, VC = |V| - N i=1 | F (E i )|.
Replication factor (denoted as RF) is the ratio between the total number of vertices in a partitioned graph and the original graph, RF =

= N i=1 |V (E i )| 1 |V| = N i=1 (| F (E i )| + |F (E i )|) 1 |V| = 1 -VC |V| + CC |V| .
For this reason it is sufficient to consider only 2 of them. We selected VC and CC. 2GraphX does not have a built-in function which provides values of these metrics, this is why we have implemented new GraphX functions to compute them [START_REF]Source code[END_REF].

GraphX partitioners

Six partitioning algorithms are considered in this work and described in this section. The first four partitioners given below are the GraphX built-in ones. The fifth (HybridCut) was proposed in [START_REF] Chen | PowerLyra: Differentiated Graph Computation and Partitioning on Skewed Graphs[END_REF] and implemented for GraphX by Larry Xiao [START_REF]HybridCut and HybridCutPlus partitioners[END_REF] together with its variant HybridCutPlus.

RandomVertexCut (denoted as RVC) partitioner randomly assigns edges to partitions to achieve good balance (with high probability). In particular it computes the hash value for each pair (source vertex id, destination vertex id). The hash space is partitioned in N sets with the same size. As a result, in a multigraph all edges with the same source and destination vertices are assigned to the same partition. On the contrary, two edges among the same nodes but with opposite directions belong in general to different partitions.

CanonicalRandomVertexCut (denoted as CRVC) algorithm is similar to the previous partitioner but in this case the two ids are ordered before calculating the hash value. Hence, all edges between two vertices are allocated to the same partition, regardless of their direction.

EdgePartition1D (denoted as EP1D) targets a small replication factor by keeping all outgoing edges from one vertex in a single partition. 3 More precisely, the algorithm assigns each edge based on the hash value of its source vertex id.

EdgePartition2D (denoted as EP2D) guarantees a bound on the replication factor metric and it can provide also good balance if the number of partitions is appropriately chosen. The algorithm splits a set of vertices in M = √ N subsets S 1 , S 2 , . . . , S M using a hash function. Then it is possible to partition the adjacency matrix in M 2 blocks, where the block (i, j) contains all edges with a source in S i and a destination in S j . Finally, the M 2 blocks are assigned to partitions in a round robin fashion. A vertex can then appear in at most 2M -1 ≤ 2( √ N +1) partitions. If N is chosen to be a square number, partitions will be on average well balanced.

HybridCut (denoted as HC) tries to reduce the number of frontier vertices without sacrificing the balance. The basic idea is that hubs will be cut in any case to avoid unbalance, while the edges of the nodes with small degree can be kept in the same partition. In particular, the algorithm considers the destination vertex of each edge. If this vertex is a hub, i.e. if its in-degree is larger than a given threshold, then the edge is assigned to a partition based on the hash value of its source id. Otherwise, the assignment is based on the destination vertex id value.

HybridCutPlus (denoted as HCP) is a hybrid partitioner which combines HC and EP2D. It computes the inDegree value (number of incoming edges) of the source and the destination vertices. If both values are either greater or less than a given threshold, then the EP2D partitioner is used for this edge. If the inDegree value of the source vertex is less than the threshold and the inDegree value of the destination vertex is greater than the threshold, then the hash value of the source is used to assign the edge to a random partition. If the inDegree value of the source vertex is greater than the threshold and the inDegree value of the destination vertex is less than the threshold the hash value of the destination is used to assign the edge to a random partition. The HCP partitioner targets balance and communication efficiency properties.

Statistical analysis methodology

We want to study if the partition metrics considered in literature and described above in Sec. 2.2 are good predictors for the final execution time. To this purpose, we would like ideally to design experiments where we could randomly and independently select the 5 different metrics, generate an input graph with these quintuple of characteristics, evaluate the execution time of the application of interest and finally use some statistical tool like regression models to identify the contribution of each metric to the execution time. Unfortunately, it is not possible to select first the partition metrics and then produce a graph with such characteristics. We can only use a partitioner on a given graph to produce a partition and then evaluate its metrics. Their values will be far from independent and in particular the structure of their correlation can be a function of the specific partitioner used.

For this reason, in order to obtain a variety of different configurations, we use all the 6 partitioners on the same graph. For each graph this leads to 6 different quintuples. A statistical model based on 6 points in a 5-dimension space would very likely overfit the data. Considering other graphs allow us to obtain other samples. The drawback is that if we use real datasets, they are going to differ for the number of nodes and of edges. Now, these two metrics are very likely to have an effect much more important on the execution time than the partition metrics, which are the focus of our study. 4 The contribution of the partition metrics would be dwarfed by the size change.

In order to overcome this difficulty, we generated 10 input graphs with the same size simply randomly permutating the ids of the vertices. Even if the different graphs are homeomorphic, partitions are calculated using the ids, and then each partitioner produces different subsets with different values for the metrics of interest. This approach allowed us to have 60 different quintuples for each graph: the ten permutations of the original graph multiplied by the six partitioners. In this way the risk of an overfitting model is significantly reduced.

To identify the most important metrics, we then used a linear regression model (denoted as LRM) with the 5 partition metrics as input variables (or predictors) and the execution time as output variable (or response). We want to find the most important predictors. Obviously the more predictors we have in a LRM the smaller is the residual error, but which ones are really important? For this purpose we used the best subset selection method [START_REF] James | An introduction to statistical learning[END_REF]Chapter 6]. In particular, given the small number of predictors [START_REF] Akaike | Akaikes information criterion[END_REF], we were able to consider all the possible 5! linear models for each original graph. Models with the same number of predictors can be easily compared through their R 2 value. In this way 5 models were selected, the best one was finally identified as the one with the smallest Akaike information criterion [START_REF] Akaike | Akaikes information criterion[END_REF]. 5 Once the best model was selected, we ordered its predictors by considering those that lead to the largest increase of the R 2 when added as predictors.

Experiments

All our experiments were performed on a cluster of nodes with dual-Xeon E5-2680 v2 @2.80GHz with 192GB RAM and 10 cores (20 threads). We used Spark version 1.4.0 and Spark standalone cluster as a resource manager. We used two different cluster configurations. In the first configuration, master and executor share the same machine. In the second case, one machine is dedicated to the master and one to each of ten executors. We configured 4 Spark properties: spark.executor.memory equals to 4 GB (if there are 10 executors) or 40 GB (if there is 1 executor), spark.driver.memory equals to 10 GB, spark.cores.max equals to 10, spark.local.dir points to local directory.

Two different processing algorithms were evaluated: Connected Components and PageRank with 10 iterations. PageRank exhibits the same computation and communication pattern at each stage (the PageRank of each vertex value is updated according to the same formula and then propagated to the neighbors), while Connected Components implements a label propagation mechanism, where as time goes on, less updates are required.

As datasets, we used two undirected graphs: the com-Youtube graph (1,134,890 vertices/2,987,624 edges) and com-Orkut graph (3,072,441 vertices/117,185,083 edges) from SNAP project [START_REF]Stanford Large Network Dataset Collection[END_REF].

For each experiment we, first, randomly selected a partitioner, a graph processing algorithm, and a permutation of an input graph. Second, we applied the selected partitioner to the input graph. Finally, we executed the selected graph processing algorithm. To overcome execution time variability (due to a shared network, shared distributed file system, operational system layer, etc.), every combination has been tested at least 30 times, obtaining 95% confidence intervals for the execution time whose relative amplitude is always less than 10%.

As we mentioned above, for each graph we obtained 60 different quintuples for the partition metrics. Table 1 shows their correlation matrix for com-Youtube graph. The matrix clearly identifies the two groups of metrics: those relative to partition balance (BAL, LP NSTDEV) and those related to communication (VC, CC). The high level of collinearity also explains why in the LRM some coefficients are negative [START_REF] James | An introduction to statistical learning[END_REF]Chapter 3].

We computed and selected linear regression models for com-Youtube, com-Youtube doubled, and com-Orkut graph, as it was discussed in Section 4. Table 2 and Table 3 summarize our experimental results respectively for PageRank and Connected Components. For each input graph and number of machines, they show the best LRM with the factors listed in decreasing order of importance. We note that, although a single node has enough resources to perform all the computation, the results in Table 2 and Table 3 show that the average execution time decreases as the number of machines increases. A similar result was observed in [START_REF] Hood | Performance impact of resource contention in multicore systems[END_REF] where it was due to the higher memory contention in case of a Below we illustrate our main findings about the partition metrics.

The execution time depends on the partition metrics

For PageRank the R 2 value is very high (always above 0.99) while the Root Mean Square Error (RMSE) is less than 5% of mean execution time. This leads to conclude that indeed the execution time depends on these metrics. The dependency is less stronger for Connected Components. This can be explained with the fact that PageRank roughly does the same operations at each stage, while Connected Components is a label propagation algorithm where very soon the values of most of the vertices will not be updated. The execution time likely depends on graph properties (e.g. the graph diameter) that are not well captured by these metrics. Figures 1 and2 confirm these results by comparing the experimental execution time with the one predicted using the best LRM model we found. Results are grouped by partitioners, for each of them results for the 10 graph permutations are shown. The plot for the com-Orkut is similar (we do not show it here due to lack of space).

En passant we observe that HC appears to be the best partitioner confirming the results in [START_REF] Chen | PowerLyra: Differentiated Graph Computation and Partitioning on Skewed Graphs[END_REF], it also outperforms its variant HCP.

Communication metrics are the most important

In all the LRMs, the most important metric is always a communication metric (CC or VC), even when a single machine is used and then the network is scarcely used.

This result is also confirmed by comparing the 5 single-predictor LRMs: the LRMs using VC or CC have larger R 2 value. Figure 3 shows how LRM models using only VC or CC are able to produce quite good predictions at least of the relative performance of the different partitioners. On the contrary, the corresponding plot for balance metrics (BAL, LP, NSTDEV) shows an almost horizontal line.

As a side note, we tried to find a setting where balance metrics would have become the most important ones. In particular we considered i) a single machine, ii) very limited memory for each partitioner (70MB), iii) a modified version of PageRank where each computation is effectuated 1000 more times. The purpose was to make communication among executors as fast as possible, while increasing the computational burden on each executor. Even in this settings, communication metrics appeared to be the most important ones.

Results are robust to the partitioners considered

One of the purposes of our analysis is to be able to rank a priori different partitioners before carrying on the experiments. We evaluated how much our results are sensitive to the specific set of partitioners used to tune the LRM. To this purpose, we carried on the same analysis using only 5 partitioners (we did not use HC). We then used the new LRM model to predict the execution time for the HC partitioner. Figure 4 shows that predicted execution time for HC is underestimated. However, the order of the partitioners in terms of execution time is the correct one. Then the LRM model correctly estimates that using HC the application would finish faster. One could argue that the LRM has been trained using HCP that is a variant of HC. We then performed the same set of experiments removing both HC and HCP. In this case the predicted executed time is really off, but the LRM still correctly predicts the ranking of the different partitioners.

There is space for better partition metrics

While the results above are encouraging, the following experiment suggests that these partition metrics do not capture all the relevant features.

We considered a version of com-Youtube graph where we doubled each edge so that both directions are present in the input graph (and then PageRank correctly computes the PageRank of the original undirected graph). We denote the doubled graph as com-Youtube doubled.

The execution time of both algorithms on com-Youtube doubled is larger, roughly 30% (see Table 4). Nevertheless, the communication partition metrics that we identified as the most important ones may fail to detect the change. Indeed, this is evident if we compare the CRVC partitioner with the other ones. Indeed, as shown in Table 4 for a specific input graph, the CRVC partitioner provides the same values for all the partition metrics, and in particular for VC and CC metrics, both for com-Youtube doubled and com-Youtube, since it partitions edges based on the vertex (either source or destination) with the smallest id. For the other partitioners instead, the metrics VC and CC are different. The difference is evident if we compute the LRM including or not CRVC. In the first case the R 2 value of the LRM for PageRank significantly decreases from 0.99 for com-Youtube (see Table 2) to 0.82. In the second case the LRMs are equally good.

We think a more meaningful communication metric could be defined that would permit to identify the difference between CRVC and the other partitioners. We plan to investigate this aspect in the future.

Conclusion

We used linear regression models with partitioning metrics as predictors and the average execution time for different graph processing algorithms as the observed values. The obtained models confirmed that there is actual dependency between these quantities. More importantly, the most important metrics are CC and VC, the both are an indicator of the amount of communication among the executors. On the contrary, the metrics that quantify load unbalance across the executors are less important. This conclusion holds whether communication is inter-machine or intra-machine, i.e. both if the network is used or not. Our results are robust to the original set of partitioners used to train the model, and the model can correctly rank other partitioners. At the same time, there is probably still space to define more useful metrics to evaluate partions.

  i )| 1 |V| . Communication cost (denoted as CC) is the total number of frontier vertices among all the partitions, CC = N i=1 |F (E i )|.It should be noticed that RF, CC, and VC metrics are linear combinations of each others: RF

Figure 1 :

 1 Figure 1: Execution time for PageRank algorithm on com-Youtube graph: experimental results vs best LRM predictions

Figure 2 :

 2 Figure 2: Execution time for Connected Components algorithm on com-Youtube graph : experimental results vs best LRM predictions

Figure 3 :

 3 Figure 3: Execution time for PageRank algorithm on com-Youtube graph: experimental results vs LRM predictions using communication metrics

Figure 4 :

 4 Figure 4: Prediction for HC partitioner (using com-Youtube graph, and PageRank algorithm)

Table 1 :

 1 Correlation matrix for partition metrics com-Youtube It is possible that the effect has the same cause here.

		BAL	LP	NSD	VC	CC
	BAL	1	0.9332 0.9445 0.2070 -0.2448
	LP		1	0.9799	0.0993 -0.3177
	NSD			1	0.1275 -0.2667
	VC				1	0.8737
	CC					1
	single machine.

Table 2 :

 2 Results for PageRank algorithm

	Graph name Cluster configuration	Mean time (ms)			Metrics (ordered by importance)	RMSE	R 2
	com-Youtube	single machine 1 master + 10 executors		12,343 9,068	metrics coefficients 1.5507 0.3156 0.42743 -0.49504 0.33548 VC CC BAL LP NSD metrics VC CC BAL LP NSD coefficients 1.2816 0.27261 0.50989 -0.64857 0.49891	111.83 0.996 125.47 0.993
	com-Orkut	single machine 1 master + 10 executors		122,929 94,738	metrics coefficients 2.312 CC metrics CC coefficients 1.3285 8.3243 LP 5.039 LP	VC 12.942 NSD -11.486	VC 12.881	1199.4 0.998 2222.8 0.989
			10.5 11	observed predicted by best LRM
			10
		Execution time (s)	7.5 8 8.5 9 9.5
			7
			6.5
			6
				H C	H C P	E P 1 D	E P 2 D	R V C	C R V C

Table 3 :

 3 Results for Connected Components algorithm

	Graph name Cluster configuration	Mean time (ms)	Metrics (ordered by importance) RMSE	R 2
	com-Youtube	single machine 1 master + 10 executors	12,429 9,635	metrics coefficients 1.7884 -1.3177 VC LP metrics VC LP coefficients 2.1131 -1.7528	NSD 1.3861 NSD 2.1043	1119.2 0.499 922.27 0.681
	com-Orkut	single machine 1 master + 10 executors	86,253 71,952	metrics coefficients 1.2039 CC metrics CC coefficients 1.0854 15.484 LP 9.771 LP	NSD -8.2634 NSD -17.96	2628.2 0.959 2235,6 0.973

Table 4 :

 4 Metrics and execution time of PageRank for HC and CRVC

	Graph name	Average time (ms)	HC VC	CC	CRVC time (ms) Average VC	CC
	com-Youtube	6468	0.304 1.091	10327	0.452 1.831
	com-Youtube doubled	8121	0.702 2.204	13404	0.452 1.831

It is also possible to assign multiple partitions to the same executor, but we do not consider this possibility.

We decided to sacrifice RF because including RF leads easily to negative coefficients in the linear regression model.

Still the vertex can appear in multiple partitions because of its incoming edges.

Remember that our goal is choosing the best partitioner for a given input graph of interest, whose size in terms of number of nodes and edges is out of our control.

We looked at other indices (Bayesian Information Criteria, sample-size corrected AIC, R 2 adjusted) and in most cases the results did not differ.

Acknowledgements

This work was partly funded by the French Government (National Research Agency, ANR) through the Investments for the Future Program reference #ANR-11-LABX-0031-01.