Multivariate decompositions with commuting projections

S.-H. dec. is an example of multivariate decompositions obtained with (a class of) commuting projections P 1 , . . . , P d [START_REF] Kuo | On decompositions of multivariate functions[END_REF]), here orthogonals:

P j (g)(x) = g(x)dµ j (x j ) = E[g(X )|X -j = x -j ]
The form of the decomposition is simply obtained by expansion: The non-overlapping condition is written here P i (g I ) = 0 for all i ∈ I. We find again that Π I is an orthogonal projection.

O. Roustant (EMSE) Sobol-Hoeffding decomposition: bounds and extremes 2018 December 11 7 / 50

An example: separable functions Consider g(x) = f 1 (x 1 ) . . . f d (x d ), and denote m j = E(X j ). Then:

g I (x I ) = i∈I (f i (x i ) -m i ) j / ∈I m j g tot I (x) = i∈I (f i (x i ) -m i ) j / ∈I f j (x j )
Proof. The Sobol-Hoedding decomposition is obtained by expanding:

g(x) = ((f 1 (x 1 ) -m 1 ) + m 1 ) . . . ((f d (x d ) -m d ) + m d )
For each bracket, for g I , choose (f i (x i ) -m i ) if i ∈ I, and m j otherwise for g tot I , choose

(f i (x i ) -m i ) if i ∈ I O.
ν i = ∂g(x) ∂x i 2 dµ(x), ν I = ∂ |I| g(x) ∂x I

Usage for screening

Assume that:

g is continuous on ∆ = [0, 1] d for all i, the support of µ i contains [0, 1] Variable screening If either D tot i = 0 or ν i = 0, then X i is non influential O.

Usage for screening

Assume that:

g is continuous on ∆ = [0, 1] d for all i, the support of µ i contains [0, 1]

Variable screening

If either D tot i = 0 or ν i = 0, then X i is non influential

Interaction screening

If either D tot i,j = 0 or ν i,j = 0, then (x i , x j ) → g(x) is additive

Total interactions can be visualized on the FANOVA graph, where the edge size is proportionnal to the index value.

Illustration on a toy example

8D g-Sobol function, with uniform inputs on [0, 1]: 1, 4.5, 9, 99, 99, 99, 99).

g(x) = 8 j=1 |4x j -2| + a j 1 + a j with a = c(0,
O. Roustant (EMSE) Sobol-Hoeffding decomposition: bounds and extremes 2018 December 11 Illustration on a toy example 8D g-Sobol function, with uniform inputs on [0, 1]: 1, 4.5, 9, 99, 99, 99, 99).

g(x) = 8 j=1 |4x j -2| + a j 1 + a j with a = c(0,

Illustration on a toy example

A 6D block-additive function, with uniform inputs on [-1, 1]:

g(x) = cos([1, x 1 , x 2 , x 3 ] β) + sin([1, x 4 , x 5 , x 6 ] γ))
with β = (-0.8, -1.1, 1.1, 1) and γ = (-0.5, 0.9, 1, -1.1) .

O.

Illustration on a toy example

A 6D block-additive function, with uniform inputs on [-1, 1]:

g(x) = cos([1, x 1 , x 2 , x 3 ] β) + sin([1, x 4 , x 5 , x 6 ] γ))
with β = (-0.8, -1.1, 1.1, 1) and γ = (-0.5, 0.9, 1, -1.1) .

Figure:

1st order analysis (left) and 2nd order analysis (right) with 10 5 simulated data

Poincaré inequality

Poincaré inequality (1-dimensional case)

A distribution µ satisfies a Poincaré inequality if the energy in L 2 (µ) sense of any centered function is controlled by the energy of its derivative:

For all h in L 2 (µ) such that h(x)dµ(x) = 0, and h (x) ∈ L 2 (µ):

h(x) 2 dµ(x) ≤ C(µ) h (x) 2 dµ(x)
The best constant is denoted C P (µ).

O. If µ i and µ j admit a Poincaré inequality, then: If µ i and µ j admit a Poincaré inequality, then:

D i ≤ D tot i ≤ C(µ i )ν i , D i,j ≤ D tot i,j ≤ C(µ i )C(µ j )ν i,j O.
D i ≤ D tot i ≤ C(µ i )ν i , D i,j ≤ D tot i,j ≤ C(µ i )C(µ j )ν i,j
Proof 1. Denote g tot i (x) := J⊇{i} g J (x J ). Then: If µ i and µ j admit a Poincaré inequality, then:

∂g(x) ∂x i = ∂g tot i (x) ∂x i O. Roustant (
D i ≤ D tot i ≤ C(µ i )ν i , D i,j ≤ D tot i,j ≤ C(µ i )C(µ j )ν i,j
Proof 1. Denote g tot i (x) := J⊇{i} g J (x J ). Then:

∂g(x) ∂x i = ∂g tot i (x) ∂x i D tot i = Var(g tot i (x)) = g tot i (x) 2 dµ(x) ≤ C(µ i ) ∂g tot i (x) ∂x i

Link between total Sobol indices and DGSM

Theorem [START_REF] Lamboni | Derivative-based global sensitivity measures: General links with Sobol' indices and numerical tests[END_REF], [START_REF] Roustant | Crossed-derivative based sensitivity measures for interaction screening[END_REF] If µ i and µ j admit a Poincaré inequality, then: If µ i and µ j admit a Poincaré inequality, then:

D i ≤ D tot i ≤ C(µ i )ν i , D i,j ≤ D tot i,j ≤ C(µ i )C(µ j )ν i,j O.
D i ≤ D tot i ≤ C(µ i )ν i , D i,j ≤ D tot i,j ≤ C(µ i )C(µ j )ν i,j
Proof 2. Denote g tot i,j (x) := J⊇{i,j} g J (x J ). Then: If µ i and µ j admit a Poincaré inequality, then:

∂ 2 g(x) ∂x i ∂x j = ∂ 2 g tot i,j (x) ∂x i ∂x j O. Roustant (EMSE) Sobol-
D i ≤ D tot i ≤ C(µ i )ν i , D i,j ≤ D tot i,j ≤ C(µ i )C(µ j )ν i,j
Proof 2. Denote g tot i,j (x) := J⊇{i,j} g J (x J ). Then:

∂ 2 g(x) ∂x i ∂x j = ∂ 2 g tot i,j (x) ∂x i ∂x j D tot i,j = Var(g tot i,j (x)) = g tot i,j (x) 2 dµ(x) ≤ C(µ i ) ∂g tot i,j (x) ∂x i 2 dµ(x) ≤ C(µ i )C(µ j ) ∂ ∂x j ∂g tot i,j (x) ∂x i

Getting optimal Poincaré constants on intervals

Assume that dµ 1 (t)/dt = e -V (t) > 0 on a bounded interval [a, b]. Then, the smallest Poincaré constant C(µ 1 ) is obtained by solving a spectral problem:

Lf := f -V f = -λf with f (a) = f (b) = 0 Comments.
For some (rare) pdf, C(µ 1 ) can be computed semi-analytically.

For many other ones, a finite element method can be used.

Adaptations are possible for unbounded intervals and pdf vanishing at the boundaries.

See technical details in [START_REF] Roustant | Poincaré inequalities on intervals -application to sensitivity analysis[END_REF].

O. 1 output: maximal annual overflow (in meters), denoted by S:

Uniform [a, b] (b -a) 2 /π 2 cos π(x-a) b-a N (µ, σ 2 ) R σ 2 x -µ [r n,i , r n,i+1 ] 1/(n + 1) H n+1 (x) [a, b] related to Kummer hypergeom. func. Db. exp. e -|x| dx/2 R 4 × (*) [a, b], ab > 0 1 4 + ω 2 -1 e x/2 cos(ωx + φ) (*, **) [a, b], ab ≤ 0 > 1 4 + ω 2 -1 e |x|/2 × trig. spline Logistic e x (1+e x ) 2 dx R 4 × Triangular [-1, 1] ≈ 0.1729 linked to Bessel J 0 (*) For the truncated Exponential on [a, b] ⊆ R + , we use ω = π/(b -a) (**) If a < 0 < b,
S = Z v + H -H d -C b with H =   Q BK s Zm-Zv L   0.6
where H is the maximal annual height of the river (in meters).

O.

Principle

Without loss of generality, assume g 0 = 0. Define:

F 1 = {g ∈ L 2 (µ) s.t. g = g 1 } functions depending exactly on x 1 F tot 1 = {g ∈ L 2 (µ) s.t. g = g tot 1 }
functions depending at least on x 1

Notice that g 1 and g tot 1 are obtained from g by orthogonal projection

g 1 = Π F 1 (g) = E[g(X )|X 1 = .] g tot 1 = Π F tot 1 (g) = g -E[g(X )|X 2 = ., . . . , X d = .] Hence, D 1 = Π F 1 (g) 2 and D tot 1 = Π F tot 1 (g) 2 .
Lower bounds of D 1 , D tot 1 are obtained by projecting onto subspaces of F 1 , F tot

Main result

Let φ 1 , . . . , φ m be orthonormal functions in F tot 1 . Then:

D tot 1 ≥ m j=1 g(x)φ j (x)dµ(x) 2 with equality iff g has the form g(x) = m j=1 α m φ m (x) + h(x 2 , . . . , x m
). If all the φ j 's belong to F 1 then the lower bound is for D 1 .

Proof.

D tot 1 = g tot 1 2 = Π F tot 1 (g) 2 ≥ Π φ 1 ,...,φm (g) 2 = m j=1 ( g, φ j ) 2
Equality is when g tot 1 = Π φ 1 ,...,φm (g), leading to the condition above. Same arguments when all the φ j 's are in

F 1 O. Roustant (EMSE)
Sobol-Hoeffding decomposition: bounds and extremes 2018 December 11

Tensor-based lower bounds

For all j, let ψ j,0 = 1, ψ j,1 , . . . , ψ j,n j -1 be orthonormal functions in L 2 (µ j ).

Consider tensors, i.e. separable functions:

φ (x) = d j=1 ψ j, j (x j ) where = ( 1 , . . . , d ) is a multi-index.
Let T 1 = { s.t. 1 ≥ 1}, the set of tensors φ involving x 1 . Then:

D tot 1 (f ) ≥ ∈T 1 f (x)φ (x)ν(dx) 2 with equality iff f has the form f (x) = ∈T 1 α φ (x) + g(x 2 , . . . , x d ).
O. Roustant (EMSE) Sobol-Hoeffding decomposition: bounds and extremes 2018 December 11 28 / 50

Tensor-based lower bounds

As an illustration, if µ i admit the first two moments, denote:

ψ i (x) = (x i -m i )/s i
where m i is the mean and s i the s.d.

Then ψ 1 , ψ 1 ψ 2 , . . . , ψ 1 ψ j are orthonormal functions of F tot 1 . Hence:

D tot 1 ≥ g(x)ψ 1 (x)dµ(x) 2 lower bound for D 1 + m j=2 g(x)ψ 1 (x)ψ j (x)dµ(x)

Derivative-based lower bounds

All the integrals above can involve derivatives by integrating by part. But this often induce weights; Here is a partial solution to avoid weights.

Assume that µ j is continuous with pdf p j ∈ H 1 (µ j ) vanishing at the boundaries but not inside, and such that p j ≡ 0 and p j /p j ∈ L 2 (µ j ). Denote:

Z j (X j ) = (ln p j ) (X j ), I j = Var(Z j (X j ))
Then:

D tot 1 ≥ I -1 1 c 2 1 lower bound for D 1 + I -1 1 d j=2 I -1 j c 2 1,j with c 1 = g(x)Z 1 (x 1 )dµ(x) = - ∂g(x) ∂x 1 dµ(x) c 1,j = g(x)Z 1 (x 1 )Z j (x j )dµ(x) = - ∂g(x) ∂x 1 Z j (x j )dµ(x) = ∂ 2 g(x) ∂x 1 ∂x j dµ(x)

Derivative-based lower bounds: examples

For normal variables N(m j , s 2 j ):

D tot 1 ≥ s 2 1 ∂g(x) ∂x 1 dµ(x) 2 lower bound for D 1 + s 2 1 d j=2 s 2 j ∂ 2 g(x) ∂x 1 ∂x j dµ(x) 2 Dist. name Support p Z I Normal R 1 s √ 2π exp -1 2 (x-m) 2 s 2 -(X -m)/s 2 1/s 2 Laplace R 1 2s exp |x-m| s -sgn(X -m)/s 1/s 2 Cauchy R 1 π s (x-x 0 ) 2 +s 2 -2(x-x 0 ) (x-x 0 ) 2 +s 2 1/(2s 2 ) O. Roustant (EMSE)
Sobol-Hoeffding decomposition: bounds and extremes 2018 December 11 31 / 50

Improvements on existing works

According to results given in the review [START_REF] Kucherenko | Derivative-Based Global Sensitivity Measures[END_REF],

For normal distributions, we improve on:

D tot 1 ≥ D 1 ≥ s 2 1 ∂g(x) ∂x 1 dµ(x) 2 . O.

Improvements on existing works

According to results given in the review [START_REF] Kucherenko | Derivative-Based Global Sensitivity Measures[END_REF],

For normal distributions, we improve on:

D tot 1 ≥ D 1 ≥ s 2 1 ∂g(x) ∂x 1 dµ(x) 2 .
For uniforms on [0, 1] using the orthonormal function obtained from x m 1 , and an integration by part, we obtain:

D tot 1 ≥ D 1 ≥ 2m + 1 m 2 (g(1, x -1 ) -g(x))dx -w (m+1) 1 2
where w

(m+1) 1 = ∂g(x) ∂x 1 x m+1 1 dx.

Improvements on existing works

According to results given in the review [START_REF] Kucherenko | Derivative-Based Global Sensitivity Measures[END_REF],

For normal distributions, we improve on:

D tot 1 ≥ D 1 ≥ s 2 1 ∂g(x) ∂x 1 dµ(x) 2 .
For uniforms on [0, 1] using the orthonormal function obtained from x m 1 , and an integration by part, we obtain:

D tot 1 ≥ D 1 ≥ 2m + 1 m 2 (g(1, x -1 ) -g(x))dx -w (m+1) 1 2
where w

(m+1) 1 = ∂g(x) ∂x 1 x m+1 1 dx.

When using derivatives and other numerical considerations

We must compute squared integrals θ = ( h(x)dµ(x)) 2 , when h has the form:

h dir = gφ 1 ,

When using derivatives and other numerical considerations

We must compute squared integrals θ = ( h(x)dµ(x)) 2 , when h has the form:

h dir = gφ 1 , gφ 1 φ j , . . ., or h der = ∂g ∂x i , ∂g ∂x j Z j , ... for centered function φ 1 , φ j , Z j .
The sample estimate θ = 1 n n i=1 h(X i ) 2 , with X 1 , . . . , X n i.i.d. ∼ µ, verifies:

θ ≈ N θ, 4θ n Var µ (h)
Hence, for one squared integral, using the derivative form can reduce estimation error when h der is less variable than h dir .

Partial conclusions

Lower bounds of a (convex comb. of) ANOVA term g I can be obtained by projection onto subspaces of its ANOVA space {g ∈ L 2 (µ) s. Using derivative-based inequalities may be useful when the derivative is less variable than the function itself.

Multivariate dependence

Denote F a multivariate cdf,

F (x) = P(X 1 ≤ x 1 , . . . , X d ≤ x d )
Assume that F is in the domain of attraction of a max-stable distribution H i.e. there exist vector sequences a n > 0, b n s.t. for indep. samples X 1 , . . . , X n of F

P max n k =1 (X k 1 ) -b n,1 a n,1 ≤ x 1 , . . . , max n k =1 (X k d ) -b n,d a n,d ≤ x d → n→∞ H (x)
O. Roustant (EMSE) Sobol-Hoeffding decomposition: bounds and extremes 2018 December 11

Multivariate dependence

Denote F a multivariate cdf,

F (x) = P(X 1 ≤ x 1 , . . . , X d ≤ x d )
Assume that F is in the domain of attraction of a max-stable distribution H i.e. there exist vector sequences a n > 0, b n s.t. for indep. samples X 1 , . . . , X n of F

P max n k =1 (X k 1 ) -b n,1 a n,1 ≤ x 1 , . . . , max n k =1 (X k d ) -b n,d a n,d ≤ x d → n→∞ H(x)
In the univariate case, H is a generalized extreme value distribution, summarizing the three types Fréchet, Weibull, Gumbel

Multivariate dependence

Denote F a multivariate cdf,

F (x) = P(X 1 ≤ x 1 , . . . , X d ≤ x d )
Assume that F is in the domain of attraction of a max-stable distribution H i.e. there exist vector sequences a n > 0, b n s.t. for indep. samples X 1 , . . . , X n of F

P max n k =1 (X k 1 ) -b n,1 a n,1 ≤ x 1 , . . . , max n k =1 (X k d ) -b n,d a n,d ≤ x d → n→∞ H(x)
In the univariate case, H is a generalized extreme value distribution, summarizing the three types Fréchet, Weibull, Gumbel

In the multivariate case, the margins are gevd, and the multivariate dependence is characterized by a multivariate function extreme value copula, stable tail dependence function, . . .

Asymptotic independence and tail dependograph

Let A, B a partition of {1, . . . , d} X A and X B are asymptotically independent

⇔ H(x) if of the form H(x) = H A (x A )H B (x B ) ⇔ (u) if of the form (u) = (u A ) + (u B ) ⇔ ∀i ∈ A, ∀j ∈ B,

Inference

The formula (u) = lim z→+∞ z 1 -F F -1 1 (u 1 /z), . . . , F -1 d (u d /z) leads to the natural estimator ( [START_REF] Huang | Statistics of bivariate extremes[END_REF])

ˆ k ,n (u) = n k 1 - 1 n n s=1 1 X (1) s < X (1) n-[ku 1 ]+1,n , . . . , X (d) s < X (d) n-[ku d ]+1,n = n k 1 - 1 n n s=1 1 u 1 < R(1) s , . . . , u d < R(d) s = n k - 1 k n s=1 d t=1 1{u t < R(t) s }
separable function with:

X (t) 1,n , . . . , X (t)
n,n : sorted data (asc. order) for coordinate t R(t)

s := n-R (t) s +1 k , where R (t) s is the rank of X (t) s among X (t) 1 , . . . , X (t)
n .

Inference

Let µ = µ 1 ⊗ • • • ⊗ µ d a measure on [0, 1] d (without special link with F ).
As a sum of separable functions, the whole Sobol-Hoeffding decomposition of the stdf estimator can be computed in closed form, and in particular

ˆ tot k ,n;{i,j} (u) = - 1 k n s=1 d t=1 1{u t < R(t) s } -1 {t∈{i,j}} µ t R (t)
s and the tail dependograph as well where Y ;I is some Gaussian process.

D tot {i,j} ( ˆ k ,n ) = 1 k 2 n s=1 n s =1 d t=1 µ t R(t) s ∧ R(t) s -1 {t∈{i,j}} µ t R(t) s µ t R (

  I d = (P 1 + (I d -P 1 )) . . . (P d + (I d -P d )

  Partial variances: D I = Var(g I (X I )), and Sobol indices S I = D I /D Partial variances: D I = Var(g I (X I )), and Sobol indices S I = D I /D

Figure :

 : Figure: 1st order analysis (left) and 2nd order analysis (right) with 10 5 simulated data O. Roustant (EMSE) Sobol-Hoeffding decomposition: bounds and extremes 2018 December 11 11 / 50

Figure :

 : Figure: Poincaré constant of µ = N (0, 1) truncated on I = [-b, b], vs µ(I)

  t. g = g I } → Illustrated on main and total effects, but very general! a (convex comb. of) ANOVA term g I can be obtained by projection onto subspaces of its ANOVA space {g ∈ L 2 (µ) s.t. g = g I } → Illustrated on main and total effects, but very general! Tensors are used to get lower bounds as a sum of squared integrals → Chaos polynomials or more general tensors O. Roustant (EMSE) Sobol-Hoeffding decomposition: bounds and extremes 2018 December 11 Partial conclusions Lower bounds of a (convex comb. of) ANOVA term g I can be obtained by projection onto subspaces of its ANOVA space {g ∈ L 2 (µ) s.t. g = g I } → Illustrated on main and total effects, but very general! Tensors are used to get lower bounds as a sum of squared integrals → Chaos polynomials or more general tensors Integration by part modify lower bounds into derivative-based forms → Specific choices of subspaces remove weights for specific pdfs Oa (convex comb. of) ANOVA term g I can be obtained by projection onto subspaces of its ANOVA space {g ∈ L 2 (µ) s.t. g = g I } → Illustrated on main and total effects, but very general! Tensors are used to get lower bounds as a sum of squared integrals → Chaos polynomials or more general tensors Integration by part modify lower bounds into derivative-based forms → Specific choices of subspaces remove weights for specific pdfs

t) s .

 s As the terms of S.-H. decomposition are obtained by linear operation, inference properties of the stdf transfer to its ANOVA terms... of S.-H. decomposition are obtained by linear operation, inference properties of the stdf transfer to its ANOVA terms... Consider the usual assumptions for stdf inference, with corresponding valid sequences k = k (n). Then, for all I ⊆ {1, . . . , d},sup u I ∈[0,1] |I| | ˆ k ,n;I (u I ) -I (u I )| P -→ 0 . √ k ˆ k ,n;I (u I ) -I (u I ) d -→ Y ;I (u I )

  ... and hence to the tail dependographD I ( ˆ k ,n ) P -→ D I ( ) If D I ( ) > 0, then D I ( ˆ k ,n ) is asympt. normal with rate √ k If D I ( ) = 0, then D I ( ˆ k ,n ) is asympt.χ 2 type with rate k (The same is true for D tot I ) Inference (A piece of intuition about asymptotic distribution) ˆ k ,n;I (u I ) = I (u I ) + 1 √ k Y ,I (u I ) + . . . of intuition about asymptotic distribution) ˆ k ,n;I (u I ) = I (u I ) + 1 √ k Y ,I (u I ) + . . . If I ≡ 0, ˆ 2 k ,n;I (u I )dµ(u) D I ( ˆ k ,n ) of intuition about asymptotic distribution) ˆ k ,n;I (u I ) = I (u I ) + 1 √ k Y ,I (u I ) + . . . If I ≡ 0, ˆ 2 k ,n;I (u I )dµ(u) D I ( ˆ k ,n ) ;I (u I )dµ(u) D I ( ˆ k ,n ) maxima temperatures at 21 French cities during 1946 -2000.

Figure :

 : Figure: Estimated tail dependograph: complete, 30 largest values, 9 largest

Figure :

 : Figure: Estimated tail dependograph: complete, 30 largest values, 9 largest

Link between total Sobol indices and DGSM Theorem [Lamboni et al., 2013], [Roustant et al., 2014]

	Roustant (EMSE)	Sobol-Hoeffding decomposition: bounds and extremes	2018 December 11	15 / 50

Link between total Sobol indices and DGSM Theorem [Lamboni et al., 2013], [Roustant et al., 2014]

	Roustant (EMSE)	Sobol-Hoeffding decomposition: bounds and extremes	2018 December 11	16 / 50

Link between total Sobol indices and DGSM Theorem [Lamboni et al., 2013], [Roustant et al., 2014]

	EMSE)	Sobol-Hoeffding decomposition: bounds and extremes	2018 December 11	16 / 50

Link between total Sobol indices and DGSM Theorem [Lamboni et al., 2013], [Roustant et al., 2014]

	Roustant (EMSE)	Sobol-Hoeffding decomposition: bounds and extremes	2018 December 11	17 / 50

  This improves on the known lower bound which has the same form, with the smaller multiplicative constant 2m+1 (m+1) 2 .

O.

Roustant

(EMSE) Sobol-Hoeffding decomposition: bounds and extremes 2018 December 11 32 / 50

  This improves on the known lower bound which has the same form, with the smaller multiplicative constant 2m+1 (m+1) 2 .

	Results on the application Results on the application
	1.0 1.0 1.0 1.0 1.0
	0.8 0.8 0.8 0.8 0.8		Optimal bound Total Sobol index Optimal bound Total Sobol index
			Lower bound
	0.6 0.6 0.6 0.6 0.6
	0.4 0.4 0.4 0.4 0.4
	0.2 0.2 0.2 0.2 0.2
	0.0 0.0 0.0 0.0 0.0
	Q Q Q Q Q	Ks Ks Ks Ks Ks	Zv Zm Hd Cb Zv Zm Hd Cb Zv Zm Hd Cb Zv Zm Hd Cb Zv Zm Hd Cb	L L L L L	B B B B B
					2018 December 11 2018 December 11	32 / 50 33 / 50

N.B. Better bounds are obtained by adding orth. funct. of the form ψ 1 ψ j . O. Roustant (EMSE) Sobol-Hoeffding decomposition: bounds and extremes Figure: Results obtained with orth. 1st order pol. tensors ψ 1 , ψ 1 ψ 2 , . . . , ψ 1 ψ 8 O. Roustant (EMSE) Sobol-Hoeffding decomposition: bounds and extremes Figure: Results obtained with orth. 1st order pol. tensors ψ 1 , ψ 1 ψ 2 , . . . , ψ 1 ψ 8 O. Roustant (EMSE) Sobol-Hoeffding decomposition: bounds and extremes 2018 December 11

  gφ 1 φ j , . . .,

	or	h der =	∂g ∂x i	,	∂g ∂x j	Z j , ...

for centered function φ 1 , φ j , Z j . O. Roustant (EMSE)

Sobol-Hoeffding decomposition: bounds and extremes 2018 December 11

dµ(x) O. Roustant (EMSE) Sobol-Hoeffding decomposition: bounds and extremes 2018 December 11

O. Roustant (EMSE)Sobol-Hoeffding decomposition: bounds and extremes 2018 December 11

dµ(x) = C(µ i )ν i O. Roustant (EMSE)Sobol-Hoeffding decomposition: bounds and extremes 2018 December 11

dµ(x) = C(µ i )C(µ j )ν i,j O. Roustant (EMSE)Sobol-Hoeffding decomposition: bounds and extremes 2018 December 11

O. Roustant (EMSE)Sobol-Hoeffding decomposition: bounds and extremes 2018 December 11

Joint work with C. Mercadier O. Roustant (EMSE) Sobol-Hoeffding decomposition: bounds and extremes 2018 December 11

Acknowledgements

Part of this research was conducted within the frame of the Chair in Applied Mathematics OQUAIDO, gathering partners in technological research (BRGM, CEA, IFPEN, IRSN, Safran, Storengy) and academia (CNRS, Ecole Centrale de Lyon, Mines Saint-Etienne, University of Grenoble, University of Nice, University of Toulouse) around advanced methods for Computer Experiments. O. Roustant (EMSE) Sobol-Hoeffding decomposition: bounds and extremes 2018 December 11 49 / 50 Part V References -Thank you for your attention!

Asymptotic independence and tail dependograph

Let A, B a partition of {1, . . . , d} X A and X B are asymptotically independent

Asymptotic independence and extremal coefficients

The extremal coefficients θ I ( ) are defined by

Sobol-Hoeffding decomposition: bounds and extremes 2018 December 11

Asymptotic independence and extremal coefficients

The extremal coefficients θ I ( ) are defined by

Asymptotic independence and extremal coefficients

The extremal coefficients θ I ( ) are defined by

Equivalently θ I ( ) = (1 I ) , and in particular 1 Asymptotic dependence

Sobol-Hoeffding decomposition: bounds and extremes 2018 December 11

Illustration: Revealing asymptotic dependence for asymmetric models Consider a 4-dim. random vector X with standard Gumbel margins, and s.t.d.f. built as a mixture of independence and logistic:

Inference

As the terms of S.-H. decomposition are obtained by linear operation, inference properties of the stdf transfer to its ANOVA terms...

Consider the usual assumptions for stdf inference, with corresponding valid sequences k = k (n). Then, for all I ⊆ {1, . . . , d},

where Y ;I is some Gaussian process.

O.


Franck Cappello

Amina Guermouche

Marc Snir
email: snir@illinois.edu

On communication determinism in parallel HPC applications (Invited Paper)

Current fault tolerant protocols for high performance computing parallel applications have two major drawbacks: either they require to restart all processes even in the case of only a single process failure or they have a high performance overhead in fault free situation. As a consequence none of existing generic fault tolerant protocols matches needs of HPC applications and surprisingly, there is no fault tolerant protocol dedicated to them. One way to design better fault tolerant protocols for HPC applications is to explore and take advantage of their specific characteristics. In particular we suspect that most of them present some form of determinism in communication patterns. Communication determinism can play an important role in the design of new fault tolerant protocols by reducing their complexity. In this paper, we explore the communication determinism in 27 HPC parallel applications that are representative of production workloads in large scale centers. We show that most of these applications have deterministic or send-deterministic communication patterns.

I. INTRODUCTION

Exploring and designing efficient fault tolerant protocols for parallel applications has always been a major research topic of the HPC community [START_REF] Rao | Egida: An extensible toolkit for low-overhead fault-tolerance. Fault-Tolerant Computing[END_REF], [START_REF] Bouteiller | Mpich-v project: A multiprotocol automatic fault-tolerant mpi[END_REF], [START_REF] Gao | Application-transparent checkpoint/restart for mpi programs over infiniband[END_REF], [START_REF] Hursey | The design and implementation of checkpoint/restart process fault tolerance for Open MPI[END_REF]. Fault tolerant protocols are essential for rollback recovery and they are also useful for migration and preemptive scheduling.

The role of fault tolerant protocols is to capture information of a distributed execution, during fault free execution, in such a way that the execution could be restarted and eventually terminates successfully producing a correct result. Several decades of research have generated a large variety of fault tolerant protocols [START_REF] Elnozahy | A survey of rollback-recovery protocols in messagepassing systems[END_REF]. In classifications, they are often grouped in different families: a) coordinated checkpointing [START_REF] Chandy | Distributed snapshots: determining global states of distributed systems[END_REF], [START_REF] Koo | Checkpointing and rollbackrecovery for distributed systems[END_REF], b) uncoordinated checkpointing, possibly with message logging [START_REF] Alvisi | Message logging: Pessimistic, optimistic, causal, and optimal[END_REF], [START_REF] Johnson | Sender-based message logging[END_REF], [START_REF] Bouteiller | Reasons for a pessimistic or optimistic message logging protocol in mpi uncoordinated failure recovery[END_REF] and c) communication induced checkpoint Supported by INRIA-Illinois Joint Laboratory on PetaScale Computing. Supported in part by NSF grant 0833128 [START_REF] Quaglia | On the no-z-cycle property in distributed executions[END_REF]. These protocols do not fit the HPC context because they have a high overhead in fault free situation (message logging protocols), require to restart all processes in case of even a single process failure (coordinated checkpointing) or may lead to restart the execution at its beginning (uncoordinated checkpointing without message logging). These weaknesses are significant enough to motivate further research. One key differentiator between existing fault tolerant protocols is the level of determinism they assume on applications. Coordinated checkpointing and uncoordinated checkpointing without message logging were designed to work in the extreme situation where all processes actions can be non deterministic. These fault tolerant protocols compute checkpoints without considering any assumption on the execution following a restart after a failure. This execution after restart can be completely different and incomparable to the original execution before the occurrence of the fault. As a consequence these protocols need to restart all processes. Message logging protocols improve this aspect and do not need to restart all processes because they assume more determinism: the application code for every process should be piecewise deterministic [START_REF] Strom | Optimistic recovery in distributed systems[END_REF]. One way to model the execution of the piecewise deterministic process is to consider that its state is only influenced by deterministic local actions and message receptions. The later may be the only non-deterministic events and all non-deterministic actions are modeled as non-deterministic receptions. One property of the piecewise deterministic assumption is that for any given process across multiple executions (or re-executions), sends are guaranteed to be identical if non-deterministic actions are replayed identically. This property is used by message logging protocols to avoid restarting all processes 1 . Thus considering more application determinism helped to solve some issues of coordinated checkpointing protocols. However message logging protocols have other limitations: they have a high performance overhead in fault free situation because they need to copy all messages and store message copies and message receptions order of every process on stable storage.

In principle, if communications in parallel HPC applications were deterministic, the overhead of fault tolerant protocols could be drastically reduced since most of it is due to guarantee consistent restarts and correct termination of distributed executions in the presence of non-determinism. The initial results in [START_REF] Cappello | Revisiting fault tolerant protocols for hpc applications[END_REF] indicates this to be the case. However, in many HPC MPI applications, programmers relax execution determinism to increase concurrency, reduce coordination overhead and improve scheduling. Resulting codes ofter use minimal synchronization and allow a high degree of asynchrony. The objective of this paper is to analyze a significant panel of HPC applications and explore the existence of determinism sources that could be exploited later to design better fault tolerance protocols 2 . More precisely we want to explore the existence of determinism in communication patterns since HPC applications rarely feature non determinism in their computation sections 3 .

Next sections introduce two forms of communication determinism that we call "Communication-determinism" and "Send-determinism". We suspect that these forms of communication determinism is dominant in HPC applications. The rest of the paper will demonstrate this hypothesis. Before introducing the analyzed applications, we present in section III formal definitions of Communication and Send-Determinism. Then we introduce, in section IV, 27 application codes and benchmarks considered in this study. Section V gives the most significant send-deterministic and non-deterministic communication patterns found in the 27 applications codes. Section VI presents the result of this analysis.

II. COMMUNICATION DETERMINISM AND SEND-DETERMINISM IN MPI HPC APPLICATIONS

Many MPI HPC applications follow a SPMD programming style where computational algorithms are implemented by separate dependent communications and 1 to ensure a compatible re-execution, all restarting processes receptions are controlled to make sure that same emissions are replayed until reaching the point where the fault occur or the last communication before the end of the application. 2 Proposing new fault tolerant protocol is out-of-scope of this paper. 3 while an application may use random number generation, it generally set deterministically the seed at the beginning of the execution leading to the generation of a deterministic sequence of numbers. computations sections, although communications and computations may overlap during the execution. Typically, MPI HPC applications use point to point and collective communications. We assume in this paper that collective communications are implemented atop point-to-point communications and the implementation is determinisitic. The last condition is satisfied by the MPI implementations we are familiar with. Therefore, w.l.o.g., we can restrict our attention to point-to-point communication. We ignore in this paper MPI features that are rarely used, such as MPI_Cancel, or features that can be handled by extending the arguments in the paper, such as I/O or one-sided communication.

For sake of clarity and precision, we present relevant point-to-point communication functions of MPI. In MPI, point to point communications use mainly two functions: MPI_Send() for send and MPI_Recv() for reception. int MPI_Send( void *buf, int count, MPI_Datatype datatype, int dest, int tag, MPI_Comm comm) sends the message specified by the buffer address buf and the arguments datatype and count to the process dest. The message has the tag tag. The communicator argument comm specifies a communication domain that is both a name space (dest is interpreted relative to comm) and a separate communication plane (communications using different communicators do not interfere with each other). int MPI_Recv(void *buf, int count, MPI_Datatype datatype,int source, int tag, MPI_Comm comm, MPI_Status *status) is the function used to perform receives. The process receives a messages with the tag tag from the sender source. The received data is placed in the buffer indicated by buf, count and datatype. The argument status returns information on the received message.

The source argument and the tag argument can either or both have special "don't care" values that match any source and/or any tag. Messages are ordered; a receive will receive the first sent message from the specified source that has a matching tag; if no source is specified, it will receive the first matching message from some source.

These are blocking calls: The progression of the calling code is stopped until the communication is completed 4 . These two functions have non-blocking variants: MPI_Isend() and MPI_Irecv(). These functions return before the completion of the communication, in order to allow its overlapping with computation; they have an extra request argument that is used to return a handle to the ongoing communication. The programmer can use other functions to test the status of the communication. MPI_Test (mpi_Request *request, int *flag, MPI_Status *status) tests for the completion of the communication indicated by request; if the communication is complete, then information on the received is returned in status. MPI_Wait (mpi_Request *request, int *flag, MPI_Status *status) blocks until the communication is complete. In addition, functions are provided to wait for the completion of any or all requests in a list. MPI applications can exhibit nondeterminism in multiple ways. In particular, the MPI process can be multithreaded, with a nondeterministic concurrent execution; the timing of a message arrival can impact execution logic, when MPI_Test is used; and the relative arrival time of messages sent by different sources to the same destination can affect the matching of sends to receives, when the receives specify a "don't care" source. However, HPC codes tend to eschew non-determinism, or use it in very limited form. A typical form of nondeterminism is code where multiple concurrent incoming messages will be handled in the (non-deterministic) order they arrive. (This can be done by posting multiple non blocking receives then using an MPI_Waitany call to complete one of the posted receives.) The execution logic of the receiving process will depend on the arrival order of the messages. However, it is typically the case that the outcome of the computation will node depend on the arrival order: The execution of the process will consist of multiple sections, where there is non-determinism within each section, but the state at the end of each section is deterministic. Furthermore, non-determinism "does not propagate": the order and content of sent messages is not affected by receive order. We call such execution send-deterministic.

Verifying that communication patterns in MPI HPC applications are deterministic or "send-deterministic" is important since this will allow the design of new fault tolerant protocols. The objective here is not to define and demonstrate the correctness and properties of some new protocols but to introduce their main principles: Deterministic or Send-deterministic communication patterns means that message sends for every process are identical across executions, from a given set of input parameters. Thus instead of logging all messages, as needed by message logging protocols, to reconstruct the state of failed processes, these messages could be reconstructed from the re-execution of their senders. Since deterministic and send-deterministic communication patterns guarantee that same messages will be resent by the re-execution of their senders, the state of a process can be re-built from its re-execution and the re-execution of the processes that sends it a message. Thus, in principle, messages do not need to be logged to lead a process to a point where it execution is coherent with the other processes. Note, however, that some precautions are needed to avoid a domino effect (forcing the whole execution to restart from the beginning). In [START_REF] Cappello | Revisiting fault tolerant protocols for hpc applications[END_REF], we present two new fault tolerant protocols based on these principles.

Before checking the existence of determinism and send-determinism in MPI HPC applications, let us define more formally these notions.

III. FORMAL DEFINITION OF DETERMINISM TYPES IN PARALLEL COMPUTATION

Characterizing MPI HPC applications with regards to the determinism of their communication patterns requires a rigorous definition of the introduced notions: determinism and send-determinism. We first give the formal definition of communication determinism and then we give the formal definition of send-determinism. These definitions will be used as the basis for the analysis of the communication patterns in the MPI HPC 27 applications considered in this work. In the following, we focus only on the determinism in communication operations and do not address determinism of the entire computation.

We model a parallel computation as consisting of a finite set of processes and a finite set of channels connecting any (ordered) pair of processes. Channels are assumed to have infinite buffers, to be error-free, and to deliver messages in the order sent.

We assume we have a set P of n processes. Each process p ∈ P is defined by a set of states C p , and initial state c p 0 , a set of events V p and a partial transition function

F p : C p ×V p -→ C p ; if p = q then V p ∩V q = ∅. A state c ∈ C p is final if there is no valid transition from state c.
A channel is modeled as an unbounded queue. A send event send(p, q, m) ∈ V p adds message m to the tail of channel (p, q); a receive event receive(q, p, m) ∈ V q deletes the message m from the head of channel (p, q).

A computation consists of a sequence E = e 1 , e 2 , . . . of events. A system configuration consists of a vector of states C = {c i , 1 ≤ i ≤ n, , c i,j , 1 ≤ i, j ≤ n} where c i is the state of process i and c i,j is the state of channel (i, j). The sequence of configurations resulting from execution E is defined as follows.

The initial configuration of the system has each process in its initial state and each channel is empty.

Let C i = {c i p , c i p,q } be the configuration at step i. Assume that e i ∈ V p . Then,

c i+1 q = c i q if q = p F p (c i p , e i ) if q = p
In addition, if e i = send(p, q, m) then message m is appended to channel (p, q) and if e i = receive(q, p, m) then message m is deleted from channel (p, q).

The execution is valid if

• F p (c i p , e i
) is defined at each step • If the event receive(q, p, m) occurs then message m is at the head of channel (p, q). • The last state of each process is final. We denote by E the set of valid executions. For each execution E we denote by E|p the subsequence of E consisting of events in V p and by C(E)|p the sequence of states of process p after each event in E|p. Then Definition 1: A parallel computation is deterministic if, for each p, C(E)|p is the same for any E ∈ E .

The computation is communication-deterministic if, for each p, E|p contains the same subsequence of communication events (sends and receives) for any E ∈ E .

The computation is send-deterministic if, for each p, E|p contains the same subsequence of send events for any E ∈ E .

This abstract model can be related to MPI programs as follows:

Program transitions are modeled by noncommunication events; if the program is determinisitc, then each state is associated with a unique valid event.

The execution of an MPI send is a transition caused by the event that appends the sent message to the channel connecting the sender to the receiver. The transition function can be defined so that this event is the only valid event for a given process state.

The arrival of a message is modeled by a receive event. One can think of this event as moving a message from the communication channel into a process buffer -the process state is changed to represent the new buffer content. The actual execution of the MPI receive operation is an internal program transition.

Before analyzing the nature of communication determinism in HPC applications considered in this study, we introduce the methodology as well as the applications.

IV. METHODOLOGY AND ANALYZED APPLICATIONS

The methodology followed in this study consisted in a static analysis of the application codes. While automatic analysis would be needed for analyzing a large number of applications, we decided to use a manual analysis. To analyze the nature of the determinism in these applications, based on the formal definitions of communication-deterministic and send-deterministic program, we decomposed each of them as a sequence of computation and communication patterns (computation patterns are deterministic). For every application, we analyzed the deterministic nature of all communication patterns and characterize the communication determinism of the application as following: 1) an application is communication-deterministic if and only if it contains only communication-deterministic communication patterns, 2) an application is send-deterministic if and only if it contains at least one send-deterministic communication pattern and no non-deterministic communication pattern, 3) an application is non-deterministic if and only if it contains at least one non-deterministic communication pattern. We considered 27 applications, benchmarks or kernels that we consider representative of a large fraction of applications run on production platforms. They are well recognized benchmarks, kernels used by many applications, codes of general science applications or they will serve for the evaluation of the next generation supercomputers (of 10+ Petaflops). Several codes are mixing Fortran 90 and MPI. We assume in the following that Fortran 90 is not a source of non determinism in the order of message reception.

The first set of applications belongs to the NERSC-6 benchmarks [START_REF] Antypas | Nersc-6 workload analysis and benchmark selection process[END_REF]. The NERSC center serves 3,000 scientists throughout the USA working on problems in combustion, climate modeling, fusion energy, materials science, physics, chemistry and computational biology. The suite consists of seven applications, spanning a wide range of science disciplines, algorithms, concurrencies and scaling methods. In this study, we considered the following applications: CAM, GTC, IMPACT, MAE-STRO, PARATEC. CAM is the dominant computational code of the fully coupled CCSM3 climate model and is considered as one of the most demanding codes of the climate workload. CAM calculations relies on spectral (FFTs) and structured grids Methods. It uses many collective communications as well as blocking and non blocking point-to-point communications. GTC is used for fusion energy research, like for studying neoclassical and turbulent transport in tokamaks and for investigating hot particle physics. It uses particle and structured Grids methods. Its code includes collective communications, MPI_SendRecv() and point-topoint communication with the use of wildcards. IMPACT is used for accelerator sciences and represents typical beam dynamics simulation workload. IMPACT code was used for the design of many linear accelerators. It relies on spectral (FFTs), particle and structured grids methods. IMPACT mainly uses point-to-point communications in blocking and non-blocking versions. MAESTRO is an astrophysics code used to simulate the long time evolution leading up to a supernova explosion. It relies on explicit and implicit solvers and uses sparse linear algebra, structured and unstructured grids methods. MAESTRO uses BoxLib, a library of Fortran90 modules that facilitates development of block-structured finite difference algorithms. It uses collective communications as well as non-blocking point-to-point communications. PARATEC (PARAllel Total Energy Code) performs quantum mechanical total energy calculations. It represents typical material science workloads run at NERSC. It uses dense linear algebra, spectral methods (FFTs) and structured grids methods. PARATEC uses collective communications, blocking and non-blocking point-to-point communications.

The second set of applications is included in the Sequoia Benchmarks [15]. This suite of codes is an evolution of the ASCI Purple benchmarks. Several codes have been rewritten and added to the suite. The Sequoia supercomputer will run codes for predicting stockpile performance. We analyzed the five application codes from the suite: Sphot, UMT, AMG, IRS and lammps. We also analyzed IOR, an I/O benchmark using MPI calls. SPhot implements Monte Carlo photon transport on a small, 2D structured mesh. Its algorithm is embarrassingly parallel. UMT performs 3D, deterministic, multi-group, photon transport on an unstructured mesh. It uses a transport algorithm solving the first-order form of the time-dependent Boltzmann transport equation. AMG uses algebraic multi-grid algorithms to solve large, sparse linear systems derived from implementing physics simulations on unstructured or block structured meshes. IRS iteratively solves a 3D radiation diffusion equation set on a block-structured mesh. It represents an important physics package representative of computation and memory access patterns used in several production physics packages. lammps can simulate a wide variety of different particle systems. lammps uses nearest neighbors and reduction operations. SPhot, UMT, AMG IRS and lammps use collective communications, blocking and non-blocking point-to-point communications. IOR is used for testing the performance of parallel filesystems from various HPC access patterns. We will discuss the nature of its specific communication pattern in next section.

The third set of applications concerns QCD calculation and is part of the USQCD suite. QCD simulations play an important role to better understand the fundamental laws of physics. USQCD suite contains software enabling high performance lattice quantum chromodynamics computations across a variety of architectures. Lattice QCD calculations allow understanding the results of particle and nuclear physics experiments in terms of QCD, the theory of quarks and gluons.The application packages used by USQCD are: the Chroma Code, the Columbia Physics System (CPS), FermiQCD and the MIMD Lattice Collaboration (MILC) Code. We analyzed CPS and MILC. CPS is a large set of lattice QCD codes that are essential for QCD applications. CPS uses collective MPI communications, blocking and non-blocking point-to-point communications. The MILC Code is a set of codes developed by the MIMD Lattice Computation (MILC) collaboration for doing simulations of four dimensional lattice gauge theory on parallel machines. MILC uses collective MPI communications, blocking and non-blocking point-to-point communications with wildcards.

SPECFEM3D [START_REF] Carrington | High-frequency simulations of global seismic wave propagation using SPECFEM3D_GLOBE on 62 thousand processor cores[END_REF] is a software package that simulates southern California seismic wave propagation. SPECFEM3D won the Gordon Bell Award for Best Performance at the ACM/IEEE SuperComputing'2003. The package contains several software including a mesher and a solver. SPECFEM3D uses a new method introduced by its authors for the calculation of synthetic seismograms in 3-D earth models: the Spectral Element Method (SEM). SPECFEM3D uses collective communications and MPI_SendRecv() where source and destination and tag parameters are deterministic.

We analyzed another geophysics application, RAY2MESH [START_REF] Grunberg | Seismic ray-tracing and earth mesh modeling on various parallel architectures[END_REF], because it is available in two versions: a SPMD version and a Master-Worker version. RAY2MESH is a high performance set of software tools for seismic tomography aiming to compute a realistic geophysic Earth model. RAY2MESH uses massive quantities of seismic waves records to ray-trace waves propagation inside a global Earth mesh with thin cells. The two versions use different MPI communication patterns. The Master-Worker version uses a pull model of scheduling.

In terms of small benchmark codes, we analyzed the NAS parallel Benchmarks NPB 3.3 [START_REF] Bailey | The nas parallel benchmarks-summary and preliminary results[END_REF]. These programs are derived from computational fluid dynamics (CFD) applications to help evaluate the performance of parallel supercomputers. The benchmark suite consist of five kernels and three pseudo-applications. In this study, we analyzed LU, MG, BT, SP, CG, EP; FT and DT. These programs use a large variety of MPI function calls including collective communications, blocking and nonblocking point-to-point communications.

In addition to applications run in production center, we analyzed several kernels used in many applications: 1) SUMMA, the matrix product kernel of the ScaLAPACK, 2) Linpack, the kernel used for the top500, 3) a Jaccobi kernel, 4) a nbody kernel.

V. EXAMPLES OF COMMUNICATION PATTERNS

In this section, we present the most common communication patterns found on analyzed applications. We focus only on patterns that are not communicationdeterministic. This includes send-determinitic and nonderterministic communication patterns.

A. Send-deterministic communication patterns

One of the most common send-deterministic patterns uses a sequence of MPI_Irecv and MPI_Isend followed by an MPI_Waitall. This pattern is send-deterministic because an Irecv can complete before another one that was posted before. But the fact that the source is known ensures that every data is received where it should be (there is no risk of confusion), for example: Note that there may be some variations, for example an MPI_Send instead of an MPI_Isend, or n*MPI_Wait instead of MPI_Waitall (where n is the number of Irecv). Another pattern, similar to the former, contains MPI_WaitANY instead of MPI_WaitAll or MPI_WAIT. It is found in the Sequoia lammp benchmark. MPI_Waitany waits for completion of one over several requests.

VI. RESULTS ANALYSIS

In this section, we present the result of our communication determinism analysis over 27 HPC applications. Table 1 presents, for every application, the diversity of used communication patterns and the class to which every application belongs.

From these results, we can draw some conclusions. First very few applications have non-deterministic communication patterns. This is an important result since this makes relevant the research for new fault tolerant protocols exploiting the determinism in HPC applications communication patterns. Second a significant portion of the applications have send-deterministic communication patterns where programers introduce more asynchrony in their applications. We believe that this kind of applications will become prevalent in the near future because the multi-core trend pushes architectures and algorithms for more parallelism and less synchrony. As a consequence, it may be irrelevant to design fault tolerant protocols uniquely for applications featuring deterministic communication patterns. Third, some applications, like Sequoia-AMG present a high degree of communication determinism. However because they have some non-deterministic communication patterns, fault tolerant protocols designed for communication-deterministic and send-deterministic communication patterns will not work for these applications. If the number of non-deterministic communication patterns is small, like for Sequoia-AMG, then it may be worth to reprogram these communication patterns to turn then communication-deterministic or send-deterministic. Fourth, a large portion of the applications use collective communications. Some applications are using them extensively. In this study, we have considered their internal implementation as communication deterministic. If there is a need to inject more asynchrony and relax determinism for performance purpose, our advise if to use send-deterministic implementations and avoid non-deterministic ones. This will make sure that specific fault tolerant protocols designed for communication-deterministic and send-deterministic communication patterns stay relevant for a large fraction of applications.

VII. CONCLUSION

In this paper we explored the determinism nature of communication patterns in a large variety of HPC applications representative of production workloads. This research was motivated by the fact that better fault tolerant protocols could be designed if some form of determinism can be exploited in communication patterns. We analyzed 27 applications from different benchmark and application suites. As we expected, most of HPC applications feature a high degree of determinism in their communication patterns: they are either communication deterministic or send-deterministic. Very few applications have non-deterministic communication patterns. This result is significant since it makes relevant the research for new fault tolerant protocols leveraging the determinism in HPC application communication patterns. Apart from the design of new fault tolerant protocols, there are several perspectives of this research: how to turn non-deterministic communication patterns into more deterministic ones while maintaining performance? New developments in MPI allow several threads to make communication calls concurrently. Applications using this possibility will have less deterministic communication patterns compared to communication-deterministic and send-deterministic communication patterns. However, the outcome of any communication pattern may still be deterministic regardless of the order of communication operations inside each pattern. We will analyze and characterize this form of determinism and explore new fault tolerance protocols from it.

TABLE I :

 I Communication determinism in parallel applications. 'x' is used to denote the presence of the pattern in the application. Numbers give the total occurrence of the patterns in the application

Note that there are several versions of these functions with different semantic of the termination.

where tag [i] is the communication tag. In this application, the tag is different for every sender. This tag avoids confusions since every communication has a different tag. The tag can also be ANY_TAG, as in Sequoia IOR (with an MPI_Recv): This pattern counts the number of processes running on the same host (hosting the process with rank 0). The reception order has no impact on the final count. Thus we consider this pattern as send-deterministic. Note however, that there is another condition to satisfy in this specific pattern: the number of processes hosted with process rank 0 should be deterministic for a given set of input parameters and execution parameters.

A subtle example of send-deterministic pattern using ANY_SOURCE and ANY_TAG is described below:

.., ANY_SOURCE, ANY_TAG,...); } else MPI_Send(T[rank], ..., 0, 0, ...); MPI_Barrier(...); sort(T);

In that case, there is no way to ensure that the content of array T after the MPI_Barrier() is the same across several executions. Note however that in this application the reception order has no importance if we extend the notion of communication pattern to include operations on receive buffers preceding their utilization to influence the rest of the system. This is because: 1) the result is never sent and 2) just after the loop, T is sorted and will always lead to the same result (for a given set of input parameters). So the subtleness here is on the definition of the "communication pattern".

B. Non-deterministic communication patterns

An example of non-deterministic patterns is shown below (from Sequoia AMG benchmark): The non-determinism is due to the MPI_IProbe with the ANY_SOURCE wilcard. MPI_IProbe checks if there is an incoming messages from any source with tag1. If so, an MPI_Recv is posted. So, receive order is not deterministic. And since the following send depends on the receive, sends may also occur in a different order from one execution to another. This does not respect one of the send-determinism conditions (for one sender, messages are always sent in the same order).

Another example of non-deterministic communication pattern is used in the Master/Worker implementation of the Ray2mesh application. The Master-Worker version of Ray2Mesh implements a first serve first come scheduling where workers pull tasks from the master. A faster worker will get more tasks. So the physical time is used in this program to control the distribution of tasks to the worker. This is representative of many Master-Worker application implementations.

Once the root process receives a request from a worker, it sends a task to the worker. The non-determinism is due to MPI_Recv(..., ANY_SOURCE, ...) and the send that depends on it. This program is not send-deterministic because the send order fundamentally depends on the non-deterministic receive order.