Chen Chen

Adrien Ugon

Xun Zhang

Amara Amara

Patrick Garda

Jean-Gabriel Ganascia

Amina Kotti

Carole Philippe
email: carole.philippe@aphp.fr

Andrea Pinna

Cross Entropy-based Automatic Thresholds Setting-Up Method for Sleep Staging System

Keywords: Cross Entropy, Thresholds, Personalization, Automatic Sleep Staging Systems

come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

Sleep disorders are affecting more and more people which deteriorate quality of life and become a significant cause of morbidity and mortality [START_REF] Hossain | The prevalence, cost implications, and management of sleep disorders: an overview[END_REF]. In diagnosis and treatment of sleep disorders, sleep staging is considered to be the fundamental step. Clinical sleep staging is based on the visual analysis of an overnight polysomnography (PSG) recordings including electroencephalogram (EEG), electrooculogram (EOG) and electromyogram (EMG) by a physician under the guidance of American Academy of Sleep Medicine (AASM) manual [START_REF] Of Sleep Medicine | The AASM manual for the scoring of sleep and associated events: rules, terminology and technical specifications[END_REF]. It is a time-consuming and labor intensive task. Inter-rater variability also exists [START_REF] Danker-Hopfe | Interrater reliability for sleep scoring according to the rechtschaffen & kales and the new aasm standard[END_REF].

To release the burden of physicians, different sleep staging systems have been proposed in the past decades [START_REF] Schaltenbrand | Neural network model: application to automatic analysis of human sleep[END_REF], [START_REF] Liang | A fuzzy inference system for sleep staging[END_REF], [START_REF] Liang | A rule-based automatic sleep staging method[END_REF], [START_REF] Ugon | Fusion symbolique et données polysomnographiques[END_REF], [START_REF] Chen | Symbolic fusion: A novel decision support algorithm for sleep staging application[END_REF]. Among existing sleep staging systems, thresholds are widely applied in transforming digital parameters into linguistic or symbolic features to model inference process under the guidance of medical knowledge [START_REF] Liang | A fuzzy inference system for sleep staging[END_REF], [START_REF] Liang | A rule-based automatic sleep staging method[END_REF], [START_REF] Ugon | Fusion symbolique et données polysomnographiques[END_REF], [START_REF] Chen | Symbolic fusion: A novel decision support algorithm for sleep staging application[END_REF].

In clinical practice, boundaries of linguistic or symbolic features are very flexible. Physicians may adjust the boundaries for each linguistic or symbolic feature according to their experience and patient information. While, in sleep staging systems, thresholds are used to determine the boundaries in differentiating different linguistic or symbolic features. However, as far as we know, there are no fully satisfying automatic setting-up thresholds methods in existing sleep staging systems. Manually predefined values of thresholds have been widely applied in [START_REF] Liang | A fuzzy inference system for sleep staging[END_REF], [START_REF] Liang | A rule-based automatic sleep staging method[END_REF], [START_REF] Ugon | Fusion symbolique et données polysomnographiques[END_REF], [START_REF] Chen | Symbolic fusion: A novel decision support algorithm for sleep staging application[END_REF] due to the following reasons: 1) To build a mathematical model or a threshold function to differentiate boundary requires a set of data with sufficient quantity and adequate quality; 2) There is lack of uniformity between subjects and thresholds variability exists [START_REF] Maggard | Automation of sleep staging[END_REF].

In this paper, we propose an automatic thresholds settingup method by using Cross Entropy [START_REF] Rubinstein | Optimization of computer simulation models with rare events[END_REF] to overcome the limitations of existing thresholds setting-up problems. With the stable and rapid searching ability in selecting the optimal thresholds combination, it can also be adapted to different person to provide personalized thresholds which is suitable to high inter-subject variability sleep staging analysis.

This paper is organized as follows. In the next section, a Symbolic Fusion-based sleep staging system is introduced with discussion on the limitations of thresholds setting-up problem. Then, an automatic thresholds setting-up method is proposed based on Cross Entropy. Followed by results and comparison with manual approaches in section four. Lastly, conclusion and further work are provided.

II. SYMBOLIC FUSION-BASED SLEEP STAGING SYSTEM

In [START_REF] Ugon | Fusion symbolique et données polysomnographiques[END_REF], [START_REF] Chen | Symbolic fusion: A novel decision support algorithm for sleep staging application[END_REF], a sleep staging system is proposed using symbolic fusion which consists three levels: data fusion, feature fusion and decision fusion as shown in Figure 1. In data fusion, digital parameters are extracted from PSG recordings to maximize the useful information and to minimize noise and artifacts. In feature fusion, digital parameters are transformed into feature parameters. It simplifies the interpretation of digital parameters, and also performs normalization, reduction and matching of digital parameters. In decision fusion, inference method is used to fulfill sleep staging on the basis of feature parameters. In [START_REF] Ugon | Fusion symbolique et données polysomnographiques[END_REF], [START_REF] Chen | Symbolic fusion: A novel decision support algorithm for sleep staging application[END_REF], thresholds are used for translating nine digital parameters into 24 features parameters as shown in Table I. However, to define these thresholds, manual interpretation was performed. It is a time-consuming process; individual variability also exists.

III. AUTOMATIC THRESHOLDS SETTING-UP METHOD

To overcome the limitations in [START_REF] Ugon | Fusion symbolique et données polysomnographiques[END_REF], [START_REF] Chen | Symbolic fusion: A novel decision support algorithm for sleep staging application[END_REF], Cross Entropybased Automatic Thresholds Setting-Up (CEATS) is proposed as shown in Figure 2. It is dedicated to (1) provide optimal thresholds for sleep staging system (2) release the burden of manual interpretation and (3) allow individual flexibility. Cross Entropy is introduced firstly, followed by details of CEATS.

A. Introduction on Cross Entropy

Cross Entropy (CE) was initially proposed to estimate probabilities of rare events for complex stochastic networks by Rubinstein [START_REF] Rubinstein | Optimization of computer simulation models with rare events[END_REF] in 1997. It has been extended to solve combinatorial optimization problems in 1999 [START_REF] Rubinstein | The cross-entropy method for combinatorial and continuous optimization[END_REF], which turned out to be an effective method.

CE is briefly introduced as follows. In solving the combinatorial optimization problem, a maximization problem can be described as shown in Equation 1.

γ * = max x∈χ [S (x)] (1)
γ * represents the maximal value of S on the domain space χ. To proceed with CE, f (•; v) is defined as a family of Probability Density Functions (PDFs) on χ, with respect to some base measure v. Then γ * can be estimated by (γ) defined in Equation 2. Setting-Up Method

(γ) = P u (S (X) ≥ γ) = E u I {S(X)≥γ} ( 2
)
where X is a random vector generated by PDFs with parameter v in f (x, v). P u is the probability of the state {S(X) ≥ γ}, E u is the corresponding expectation operator and I(•) is the indicator function, i.e., I {S(X)≥γ} = 1 only if S(X) ≥ γ, otherwise, it equals to zero.

Based on the important sampling: take N random samples X = (X 1 , X 2 , ..., X N ) from an important sampling density g on χ, the unbiased estimator (γ) of (γ) can be defined as shown in Equation 3.

(γ) = 1 N N i=1 I {S(Xi)≥γ} f (X i ; v) g(X i ) = (γ) (3)
From Equation 3, only one sample suffices to estimate (γ) since it is true for all i. While it is difficult to directly compute g because this g depends on the unknown parameter . Moreover, it is convenient to choose a g in the family of densities f (•; v) . The CE method solves this efficiently by finding the minimal Kullback-Leibler distance as shown in Equation 4 which defines the distance between g and f (•; v).

D(g, f ) = g(x)lng(x)dx -g(x)lnf (x; v)dx (4)
In solving combinatorial optimization problems, CE creates a sequence of f (•; v 1 ), f (•; v 2 ), ... of PDFs that are driven in the direction of the theoretically optimal density f (•; v * ). f (•; v * ) corresponds to the degenerate density of the optimal solution. In each iteration, it generates a set of samples and the elite samples (in terms of solution quality) would be selected to update the parameters of the PDF f (x; v) parameterized by v. Since the elite samples are selected in each iteration, γ would be improved and can converge quickly to the optimal solution γ * .

B. Cross Entropy-based Automatic Thresholds Setting-Up Method

We propose CE-based method to set up the personalized thresholds for a Symbolic Fusion-based sleep staging system [START_REF] Ugon | Fusion symbolique et données polysomnographiques[END_REF], [START_REF] Chen | Symbolic fusion: A novel decision support algorithm for sleep staging application[END_REF]. As shown in the Figure 2, a training set is initially selected. Clinical sleep staging is performed by an expert to classify the training set into different sleep stages ("Expert Analysis"). In parallel, automatic analysis is processed ("Data Fusion", "Feature Fusion" and "Decision Fusion"). In order to assess impact of different thresholds, "Assessment" of sleep staging is adopted by comparing results of automatic and expert analysis. With assessment results, CE can generate new thresholds combinations for the next loop. This process is repeated until optimal thresholds combination is found, i.e. until the "Terminate Condition" is satisfied.

To realize automatic thresholds selection for this sleep staging system, details of CEATS procedure is described below.

• Step 1. Initialization: Define a specified mechanism to generate PDFs. Defined the sample size N, elite sampling rate ρ and the max iteration number G.

To search for the optimal thresholds combinations for the sleep staging system, normal distribution mechanism is performed to generate PDFs for each threshold which can be presented as N (µ, σ 2 ). µ is the mean of the distribution and σ is the standard deviation. Initial value of µ and σ is set according to the distribution of each digital parameters as shown in Table I. Different sample sizes and elite sampling rates are analyzed in the next section. The max iteration is 50 which takes the consideration of computational time. IV. RESULTS In this section, we present the simulation results obtained by using CEATS and compare them with results obtained with Manually Thresholds Setting-Up (MTS).

A. Subjects, Recordings and Preprocessing

Overnight PSG signals were recorded from 14 subjects (9 males and 5 females) ranging from 22 to 65 years old (µ = 42, σ = 15.8) in La Pitié-Salpêtrière hospital (AP-HP), in Paris. AHI (average number of apneas and hypopneas per hour of sleep) ranges from 0 to 35.3 (µ = 13.7, σ = 19). PSG recordings were segmented into 30s epoch and manually scored by experts into five different stages : W, N1, N2, and N3 and R according to AASM manual. By using Symbolic Fusion, nine digital parameters have been extracted from raw PSG signals as shown in [START_REF] Ugon | Fusion symbolique et données polysomnographiques[END_REF], [START_REF] Chen | Symbolic fusion: A novel decision support algorithm for sleep staging application[END_REF].

B. Evaluation of Cross Entropy Parameters

Sample sizes of 100, 500, 1000 and 10000 with elite sampling rates 0.1, 0.3, 0.5 and 0.7 were set to evaluate the performance of CEATS.

Figure 3 shows F-Measure dependence on iteration number for different sample sizes of stage W classification for one patient using CEATS. Due to the stochastic ability of Cross Entropy, F-Measure in Figure 3 are mean values of 20 independent runs. Large sample size potentially increases the sample diversity and may provide fast convergence speed. However, it also increases the computational time. For the same sample size, the increase of elite sampling rate slows down the convergence speed. According to the simulation results, sample size of 500 with elite sampling rate of 0.1 provides optimum between F-Measure and computational complexity. To evaluate impact of training set on CEATS, four values were used to balance the size of training set vs evaluation set : 5% vs 95%, 10% vs 90%, 15% vs 85% and 20% vs 80%. CEATS selects the optimal thresholds according to different training sets and pass these thresholds to the corresponding evaluation sets. Training set is randomly selected.

Figure 4 shows F-Measure of evaluation sets that use different training sets in CEATS. It shows that with reduction of training set size used in CEATS, it almost has no impact on the F-Measure for stage W, N2, N3 and R; only very slight difference for stage N1. In order to reduce the size of training sets, 5% is adopted in this paper. F-Measure is used to evaluate the impact of different thresholds on the sleep staging system. CEATS can provide personalized thresholds automatically for sleep staging system; while MTS needs manual efforts in setting-up thresholds. Comparison between CEATS and MTS [START_REF] Ugon | Fusion symbolique et données polysomnographiques[END_REF], [START_REF] Chen | Symbolic fusion: A novel decision support algorithm for sleep staging application[END_REF] is shown in Table II II, F-Measure for all the stages is increased by using CE-based method. For these 14 patients, average 40.6%, 13.1%, 32.3%, 74.1% and 51.3% F-Measure are improved for stage W, N1, N3 and R, respectively. CEATS proved to be an efficient method which can automatically set up personalized thresholds for each subject. It can be applied to realize personalized sleep staging. To our knowledge, this is the first study to explore CE in thresholds setting-up problems for the sleep staging applications and it can also be extended to other domains.

V. CONCLUSION

In this paper, a CE-based Thresholds Setting-Up method is proposed for automatic sleep staging system. It automatically provides optimal thresholds for releasing the burden of manual interpretation of thresholds and considering individual variability. Instead of building a mathematical model or a thresholds function, CE is a way to search for the optimal thresholds combination among the possible thresholds combinations space. With robustness and rapid searching ability, CE demonstrates to be an effective method in solving combinatorial thresholds setting-up problems. For the further work, proposed method can be extended to realize complete personalized sleep disorders analysis and plan to be integrated in an embedded system to realize remote sleep monitoring.

Figure 1 :

 1 Figure 1: Symbolic Fusion-based Sleep Staging System

Figure 2 :

 2 Figure 2: Cross Entropy-based Automatic Thresholds Setting-Up Method

•

  Step 2. Sampling: According to normal distribution density function N (µ, σ 2 ), a set of samples are generated. Each sample is composed of 15 different thresholds. • Step 3. Fitness function evaluation: F-Measure is used as the fitness function to evaluate the impact of different thresholds on the sleep staging system. It balances both precision and recall as shown in Equation 5. Precision equals to (Tp/(Tp + Fp)) and Recall equals to (Tp/(Tp + Fn)). The parameters TP, FN, and FP are respectively True Positives, False Negatives, and False Positives, which are used to quantify the quality of a classification. F -M easure = 2 * P recision * Recall P recision + Recall (5)

Figure 3 :

 3 Figure 3: Cross Entropy with Different Sample Sizes and Elite Sampling Rates

Figure 4 :

 4 Figure 4: Radar Chart of Different Training Sets

  VI. ACKNOWLEDGMENT This work was financially supported by China Scholarship Council (CSC) and under collaboration of Laboratoire d'Informatique de Paris 6 and Institut Supérieur d'Electronique de Paris. We would also like to thank Dr C. Philippe and Dr A. Kotti for the collaboration and for providing the database and help in the field of medicine.

Table I :

 I Thresholds used in Symbolic Fusion-based Sleep Staging System

	Digital Parameters	Features Parameters	Threshold
	EEG LowWaveEnergy	High -Middle -Low	2
	EEG SleepSpindles	Confident -Not Confident	1
	EEG KComplex	High -Middle -Low	2
	EEG Delta	High -Low	1
	EEG Theta	High -Low	1
	EEG Stability	Stable -Not Confident -Unstable	2
	EOG EyeMovement	High -Middle -Low -Lowest	3
	EOG Correlation	Conjugate -Disconjugate	1
	EMG Activity	High -Middle -Low	2

  Step 7. Repeat: If the terminate condition is not satisfied, repeat from Step 2 to Step 6 until one of the terminate conditions is satisfied.

	• Step 4. Selection: Rank the values of fitness function
	and select elite samples in term of F-Measure. In this
	step, a number of ρN samples with higher F-Measure are
	selected as elite samples.

• Step 5. Updating: Update the PDFs parameters µ and σ. According to the elite samples, new µ and σ are calculated. • Step 6. Check terminate conditions: If one of terminate conditions like, F-Measure reaches to a desired value or iteration reaches to the pre-defined value or the standard deviation σ is close to zero, is satisfied, then the procedure stops. •

  . It shows the increased F-Meausre value of each patient by applying CEATS in comparison to MTS method.

	Database	W	N1	N2	N3	R
	Patient 1	+0.540	+0.147	+0.193	+0.782	+0.271
	Patient 2	+0.507	+0.087	+0.152	+0.796	+0.832
	Patient 3	+0.118	+0.105	+0.308	+0.716	+0.111
	Patient 4	+0.235	+0.055	+0.495	+0.530	+0.218
	Patient 5	+0.144	+0.504	+0.430	+0.742	+0.517
	Patient 6	+0.412	+0.185	+0.339	+0.620	+0.321
	Patient 7	+0.481	+0.057	+0.283	+0.724	+0.687
	Patient 8	+0.509	+0.067	+0.418	+0.840	+0.759
	Patient 9	+0.235	+0.155	+0.181	+0.834	+0.714
	Patient 10	+0.552	+0.129	+0.464	+0.824	+0.777
	Patient 11	+0.580	+0.073	+0.225	+0.757	+0.728
	Patient 12	+0.477	+0.080	+0.451	+0.808	+0.327
	Patient 13	+0.436	+0.169	+0.359	+0.638	+0.249
	Patient 14	+0.464	+0.026	+0.223	+0.770	+0.667
	Table II: F-Measure Improvements by applying CEATS
	In Table


Anne Benoit

Aurélien Cavelan

Valentin Le Fèvre

Yves Robert

Hongyang Sun

Towards Optimal Multi-Level Checkpointing

Keywords: resilience, fail-stop errors, multi-level checkpointing, optimal pattern résilience, erreurs fatales, checkpoint multi-niveaux, schéma optimal

We provide a framework to analyze multi-level checkpointing protocols, by formally defining a k-level checkpointing pattern. We provide a first-order approximation to the optimal checkpointing period, and show that the corresponding overhead is in the order of k =1 √ 2λ C , where λ is the error rate at level , and C the checkpointing cost at level . This nicely extends the classical Young/Daly formula on single-level checkpointing. Furthermore, we are able to fully characterize the shape of the optimal pattern (number and positions of checkpoints), and we provide a dynamic programming algorithm to determine the optimal subset of levels to be used. Finally, we perform simulations to check the accuracy of the theoretical study and to confirm the optimality of the subset of levels returned by the dynamic programming algorithm. The results nicely corroborate the theoretical study, and demonstrate the usefulness of multi-level checkpointing with the optimal subset of levels.

Introduction

Checkpointing is the de-facto standard resilience method for HPC platforms at extreme-scale. However, the traditional single-level checkpointing method suffers from significant overhead, and multi-level checkpointing protocols now represent the state-of-the-art technique. These protocols allow different levels of checkpoints to be set, each with a different checkpointing overhead and recovery ability. Typically, each level corresponds to a specific fault1 type, and is associated to a storage device that is resilient to that type. For instance, a two-level system would deal with (i) transient memory errors (level 1) by storing key data in main memory; and (ii) node failures (level 2) by storing key data in stable storage (remote redundant disks).

We consider a very general scenario, where the platform is subject to k levels of faults, numbered from 1 to k. Level is associated with an error rate λ , a checkpointing cost C , and a recovery cost R . A fault at level destroys all the checkpoints of lower levels (from 1 to -1 included) and implies a roll-back to a checkpoint of level or higher. Similarly, a recovery of level will restore data from all lower levels. Typically, fault rates are decreasing and checkpoint/recovery costs are increasing when we go to higher levels:

λ 1 ≥ λ 2 ≥ • • • ≥ λ k , C 1 ≤ C 2 ≤ • • • ≤ C k , and R 1 ≤ R 2 ≤ • • • ≤ R k .
The idea of multi-level checkpointing is that checkpoints are taken for each level of faults, but at different periods. Intuitively, the less frequent the faults, the longer the checkpointing period: this is because the risk of a failure striking is lower when going to higher levels; hence the expected re-execution time is lower too; one can safely checkpoint less frequently, thereby reducing failure-free overhead (checkpointing is useless in the absence of fault). There are several natural approaches to implement multi-level checkpointing. The first option is to use independent checkpointing periods for each level. This option raises several difficulties, the most prominent one being overlapping checkpoints. Typically, we need to checkpoint different levels in sequence (e.g., writing into memory before writing onto disk), so we would need to delay some checkpoints, which might not be possible in some environments, and which would introduce irregular periods. The second option is to synchronize all checkpoint levels by nesting them inside a periodic pattern that repeats over time, as illustrated in Figure 1(a). In this figure, the pattern has five computational segments, each followed by a level-1 checkpoint. A segment is a chunk of work between two checkpoints, and a pattern consists in segments and checkpoints. The second and fifth level-1 checkpoints are followed by a level-2 checkpoint. Finally, the pattern ends with a level-3 checkpoint. When using patterns, a checkpoint at level is always preceded by checkpoints at all lower levels 1 to -1, which makes good sense in practice (e.g., with two levels, main memory and disk, one writes the data into memory before transferring it to disk).

Using periodic patterns simplifies the orchestration of checkpoints at all levels. In addition, repeatedly applying the same pattern is optimal for on-line scheduling problems, or for jobs running a very long (even infinite) time on the platform. Indeed, in this scenario, we seek the best pattern, i.e., the one whose overhead is minimal. The overhead of a pattern is the price per work unit to pay for resilience in the pattern; hence minimizing overhead is equivalent to optimizing platform throughput. For a pattern P(W ) with W units of work (the cumulated length of all its segments), the overhead H(P(W)) is defined as the ratio of the pattern's expected execution time E(P(W)) over its total work W minus 1:

Time C 1 C 2 C 3 C 1 C 1 C 2 C 1 C 1 C 1 C 2 C 3 (a) Time C 1 C 1 (b) Time C 1 C 2 C 1 C 1 C 2 C 1 C 1 C 2 (c)
H(P(W )) = E(P(W )) W -1. (1)
If there were neither checkpoint nor fault, the overhead would be zero. Determining the optimal pattern (with minimal overhead), and then repeatedly using it until job completion, is the optimal approach with Exponential failure distributions and long-lasting jobs. Indeed, once a pattern is successfully executed, the optimal strategy is to re-execute the same pattern. This is because of the memoryless property of exponential distributions: the history of failures has no impact on the solution, so if a pattern is optimal at some point in time, it stays optimal later in the execution, because we have no further information about the amount of work still to be executed. The difficulty of characterizing the optimal pattern dramatically increases with the number of levels. How many checkpoints of each level should be used, and at which locations inside the pattern? What is the optimal length of each segment? With one single level (see Figure 1(b)), there is a single segment of length W , and the Young/Daly formula [START_REF] Young | A first order approximation to the optimum checkpoint interval[END_REF][START_REF] Daly | A higher order estimate of the optimum checkpoint interval for restart dumps[END_REF] gives W opt = 2C 1 λ 1 . The minimal overhead is then H opt = √ 2λ 1 C 1 + O(λ 1 ) [START_REF] Benoit | Optimal resilience patterns to cope with fail-stop and silent errors[END_REF]. With two levels, the pattern still has a simple shape, with N segments followed by a level-1 checkpoints, and ended by a level-2 checkpoint (see Figure 1(c)). Recent work [START_REF] Di | Toward an optimal online checkpoint solution under a two-level HPC checkpoint model[END_REF] shows that all segments have same length in the optimal pattern, and provides mathematical equations that can be solved numerically to compute both the optimal length W opt of the pattern and its optimal number of segments. However, no closed-form expression is available, neither for W opt , nor for the minimal overhead H opt .

With three levels, no optimal solution is known. The pattern shape becomes quite complicated. Coming back to Figure 1(a), we identify two sub-patterns ending with a level-2 checkpoint. The first sub-pattern has 2 segments while the second one has 3. The memoryless property does not imply that all sub-patterns are identical, because the state after completing the first sub-pattern is not the same as the initial state when beginning the execution of the pattern. In the general case with k levels, the shape of the pattern will be even more complicated, with different-shaped sub-patterns (each ended by a level k-1 checkpoint). In turn, each sub-pattern may have different-shaped sub-sub-patterns (each ended by a level k -2 checkpoint), and so on. The major contribution of this work is to provide an analytical characterization of the optimal pattern with an arbitrary number k of checkpointing levels, with closed-form formulas for the pattern length W opt , the number of checkpoints at each level, and the optimal overhead H opt . In particular, we obtain the following beautiful result:

H opt = k =1 2λ C + O(Λ), (2)
where Λ = k =1 λ . However, we point out that this analytical characterization relies on a first-order approximation, so it is valid only when resilience parameters C and R are small in front of the platform Mean Time Between Failures (MTBF) µ = 1/Λ. Also, the optimal pattern has rational number of segments, and we use rounding to derive a practical solution. Still, Equation ( 2) provides a lower bound on the optimal overhead, and this bound is met very closely in all our experimental scenarios.

Finally, in many practical cases, there is no obligation to use all available checkpointing levels. For instance, with k = 3 levels, one may choose among four possibilities: level 3 only, levels 1 and 3, levels 2 and 3, and all levels 1, 2 and 3. Of course, we still have to account for all failure types, which translates into the following:

• level 3: use

λ 3 ← λ 1 + λ 2 + λ 3 ;
• levels 1 and 3: use λ 1 and λ 3 ← λ 2 + λ 3 ;

• levels 2 and 3: use λ 2 ← λ 1 + λ 2 and λ 3 ;

• all levels: use λ 1 , λ 2 and λ 3 . Our analytical characterization of the optimal pattern leads to a simple dynamic programming algorithm for selecting the optimal subset of levels.

The rest of this paper is organized as follows. Section 2 surveys the related work. Section 3 is the heart of the paper and shows how to compute the optimal pattern as well as the optimal subset of levels. Section 4 is devoted to simulations assessing the accuracy of the first-order approximation. Finally, Section 5 provides concluding remarks and hints for future work.

Related work

Given the checkpointing cost and platform MTBF, classical formulas due to Young [START_REF] Young | A first order approximation to the optimum checkpoint interval[END_REF] and Daly [START_REF] Daly | A higher order estimate of the optimum checkpoint interval for restart dumps[END_REF] are well known to determine the optimal checkpointing period in the single-level checkpointing scheme. However, this method suffers from the intrinsic limitation that the cost of checkpointing/recovery grows with failure probability, and becomes unsustainable at large scale [START_REF] Ferreira | Evaluating the Viability of Process Replication Reliability for Exascale Systems[END_REF][START_REF] Bosilca | Unified model for assessing checkpointing protocols at extreme-scale[END_REF] (even with diskless or incremental checkpointing [START_REF] Plank | Diskless checkpointing[END_REF]).

To reduce the I/O overhead, various two-level checkpointing protocols have been studied. Vaidya [START_REF] Vaidya | A case for two-level distributed recovery schemes[END_REF] proposed a two-level recovery scheme that tolerates a single node failure using a local checkpoint stored on a parter node. If more than one failure occurs during any local checkpointing interval, the scheme resorts to the global checkpoint. Silva and Silva [START_REF] Silva | Using two-level stable storage for efficient checkpointing[END_REF] advocated a similar scheme by using memory to store local checkpoints, which is protected by XOR encoding. Di et al. [START_REF] Di | Toward an optimal online checkpoint solution under a two-level HPC checkpoint model[END_REF] analyzed a two-level checkpointing pattern, and proved equal-length segments in the optimal solution. They also provided mathematical equations that can be solved numerically to compute the optimal pattern length and number of segments. Benoit et al. [START_REF] Benoit | Optimal resilience patterns to cope with fail-stop and silent errors[END_REF] relied on disk checkpoints to cope with fail-stop failures and memory checkpoints coupled with error detectors to handle silent data corruptions. They derived first-order approximation formulas for the optimal pattern length and the number of memory checkpoints between two disk checkpoints. Some authors have also generalized two-level checkpointing to account for an arbitrary number of levels. Moody et al. [START_REF] Moody | Design, Modeling, and Evaluation of a Scalable Multi-level Checkpointing System[END_REF] implemented this approach in a three-level Scalable Checkpoint/Restart (SCR) library. They relied on a rather complex Markov model to recursively compute the efficiency of the scheme. Bautista-Gomez et al. [START_REF] Bautista-Gomez | FTI: High performance fault tolerance interface for hybrid systems[END_REF] designed a four-level checkpointing library, called Fault Tolerance Interface (FTI), in which partner-copy and Reed-Solomon coding are employed as two intermediate levels between local and global disks. Based on FTI, Di et al. [START_REF] Di | Optimization of multi-level checkpoint model for large scale HPC applications[END_REF] proposed an iterative method to compute the optimal checkpointing interval for each level with prior knowledge of the application's total execution time. Hakkarinen and Chen [START_REF] Hakkarinen | Multilevel diskless checkpointing[END_REF] considered multi-level diskless checkpointing for tolerating simultaneous failures of multiple processors. Balaprakash et al. [START_REF] Balaprakash | Analysis of the tradeoffs between energy and run time for multilevel checkpointing[END_REF] studied the trade-off between performance and energy for general multi-level checkpointing schemes.

While all of these works relied on numerical methods to compute the checkpointing intervals at different levels, this paper is the first one to provide explicit formulas on the optimal parameters in a multi-level checkpointing protocol (up to first-order approximation as in Young/Daly's classical result).

Computing the optimal pattern

This section computes the optimal multi-level checkpointing pattern. We first state our assumptions in Section 3.1, and then analyze the simple case with k = 2 levels in Section 3.2, before proceeding to the general case in Section 3.3. Finally, the algorithm to compute the optimal subset of levels is described in Section 3.4.

Assumptions

In this paper, we assume that failures from different levels are independent2 . For each level , the arrival of failures follows Poisson process with error rate λ . In order to deal with the interplay of failures from different levels, we make use of the following well-known properties of independent Poisson processes [START_REF] Gallager | Stochastic Processes: Theory for Applications[END_REF]Chapter 2.3].

Property 1. During the execution of a segment with length w, let X denote the time when the first level-error strikes. Thus, X is a random variable following an Exponential distribution with parameter λ , for all = 1, 2, . . . , k.

(1). Let X denote the time when the first error (of any level) strikes. We have X = min{X 1 , X 2 , . . . , X k }, which follows an Exponential distribution with parameter Λ = k =1 λ . The probability of having an error (from any level) in the segment is therefore

P (X ≤ w) = 1 -e -Λw .
(2). Given that an error (from any level) strikes during the execution of the segment, the probability that the error belongs to a particular level is proportional to the error rate of that level, i.e., P (X = X |X ≤ w) = λ Λ , for all = 1, 2, . . . , k. Moreover, we assume that error rates of different levels are of the same order, i.e., λ = Θ(Λ) for all = 1, 2, . . . , k, and that errors only strike during the computations, while checkpointing and recovery are error-free. Indeed, the durations of checkpoints and recoveries are generally small compared to the pattern length, so the probability of a failure striking during these operations is low. It has been shown in [START_REF] Benoit | Optimal resilience patterns to cope with fail-stop and silent errors[END_REF] that removing this assumption does not impact the first-order approximation of the pattern overhead.

Optimal two-level pattern

We start by analyzing the two-level pattern shown in Figure 1(b). The goal is to determine a first-order approximation to the optimal pattern length W , the number n of level-1 checkpoints in the pattern, as well as the length w i = α i W of the i-th segment, for all 1 ≤ i ≤ n, where n i=1 α i = 1.

With a single segment

We first consider a special case of the two-level pattern, in which only a single segment is present, i.e., n = 1. The result establishes the order of the optimal pattern length W opt , which will be used later for analyzing the general case. Recall that Λ = λ 1 + λ 2 and, for convenience, let us also define

C = C 1 + C 2 .
The following proposition shows the expected time of such a pattern with fixed length W .

Proposition 1. The expected execution time of a two-level pattern with a single segment and fixed length W is

E = W + C + 1 2 ΛW 2 + O(max{Λ 2 W 3 , ΛW }).
Proof. We can express the expected execution time of the pattern recursively as follows:

E = P E lost (W, Λ) + λ 1 Λ (R 1 + E) + λ 2 Λ R 2 + R 1 + E + (1 -P ) (W + C) , (3)
where P = 1 -e -ΛW denotes the probability of having a failure (either level-1 or level-2) during the execution of the pattern based on Property 1.1, and E lost (w i , Λ) denotes the expected time lost when such a failure occurs. In this case, and based on Property 1.2, if the failure belongs to level 1, which happens with probability λ 1 Λ , we can recover from the latest level-1 checkpoint (R 1 ). Otherwise, the failure belongs to level 2 with probability λ 2 Λ , and we need to first recover from the latest level-2 checkpoint (R 2 ) before restoring the level-1 checkpoint (R 1 ). In both cases, the entire pattern needs to be re-executed again. Finally, if no error (of any level) strikes, which happens with probability 1 -P , the pattern is completed after W time of execution followed by the time C to perform the two checkpoints, which are assumed to be error-free.

From [13, Equation (1.13)], the expected time lost when executing a segment of length W with error rate Λ is

E lost (W, Λ) = 1 Λ - W e ΛW -1 . (4)
Substituting Equation (4) into Equation (3) and solving for E, we get:

E = e ΛW -1 1 Λ + R 1 + λ 2 Λ R 2 + C 1 + C 2 , (5)
which is an exact formula on the expected execution time of the pattern. Now, using Taylor series to expand e ΛW = 1 + ΛW + Λ 2 W 2

2

+ O(Λ 3 W 3 ) while assuming W = Θ(Λ -x ), where RR n°8930 0 < x < 1, we can re-write Equation ( 5) as

E = W + 1 2 ΛW 2 + C 1 + C 2 + O(Λ 2 W 3 ) + ΛW + Λ 2 W 2 2 + O(Λ 3 W 3 ) R 1 + λ 2 Λ R 2 .
Since recovery costs (R 1 , R 2 ) are assumed to be constants, and error rates (λ 1 , λ 2 , Λ) are in the same order, the expected execution time can be expressed as follows:

E = W + C 1 + C 2 + 1 2 ΛW 2 + O(Λ 2 W 3 ) + O(ΛW ),
which completes the proof of the proposition.

From Proposition 1, the expected execution overhead of the pattern can be derived as

H = C W + 1 2 ΛW + O(max{Λ 2 W 2 , Λ}).
Assume that the platform MTBF µ = 1/Λ is large in front of the resilience parameters, and consider the first two terms of H: the overhead is minimized when the pattern has length W = Θ(Λ -1/2 ), and in that case both terms are in the order of Θ(Λ 1/2 ), so we have

H = Θ(Λ 1/2 ) + O(Λ). Indeed, the last term O(Λ 2 W 2 ) = O(Λ) becomes negligible compared to Θ(Λ 1/2
). Hence, the optimal pattern length W opt can be obtained by balancing the first two terms in H, which gives

W opt = 2C Λ = Θ(Λ -1/2 ), (6)
and the optimal execution overhead becomes

H opt = √ 2ΛC + O(Λ). (7)
Remarks. Unlike in single-level checkpointing, the checkpoint to roll back to in a twolevel pattern depends on which type of error strikes first. Under first-order approximation and assuming that the resilience parameters are small compared to the platform MTBF and pattern length, the formulas shown in Equations ( 6) and ( 7) reduce exactly to Young/Daly's classical result by aggregating the error rates and checkpointing costs of both levels.

With multiple segments

We now consider the general two-level pattern with multiple segments, and derive the optimal pattern parameters. As in the single-segment case, we start with a proposition showing the expected time to execute a two-level pattern with fixed parameters.

Proposition 2. The expected execution time of a given two-level pattern is

E = W +nC 1 +C 2 + 1 2 λ 1 n i=1 α 2 i + λ 2 W 2 + O(Λ 1/2 ).
Proof. We first prove the following result (by induction) on the expected time E i to execute the i-th segment of the pattern (up to the level-1 checkpoint at the end of the segment):

E i = w i + C 1 + λ 1 2 w 2 i + λ 2   w 2 i 2 + i-1 j=1 w j w i   + O(Λ 1/2 ). ( 8
)
According to the result with a single segment, we know that the optimal pattern length and hence the segment length are in the order of O(Λ -1/2 ), which implies that

E i = w i + O(1).
For the ease of analysis, we assume that there is a hypothetical segment at the beginning of the pattern with length w 0 = 0 (hence no need to checkpoint). For this segment, we have E 0 = w 0 = 0, satisfying Equation [START_REF] Di | Toward an optimal online checkpoint solution under a two-level HPC checkpoint model[END_REF]. Suppose the claim holds up to E i-1 . Then, E i can be recursively expressed as follows:

E i = P i E lost (w i , Λ) + λ 1 Λ (R 1 + E i ) + λ 2 Λ R 2 + R 1 + i-1 j=1 E j + E i + (1 -P i )(w i + C 1 ), (9)
where P i = 1 -e -Λw i denotes the probability of having a failure (either level-1 or level-2) during the execution of the segment, and E lost (w i , Λ) denotes the expected time lost when such a failure occurs. Equation ( 9) is very similar to Equation (3), except when a level-2 failure occurs we need to re-execute all the segments (up to segment i) that have been executed so far. Following the derivation of Proposition 1 and applying E j = w j + O(1) for j = 1, 2, . . . , i -1, we can derive the first-order approximation of E i as follows:

E i = w i +C 1 + 1 2 λ 1 w 2 i +λ 2 w 2 i +2λ 2 w i i-1 j=1 E j +O(Λ 1/2 ) = w i +C 1 + 1 2 λ 1 w 2 i +λ 2 w 2 i +2λ 2 w i i-1 j=1 w j +O(1) +O(Λ 1/2 ) = w i +C 1 + 1 2 λ 1 w 2 i +λ 2 w 2 i +2 i-1 j=1 w j w i +O(Λ 1/2 ). ( 10
)
Since the level-2 checkpoint at the end of the pattern is also assumed to be error-free, we can compute the expected execution time of the pattern as

E = n i=1 E i + C 2 = W + nC 1 + C 2 + 1 2 λ 1 n i=1 α 2 i + λ 2 W 2 + O(Λ 1/2 ), since n i=1 w 2 i + 2 n i=1 i-1 j=1 w j w i = ( n i=1 w i ) 2 = W 2 . RR n°8930
Theorem 1. A first-order approximation to the optimal two-level pattern is characterized by

n opt = λ 1 λ 2 • C 2 C 1 , (11)
α opt i = 1 n opt ∀i = 1, 2, . . . , n opt , (12)
W opt = n opt C 1 + C 2 1 2 λ 1 n opt + λ 2 , ( 13
)
where n opt is the number of segments, α opt i W opt is the length of the i-th segment, and W opt is the pattern length. The optimal pattern overhead is

H opt = 2λ 1 C 1 + 2λ 2 C 2 + O(Λ). (14)
Proof. For a given pattern with a fixed number n of segments, n i=1 α 2 i is minimized subject to

n i=1 α i = 1 when α i = 1
n for all i = 1, 2, . . . , n. Hence, we can derive the expected execution overhead from Proposition 2 as follows:

H = nC 1 + C 2 W + 1 2 λ 1 n + λ 2 W + O(Λ). (15)
For a given n, the optimal work length can then be computed from Equation [START_REF] Plank | Diskless checkpointing[END_REF], and it is given by

W opt = nC 1 +C 2 1 2 λ 1 n +λ 2
. In that case, the execution overhead becomes

H = 2 λ 1 n + λ 2 (nC 1 + C 2 ) + O(Λ), (16)
which is minimized as shown in Equation ( 14) when n satisfies Equation [START_REF] Hakkarinen | Multilevel diskless checkpointing[END_REF]. Indeed,

2 λ 1 n opt + λ 2 n opt C 1 + C 2 = 2λ 1 C 1 + 2λ 2 C 2 + 4 √ λ 1 λ 2 C 1 C 2 = ( √ 2λ 1 C 1 + √ 2λ 1 C 1 ) 2 .
In practice, since the number of segments can only be a positive integer, the optimal solution is either max(1, n opt ) or n opt , whichever leads to a smaller value of the convex function H as shown in Equation ( 16). [START_REF] Di | Toward an optimal online checkpoint solution under a two-level HPC checkpoint model[END_REF] with

Remarks. Consider the example given in

C 1 = R 1 = 20, C 2 = R 2 = 50, λ 1 = 2.
78 × 10 -4 and λ 2 = 4.63 × 10 -5 . The optimal solution3 provided by [START_REF] Di | Toward an optimal online checkpoint solution under a two-level HPC checkpoint model[END_REF] gives n opt = 3.83, W opt = 1362.49 and H opt = 0.1879, while Theorem 1 suggests n opt = 3.87, W opt = 1378.27 and H opt = 0.1735, which is quite close to the exact optimum. The difference in overhead is due to the negligence of lower-order terms in the first-order approximation. We point out that the solution provided by [START_REF] Di | Toward an optimal online checkpoint solution under a two-level HPC checkpoint model[END_REF] relies on numerical methods to solve rather complex mathematical equations, whose convergence is not always guaranteed, and it is only applicable to two levels. Our result, on the other hand, is able to provide fast and good approximation to the optimal solution when the error rates are sufficiently small, and it can be readily extended to an arbitrary number of levels, as shown in the next section.

Optimal k-level pattern

In this section, we derive the first-order approximation to the optimal k-level pattern by determining its length W , the number N of level-checkpoints for all 1 ≤ ≤ k, as well as the positions of all checkpoints in the pattern.

Observations

Before analyzing the optimal pattern, we make several observations. First, we can obtain the orders of the optimal length and pattern overhead as shown below (recall that Λ = k =1 λ ). Observation 1. Consider the simplest k-level pattern with a single segment of length W . We can conduct the same analysis as in Section 3.2.1 to show that the optimal pattern length satisfies W opt = Θ(Λ -1/2 ), and the corresponding overhead satisfies

H opt = Θ(Λ 1/2 ).
From the analysis of the two-level pattern, we can also observe that the overall execution overhead of any pattern comes from two distinct sources defined below.

Observation 2. There are two types of execution overheads for a pattern:

(1). Error-free overhead, denoted as o ef , is the total cost of all the checkpoints placed in the pattern. For a given set of checkpoints, the error-free overhead is completely determined regardless of their positions in the pattern.

(2). Re-executed fraction overhead, denoted as o re , is the expected fraction of work that needs to be re-executed due to errors. The re-executed fraction overhead depends on both the set of checkpoints and their positions.

For example, in the two-level pattern with n level-1 checkpoints and given values of α i for all i = 1, 2, . . . , n, the two types of overheads are given by o ef = nC 1 + C 2 and o re =

1 2 f 1 n i=1 α 2 i + f 2 , where f = λ Λ for = 1, 2.
Assuming that checkpoints at all levels have constant costs and that the error rates at all levels are in the same order, then both o ef and o re can be considered as constants, i.e., o ef = O(1) and o re = O(1).

A trade-off exists between these two types of execution overheads, since placing more checkpoints generally reduces the re-executed work fraction when an error strikes, but it can adversely increase the overhead when the execution is error-free. Therefore, in order to achieve the best overall overhead, a resilience algorithm must seek an optimal balance between o ef and o re .

For a given pattern with fixed overheads o ef and o re , we can make the following observation based on Propositions 1 and 2, which partially characterizes the optimal pattern. Observation 3. For a given pattern (with fixed o ef and o re ), the expected execution time is given by

E = W + o ef error-free execution time + ΛW expected # errors • o re W re-executed work in case of error + O(Λ 1/2 ), ( 17
)
and the optimal pattern length and the resulting expected execution overhead of the pattern are

W opt = o ef Λ • o re , (18)
H opt = 2 Λ • o ef • o re + O(Λ). (19)
s (4) s Equation ( 19) shows that the trade-off between o ef and o re is manifested as the product of the two terms. Hence, in order to determine the optimal pattern, it suffices to find the pattern parameters (e.g., n and α i ) that minimize o ef • o re .

(3) 1 s (3) 2 c 4 c 3 c 3 s (2) 1,1 s (2) 1,2 s (2) 1,3 s (2) 1,4 c 2 s (2) 2,1 s (2) 2,2 c 2 c 2 c 2 c 2 c 1 c 1 c 1 s (1) 1,3,1 s (1) 1,3,2 s (1) 1,3,3 c 4

Analysis

We now extend the analysis to derive the optimal multi-level checkpointing patterns. Generally, for a k-level pattern, each computational segment s ( ) i k-1 ,...,i can be uniquely identified by its level as well as its position i k-1 , . . . , i within the multi-level hierarchy. For instance, in a four-level pattern, the segment s

(2)
1,3 denotes the third level-2 segment inside the first level-3 segment of the pattern (see Figure 2). Note that a segment can contain multiple sub-segments at the lower levels (except for bottom-level segments) and is a sub-segment of a larger segment at a higher level (except for top-level segments). The entire pattern can be denoted as s (k) , which is the only segment at level k.

For any segment s

( ) i k-1 ,...,i at level , where 1 ≤ ≤ k, let w ( ) i k-1 ,.
..,i denote its length. Hence, we have w

( +1) i k-1 ,...,i +1 = i w ( ) i k-1 ,...,i and w (k) = W . Also, let n ( ) i k-1 ,.
..,i denote the number of sub-segments contained by s

( ) i k-1 ,...,i at the lower level -1. We have n (1) i k-1 ,...,i 1 = 1 for all i k-1 , . . . , i 1 . For convenience, we further define α ( ) i k-1 ,...,i = w ( ) i k-1 ,...,i

W

as the fraction of the length of segment s ( ) i k-1 ,...,i inside the pattern, and define N to be the total number of level-segments in the entire pattern. Therefore, we have

N k = 1, N k-1 = n (k) , and in general N = i k-1 ,...,i +1 n ( +1) i k-1 ,...,i +1 .
The following proposition shows the expected time to execute a given k-level pattern.

Proposition 3. The expected execution time of a given k-level pattern is

E = W + k-1 =1 N C + C k + W 2 2   k =1 λ i k-1 ,...,i α ( ) i k-1 ,...,i 2   + O(Λ 1/2 ).
Proof. We show that the expected time to execute any segment s

(h) i k-1 ,.
..,i h at level h, where 1 ≤ h ≤ k, satisfies the following (without counting the time to execute all the checkpoints inside the segment):

E (h) i k-1 ,...,i h = w (h) i k-1 ,...,i h + W 2 2   h =1 λ i h-1 ,...,i α ( ) i k-1 ,...,i 2   + Λ [h+1,k]    w (h) i k-1 ,...,i h 2 2 + w (h) i k-1 ,...,i h i h -1 j h =1 E (h) i k-1 ,...,j h    + w (h) i k-1 ,...,i h k =h+2 Λ [ ,k] i -1 -1 j -1 =1 E ( -1) i k-1 ,...,j -1 +O(Λ 1/2 ), (20)
where Λ [x,y] = y =x λ and, if x > y, we define Λ [x,y] = 0. The proposition can then be proven by setting E = E (k) + k-1 =1 N C + C k , since checkpoints are assumed to be error-free. We now prove Equation (20) by induction on the level h. For the base case, i.e., when h = 1, consider a segment s (1) i k-1 ,...,i 1 at the first level. Following the proof of Proposition 2 (in particular, Equation ( 9)), we can express its expected execution time E

(1)

i k-1 ,...,i 1 , as E (1) i k-1 ,...,i 1 = P (1) i k-1 ,...,i 1 E lost w (1) i k-1 ,...,i 1 , Λ + λ 1 Λ R 1 + E (1) i k-1 ,...,i 1 + λ 2 Λ 2 j=1 R j + i 1 j 1 =1 E (1) i k-1 ,...,j 1 + λ 3 Λ 3 j=1 R j + i 2 -1 j 2 =1 E (2) i k-1 ,...,j 2 + i 1 j 1 =1 E (1) i k-1 ,...,j 1 . . . + λ k Λ k j=1 R j + i k-1 -1 j k-1 =1 E (k-1) j k-1 + i k-2 -1 j k-2 =1 E (k-2) i k-1 ,j k-2 + • • • + i 1 j 1 =1 E (1) i k-1 ,...,j 1 + 1 -P (1)
i k-1 ,...,i 1 w

(1)

i k-1 ,...,i 1 , (21)
where Λ = k =1 λ is the total rate of all error sources, and P

(1)

i k-1 ,...,i 1 = 1-e Λ•w (1)
i k-1 ,...,i 1 denotes the probability of having an error (from any level) during the execution of the segment. Simplifying Equation (21) and solving for E (1) i k-1 ,...,i 1 we get:

E (1) i k-1 ,...,i 1 = w (1) i k-1 ,...,i 1 + W 2 2 Λ [1,k] α (1) i k-1 ,...,i 1 2 + w (1) i k-1 ,...,i 1 k =2 Λ [ ,k] i -1 -1 j -1 =1 E ( -1) i k-1 ,...,j -1 + O(Λ 1/2 ) = w (1)
i k-1 ,...,i 1 + W 2 2 λ 1 α (1) i k-1 ,...,i 1 2 + Λ [2,k]    w (1) i k-1 ,...,i 1 2 2 + w (1) i k-1 ,...,i 1 i 1 -1 j 1 =1 E (1) i k-1 ,...,j 1    + w (1) i k-1 ,...,i 1 k =3 Λ [ ,k] i -1 -1 j -1 =1 E ( -1) i k-1 ,...,j -1 +O(Λ 1/2 ),
which satisfies Equation (20). Suppose Equation (20) holds up to any segment s

i k-1 ,...,i h at level h. Following the proof of Proposition 2 (in particular, the derivation of Equation ( 10)), we can show by induction that E (h)

i k-1 ,...,i h = w (h) i k-1 ,...,i h + O(1)
. Hence, for segment s (h+1) i k-1 ,...,i h+1 at level h + 1, we have:

E (h+1) i k-1 ,...,i h+1 = i h E (h) i k-1 ,...,i h = i h w (h) i k-1 ,...,i h + W 2 2   h =1 λ i h ,...,i α
( ) i k-1 ,...,i 2   + Λ [h+1,k] i h    w (h) i k-1 ,...,i h 2 2 +w (h) i k-1 ,...,i h i h -1 j h =1 w (h) i k-1 ,...,j h    + i h w (h) i k-1 ,...,i h k =h+2 Λ [ ,k] i -1 -1 j -1 =1 E ( -1) i k-1 ,...,j -1 + O(Λ 1/2 ) = w (h+1) i k-1 ,...,i h+1 + W 2 2   h =1 λ i h ,...,i α
( ) i k-1 ,...,i 2   + Λ [h+1,k] w (h+1) i k-1 ,...,i h+1 2 2 + w (h+1) i k-1 ,...,i h+1 k =h+2 Λ [ ,k] i -1 -1 j -1 =1 E ( -1) i k-1 ,...,j -1 + O(Λ 1/2 ) RR n°8930 = w (h+1) i k-1 ,...,i h+1 + W 2 2   h+1 =1 λ i h ,...,i α ( ) i k-1 ,...,i 2   + Λ [h+2,k]    w (h+1) i k-1 ,...,i h+1 2 2 +w (h+1) i k-1 ,...,i h+1 i h+1 -1 j h+1 =1 E (h+1) i k-1 ,...,j h+1    + w (h+1) i k-1 ,...,i h+1 k =h+3 Λ [ ,k] i -1 -1 j -1 =1 E ( -1) i k-1 ,...,j -1 + O(Λ 1/2 ).
Hence, Equation (20) also holds for any segment at level h + 1. This completes the proof of the proposition.

Proposition 3 shows that, for a given k-level checkpointing pattern, the error-free overhead o ef and the re-executed fraction overhead o re are given as follows:

o ef = k-1 =1 N C + C k , (22)
o re = 1 2 k =1 f i k-1 ,...,i α
( ) i k-1 ,...,i 2 , (23)
where f = λ Λ . According to Observation 3, it remains to find parameters of the pattern such that o ef • o re is minimized.

To derive the optimal pattern, we first consider the case where o ef is fixed, i.e., the set of checkpoints is given. The following proposition shows the optimal value of o re . Proposition 4. For a k-level checkpointing pattern, suppose the number N of checkpoints at each level is given, i.e., the error-free overhead o ef is fixed (as in Equation (22)). Then, the optimal value of the re-executed work overhead is given by

o opt re = 1 2 k-1 =1 f N + f k , ( 24
)
and it is obtained when all the checkpoints of each level are equally spaced in the pattern.

Proof. According to Equation (23), which shows the value of o re for the entire pattern, we can define the corresponding overhead for each level-h segment s (h) i k-1 ,...,i h recursively as follows:

o re s (h) i k-1 ,...,i h = f h 2 • α (h) i k-1 ,...,i h 2 + i h-1 o re s (h-1) i k-1 ,...,i h-1 , with o re s (0)
i k-1 ,...,i 0 = 0 by definition. For each segment s (h) i k-1 ,...,i h , we also define N s (h) i k-1 ,...,i h to be the total number of levelsegments it contains, with ≤ h. We will show that the optimal value o opt re s (h) i k-1 ,...,i h for the segment satisfies:

o opt re s (h) i k-1 ,...,i h = 1 2   h =1 f N s (h) i k-1 ,...,i h   α (h) i k-1 ,...,i h 2 , ( 25
)
and it is achieved when its level-checkpoints are equally spaced, for all ≤ h -1. The proposition can then be proven by setting

o opt re = o opt re s (k) , since N s (k) = N , N k = 1, and
α (k) = 1.
Now, we prove Equation ( 25) by induction on the level h. For the base case, i.e., when h = 1, we have o re s

(1)

i k-1 ,...,i 1 = f 1 2 • α (1) i k-1 ,...,i 1 2
by definition, and it satisfies Equation (25), because N 1 s

(1)

i k-1 ,...,i 1 = 1. Suppose Equation (25) holds for any segment s (h) i k-1 ,...,i h at level h. Then, for segment s (h+1) i k-1 ,...,i h+1 at level h + 1, we have:

o re s (h+1) i k-1 ,...,i h+1 = f h+1 2 • α (h+1) i k-1 ,...,i h+1 2 + i h o opt re s (h) i k-1 ,...,i h = f h+1 2 • α (h+1) i k-1 ,...,i h+1 2 + 1 2 y, (26)
where

y = i h x (h) i k-1 ,...,i h • α (h) i k-1 ,...,i h 2
, and x

(h) i k-1 ,...,i h = h =1 f N s (h) i k-1 ,...,i h . To minimize o re s (h+1)
i k-1 ,...,i h+1 as shown in Equation ( 26), it suffices to solve the following minimization problem:

minimize y = i h x (h) i k-1 ,...,i h • α (h) i k-1 ,...,i h 2 , subject to i h α (h) i k-1 ,...,i h = α (h+1) i k-1 ,...,i h+1 .
Since y is clearly a convex function of α (h) i k-1 ,...,i h , we can readily get, using Lagrange multiplier [START_REF] Boyd | Convex Optimization[END_REF], the minimum value of y as follows:

y min = 1 i h 1/x (h) i k-1 ,...,i h • α (h+1) i k-1 ,...,i h+1 2 , (27)
which is obtained at

α (h) i k-1 ,...,i h = 1/x (h) i k-1 ,...,i h j h 1/x (h) i k-1 ,...,j h • α (h+1) i k-1 ,...,i h+1 . ( 28
) Let us define z = i h 1/x (h) i k-1 ,.
..,i h . We now need to solve the following maximization problem:

maximize z = i h 1 h =1 f N s (h) i k-1 ,...,i h , subject to i h N s (h) i k-1 ,...,i h = N s (h+1) i k-1 ,...,i h+1 , ∀ = 1, . . . , h.

Again, z is a convex function of N s (h)

i k-1 ,...,i h , and it can be shown to be maximized when

N s (h) i k-1 ,...,i h = N s (h+1) i k-1 ,...,i h+1 n (h+1) i k-1 ,...,i h+1 , ∀ = 1, . . . , h, which gives α (h) i k-1 ,...,i h = 1 n (h+1) i k-1 ,...,i h+1 • α (h+1)
i k-1 ,...,i h+1 according to Equation (28). This implies that all level-checkpoints are also equally spaced inside segment s (h+1) i k-1 ,...,i h+1 , for all ≤ h. The maximum value of z in this case is

z max = 1 h =1 f N s (h+1) i k-1 ,...,i h+1
, and the optimal value of y min according to Equation ( 27) is then given by

y opt min = 1 z max α (h+1) i k-1 ,...,i h+1 2 =   h =1 f N s (h+1) i k-1 ,...,i h+1   α (h+1) i k-1 ,...,i h+1 2 .
Substituting y opt min into Equation (26), we get the optimal value of o re s (h+1)

i k-1 ,...,i h+1 as follows:

o opt re s (h+1) i k-1 ,...,i h+1 = f h+1 2 • α (h+1) i k-1 ,...,i h+1 2 + 1 2 y opt min = 1 2   h+1 =1 f N s (h+1) i k-1 ,...,i h   α (h+1) i k-1 ,...,i h+1 2 .
This shows that Equation (25) also holds for segment s (h+1)

i k-1 ,...,i h+1 at level h + 1 and, hence, completes the proof of the proposition.

We are now ready to characterize the optimal k-level pattern. The result is stated in the following theorem. Theorem 2. A first-order approximation to the optimal k-level pattern and its overhead are characterized by

W opt = 2 k-1 =1 N opt C + C k k-1 =1 λ N opt + λ k , (29)
N opt = λ C • C k λ k , ∀ = 1, 2, . . . , k -1, (30)
H opt = k =1 2λ C + O(Λ). (31)
Proof. From Observation 3, Equation ( 22) and Proposition 4, we know that the optimal pattern can be obtained by minimizing the following function:

F = o ef • o opt re = 1 2 k-1 =1 N C + C k k-1 =1 f N + f k .
We first compute the optimal number of checkpoints at each level using a two-phase iterative method. Towards this end, let us define

o ef (h) = k-1 =h N C + C k , o opt re (h) = 1 2 k-1 =h f N + f k .
In the first phase, we set initially

F (1) = o ef (1) • o opt re (1).
The optimal value of N 1 that minimizes F (1) can then be obtained by setting

∂F (1) ∂N 1 = C 1 o opt re (1) -o ef (1)
f 1 2N 2 1 = C 1 f 1 2N 1 + o opt re (2) -(N 1 C 1 + o ef (2)) f 1 2N 2 1 = C 1 o opt re (2) -o ef (2) f 1 2N 2 1 = 0, which gives N opt 1 = f 1 C 1 • o ef (2)
2o opt re [START_REF] Bautista-Gomez | FTI: High performance fault tolerance interface for hybrid systems[END_REF] . Substituting it into F (1) and simplifying, we can get the value of F after the first iteration as

F (2) = 1 2 f 1 C 1 + o ef (2) • o opt re (2) 2 .
Repeating the above process, we can get the optimal value of F after k -1 iterations as

F opt = F (k) = 1 2 k =1 f C 2 , ( 32
)
and the optimal value of N as

N opt = f C • o ef ( + 1)
2o opt re ( + 1)

, ∀ = 1, 2, . . . , k -1. ( 33
)
In the second phase, we first compute from Equation ( 33)

N opt k-1 = f k-1 C k-1 • o ef (k) 2o opt re (k) = f k-1 C k-1 • C k f k = λ k-1 C k-1 • C k λ k .
Substituting it into N opt k-2 , we obtain:

N opt k-2 = f k-2 C k-2 • o ef (k -1) 2o opt re (k -1) = λ k-2 C k-2 • N opt k-1 C k-1 + C k λ k-1 N opt k-1 + λ k = λ k-2 C k-2 • λ k-1 λ k C k-1 C k + C k λ k-1 λ k C k-1 C k + λ k = λ k-2 C k-2 • C k λ k-1 λ k • C k-1 C k + 1 λ k λ k-1 λ k • C k-1 C k + 1 = λ k-2 C k-2 • C k λ k .
Repeating the above process iteratively, we can compute the optimal values of N opt for = k -3, . . . , 2, 1, as given in Equation (30) by using values of N opt k-1 , . . . , N opt +1 . The optimal pattern length, according to Equation [START_REF] Young | A first order approximation to the optimum checkpoint interval[END_REF], can be expressed as

W opt = o ef Λ•o opt re
, which turns out to be Equation (29) with the optimal values of N opt .

The optimal overhead, according to Equations ( 19) and (32), can be expressed as

H opt = 2 √ Λ • F opt +O(Λ)
, which gives rise to Equation (31). This completes the proof of the theorem.

Since Proposition 4 shows that all the checkpoints of each level are equally spaced in the pattern, we can readily obtain the following corollary.

Corollary 1. In an optimal k-level pattern, the number of level-checkpoints between any two consecutive level-( + 1) checkpoints is given by

n opt = N opt N opt +1 = λ λ +1 • C +1 C . ( 34
)
for all = 1, . . . , k -1.

Remarks. The optimal k-level pattern derived in this section has a rational number of segments, while the optimal integer solution could be much harder to compute. In Section 4, we use rounding to derive a practical solution. Still, Equation (31) provides a lower bound on the optimal overhead, which is met very closely in all our experimental scenarios.

Optimal subset of levels

The preceding section characterizes the optimal pattern by using k levels of checkpoints. In many practical cases, there is no obligation to use all available levels. This section addresses the problem of selecting the optimal subset of levels in order to minimize the overall execution overhead.

Checkpoint cost models

So far, we have assumed that all the checkpoint costs are fixed under a multi-level checkpointing scheme. In practice, the checkpoint costs may vary depending upon the implementation, and upon the subset of selected levels. In order to determine the optimal subset, we identify the following two checkpoint cost models:

• Fixed independent costs. The checkpoint cost C at level is the cost paid to save data at level , independently of the subset of levels used. In this model, the checkpoint costs stay the same for all possible subsets.

• Incremental costs. The checkpointing cost C at level is the additional cost paid to save data when going from level -1 to . In this model, the checkpoint cost at a particular level depends on the subset of levels selected.

For example, with k = 2 levels and C 1 = 10, C 2 = 20, two subsets are possible: {1, 2} and {2}. In the fixed independent cost model, these costs will stay unchanged regardless of the subset chosen. In the incremental cost model, since C 2 is the additional cost paid after C 1 is done, when using subset {2}, i.e., only placing level-2 checkpoints in the pattern, we need to adjust its cost as C 2 = 10 + 20 = 30. In both cases, once the subset is decided, the checkpoint costs at the selected levels can be computed and therefore considered as fixed constants. The theoretical analysis presented in Section 3.3 can then be used to compute the optimal pattern.

But how to determine the optimal subset of levels? Consider again the example with k = 2 levels. In the incremental cost model, Equation (31) suggests that the optimal solution (ignoring lower-order terms) uses both levels if and only if

2λ 1 C 1 + 2λ 2 C 2 ≤ 2 (λ 1 + λ 2 ) (C 1 + C 2 ) ⇔ 0 ≤ λ 1 C 2 -λ 2 C 1 2 ,
which is always true when assuming λ 1 ≥ λ 2 and C 1 ≤ C 2 . We can easily apply the same argument to show that the optimal subset must contain all levels available as long as all checkpoint costs are positive.

In the fixed independent cost model, however, it is not clear whether all available levels should be used. Consider the same example with k = 2 levels, and define α = λ 2 λ 1 and

β = C 2 C 1 .
The optimal solution uses both levels if and only if

2λ 1 C 1 + 2λ 2 C 2 ≤ 2 (λ 1 + λ 2 ) C 2 ⇔ 4αβ ≤ (β -1) 2 ,
which is not true when α = 0.5 and β = 2. In this case, using only level-2 checkpoints leads to a smaller overhead.

Dynamic programming algorithm

In the fixed independent cost model, the optimal subset of levels in a general k-level pattern could well depend on the checkpoint costs and error rates of different levels. One can enumerate all 2 k-1 possible subsets and select the one that leads to the smallest overhead.

The following theorem presents a more efficient dynamic programming algorithm when the number k of levels is large.

Theorem 3. Suppose there are k levels of checkpoints available and their costs are fixed. Then, the optimal subset of levels to use can be obtained by dynamic programming in O(k 2 ) time.

Proof. Let S opt (h) ⊆ {0, 1, . . . , h} denote the optimal subset of levels used by a pattern that is capable of handling errors up to level h, and let H opt (h) denote the corresponding optimal overhead (ignoring lower-order terms) incurred by the pattern. Define S opt (0) = ∅ and H opt (0) = 0. Recall that Λ [x,y] = y =x λ . We can compute H opt (h) using the following dynamic programming formulation:

H opt (h) = min 0≤ ≤h-1 H opt ( ) + 2Λ [ +1,h] C h , (35)
and the optimal subset is S opt (h) = S opt ( opt ) {h}, where opt is the value of that yields the minimum H opt (h).

The optimal subset of levels to handle all k levels of errors is then given by S opt (k) with the optimal overhead H opt (k). The complexity is clearly quadratic in the total number of levels.

Simulations

In this section, we conduct a set of simulations whose goal is threefold: (i) to verify the accuracy of the first-order approximation; (ii) to confirm the optimality of the subset of levels found by the dynamic programming algorithm; and (iii) to evaluate the performance of our approach and to compare it with other multi-level checkpointing algorithms. After introducing the simulation setup in Section 4.1, we proceed in two steps. First, in Section 4.2, we instantiate the model with realistic parameters from the literature and run simulations for all possible subsets of levels and roundings. Then, in Section 4.3, we instantiate the model with different test cases from the recent work of Di et al. [START_REF] Di | Optimization of multi-level checkpoint model for large scale HPC applications[END_REF][START_REF] Di | Toward an optimal online checkpoint solution under a two-level HPC checkpoint model[END_REF] on multilevel checkpointing and compare the overheads obtained with three approaches: (a) Young/Daly's classical formula; (b) our first-order approximation formula; and (c) Di et al.'s iterative/optimal algorithm. The simulator code is publicly available at http: // perso. ens-lyon. fr/ aurelien. cavelan/ multilevel. zip , so that interested readers can experiment with it and instantiate the model with parameters of their own choice.

Simulation setup

Checkpoint and recovery costs both depend on the volume of data to be saved, and are mostly determined by the hardware resource used at each level. As such, we assume that recovery cost for a given level is equivalent to the corresponding checkpointing cost, i.e., R = C for 1 ≤ ≤ k (unless specified otherwise). This a common assumption [START_REF] Moody | Design, Modeling, and Evaluation of a Scalable Multi-level Checkpointing System[END_REF][START_REF] Di | Optimization of multi-level checkpoint model for large scale HPC applications[END_REF], even though in practice the recovery cost can be expected to be somewhat smaller than the checkpoint cost [START_REF] Di | Optimization of multi-level checkpoint model for large scale HPC applications[END_REF][START_REF] Di | Toward an optimal online checkpoint solution under a two-level HPC checkpoint model[END_REF]. All costs are fixed and independent (as discussed in Section 3.4.1).

The simulator is fed with k levels of errors and their MTBFs µ = 1/λ , as well as the resilience parameters C and R . For each of the 2 k-1 possible subsets of levels (the last level is always included), we take the optimal pattern given in Theorem 2 and Corollary 1, and then try all possible roundings (floor and ceiling) based on the optimal (rational) number of checkpoints (n opt given in Equation ( 34)). For each rounding, we compare the following three overheads:

• Simulated overhead, obtained by running the simulation 10000 times and averaging the results;

• Corresponding theoretical overhead, obtained from Equations ( 19), ( 22) and (24) using the integer solution that corresponds to the rounding;

• Theoretical lower bound, obtained from Equation (31) with the optimal rational solution.

In the following, we associate Young/Daly's classical formula, defined as W opt = 2C Λ , with the highest checkpointing level available, i.e., C = C k . Note that in this case, Young/Daly's formula and Equation (29) can be used interchangeably, and the corresponding theoretical overhead is obtained with H opt = √ 2ΛC.

Assessing accuracy of first-order approximation

In this section, we run simulations with two sets of parameters, described in Table 1. For each set of parameters, we consider all possible subsets of levels. Then, for each subset, we compute the optimal pattern length and number of checkpoints to be used at each level. We show the accuracy of our approach in both scenarios, and we confirm the optimality of the subset of levels returned by the dynamic programming algorithm.

Using set of parameters (A)

The first set of parameters (shown in set (A) of Table 1) corresponds to the Coastal platform, a medium-sized HPC system of 1104 nodes at the Lawrence Livermore National Laboratory (LLNL). The Coastal platform has been used to evaluate the Scalable Checkpoint/Restart (SCR) library by Moody et al. [START_REF] Moody | Design, Modeling, and Evaluation of a Scalable Multi-level Checkpointing System[END_REF], who provided accurate measurements for the checkpoint costs using real applications (given in the first row of Table 1). There are k = 3 levels of checkpoints. First-level checkpoints are written to the local RAMs of the nodes, and this is the fastest method (0.5s). Second-level checkpoints are also written to local RAMs, but small sets of nodes collectively compute and store parity redundancy data, which takes a little while longer (4.5s). Lastly, Lustre is used to store third-level checkpoints onto the parallel file system, which takes significantly longer time (1051s). Failures were analyzed in [START_REF] Moody | Design, Modeling, and Evaluation of a Scalable Multi-level Checkpointing System[END_REF], and the error rates are given in the second row of Table 1. Note that the error rate at level 2 is higher than those of levels 1 and 3.

Results: Table 2 and Figure 3 present the simulation results. Table 2 shows, from left to right, the subset of levels used, the number of checkpoints computed by our first-order approximation formula for each possible rounding (N 1 , N 2 , N 3 ), the corresponding optimal pattern length (W opt (s)), the simulated overhead (Sim. Ov.), the corresponding theoretical overhead (Th. Ov.), the absolute and relative differences of these two overheads (Ab. Diff. = 100 × (Sim. Ov. -Th. Ov.), and Rel. Diff. = 100 × (Sim. Ov. -Th. Ov.)/Sim. Ov.), and finally the theoretical lower bound for this subset (Th. L.B.).

With k = 3, there are four possible subsets of levels, and both the best simulated overhead and the corresponding theoretical overhead are achieved for the subset {2, 3}, with N 2 = 35 and N 3 = 1 (highlighted in bold in the table). First, the difference between the simulated and theoretical overheads is very small, with a difference < 0.7% in absolute values, and a relative difference ranging from 2.9% (for subset {1, 2, 3}) to 8.14% (for subset {3}), which shows the accuracy of the first-order approximation for this set of parameters. The simulated overhead is always higher than the theoretical one, which is expected, because the first-order approximation ignores some lower-order terms. Next, we observe that, for each subset, all roundings of the number of checkpoints yield similar overheads on this platform, and the difference between the best and worst roundings is almost negligible.

Furthermore, using the best subset ({2, 3}) improves the overhead by over 50% compared to using level-3 checkpoints alone (as in Young/Daly's result). This is indeed the subset returned by the dynamic programming algorithm, and the result matches closely the minimum theoretical lower bound. Finally, comparing our result to the one obtained by the optimal two-level algorithm by Di et al. [START_REF] Di | Toward an optimal online checkpoint solution under a two-level HPC checkpoint model[END_REF] on this best subset, we see that the simulated overheads are similar under the optimal subset, as the patterns found using both approaches share the same number of checkpoints and the pattern lengths are also almost identical.

Using set of parameters (B)

The second set of parameters correspond to the execution of the LAMMPS application on the large BG/Q platform Mira at the Argonne National Laboratory (ANL) [START_REF] Balaprakash | Analysis of the tradeoffs between energy and run time for multilevel checkpointing[END_REF]. The parameters are presented in set (B) of Table 1. In this setting, the Fault Tolerance Interface (FTI) [START_REF] Bautista-Gomez | FTI: High performance fault tolerance interface for hybrid systems[END_REF] was used, which has four checkpoint levels (k = 4): Local checkpoint; Local checkpoint + Partnercopy; Local checkpoint + Reed-Solomon coding; and PFS-based checkpoint. The MTBFs correspond to the failure rates typically observed for petascale HPC applications [START_REF] Bautista-Gomez | FTI: High performance fault tolerance interface for hybrid systems[END_REF][START_REF] Moody | Design, Modeling, and Evaluation of a Scalable Multi-level Checkpointing System[END_REF][START_REF] Di | Optimization of multi-level checkpoint model for large scale HPC applications[END_REF].

Results: Table 3 and Figure 4 present the simulation results for this set of parameters. There are 8 possible subsets of levels. As before, we observe that the theoretical overhead is always slightly smaller than the simulated one, with an absolute difference of less than 2%, and a relative difference between 6-14%, demonstrating the accuracy of the model. Again, the results are very close to the theoretical lower bound. For this platform, the simulated overheads vary from 9.68% (with optimal subset of levels {1, 3, 4} found by the dynamic programming algorithm) to 14.3% (with level-4 checkpoints alone). For a given subset of levels, the rounding does not play a significant role, as W opt is also adjusted accordingly (increased or decreased) as a result of rounding. For instance, we observe that, for subset {1, 2, 3, 4}, the numbers of checkpoints at levels 1 and 2 are halved for the third rounding compared to the first rounding in Table 3, but W opt is also reduced by 31%, so that for the same amount of work, the number of checkpoints does not change by much. We can also see that the pattern length W opt for the smallest overhead is around 10400s, but only 2450s for the largest overhead. In fact, the largest pattern lengths are obtained for the highest cumulated checkpoint cost, which turns out to be 830s for {2, 3, 4} with N 2 = 16, N 3 = 4, N 4 = 1, and for {1, 2, 3, 4} with N 1 = 24, N 2 = 8, N 3 = 4 and N 4 = 1. This is because using more checkpoints both increases the error-free overhead and reduces the time lost due to re-executions upon errors. As a consequence, and to mitigate the aforementioned overhead, the length of the pattern increases (e.g., W opt = 17000s for {2, 3, 4} and W opt = 16600s for {1, 2, 3, 4}). And the converse is also true: when using fewer checkpoints, the error-free overhead decreases and the time lost upon errors increases. In order to compensate, the pattern length decreases (e.g., W opt = 8330s for {1, 2, 3, 4} with N 1 = 6, N 2 = 3, N 3 = 3 and N 4 = 1). We note that, in this case, our first-order solution slightly outperforms the iterative method by Di et al. [START_REF] Di | Optimization of multi-level checkpoint model for large scale HPC applications[END_REF] on multi-level checkpointing (with a simulated overhead of 9.68e-2 compared to 9.75e-2). The reason is that their algorithm computes a solution under the independent checkpointing model, i.e., checkpoints at different levels are taken according to different independent periods. However, it is not clear how such a model can be implemented in practice due to the difficulties as explained in Section 1 and the different options to rollback to a checkpoint in case of a fault. Therefore, we transformed their result to a pattern-based solution by rounding the different numbers of checkpoints obtained using their algorithm to create equal number of checkpoints at level -1 between two consecutive level-checkpoints. Although the best rounding is selected here for comparison, the result can still change drasti- cally the number of checkpoints computed by their initial rational solution without changing the pattern length, thus increasing the overhead.

Comparing performance of different approaches

In this section, we conduct simulations using settings from Di et al.'s recent work on multilevel checkpointing, which comprises two cases with four levels [START_REF] Di | Optimization of multi-level checkpoint model for large scale HPC applications[END_REF] and eight cases with two levels [START_REF] Di | Toward an optimal online checkpoint solution under a two-level HPC checkpoint model[END_REF], thus covering a wide range of configurations. For each case, we compare the performance of three different approaches: (a) Young/Daly's classical formula; (b) our first-order approximation formula; and (c) Di et al.'s iterative algorithm.

Using set of parameters (C)

We first run simulations for Cases #A and #B, whose parameters are presented in Table 4. These parameters are based on the FTI multilevel checkpointing model and have been used by Di et al. [START_REF] Di | Optimization of multi-level checkpoint model for large scale HPC applications[END_REF] to evaluate the performance of their approach. Note that the recovery cost is about half that of the checkpointing cost in Case #B.

In their work, Di et al. considered independent checkpointing periods, as opposed to the nested method based on periodic patterns (as discussed in Section 1). Although they provided an optimal solution, an iterative approach was used to compute it numerically in contrast to the simple formula we propose in this paper. Recall that using independent checkpointing periods allows checkpoints at different levels to be taken simultaneously, which can hardly be done in practice. Adapting their solution to our model results in rational numbers of checkpoints, and we again use rounding to resolve this issue. We find that, using the best roundings for both approaches, their solution turns out to be very similar to ours (with the same number of checkpoints, and close periods with < 1% difference).

Results: Figure 5 presents the overheads for both cases. First, we observe that Di et al.'s optimal iterative algorithm has almost identical performance to our solution, with a simulated overhead around 45% for Case #A and 140% for Case #B under both approaches. However, using Young/Daly's formula to checkpoint only at the highest level yields significantly worse overheads (around 90% for Case #A and 170% for Case #B). Overall, our solution is as good as Di et al.'s optimal numerical one (but has much less complexity), and it is up to 45% better than Young/Daly's formula in Case #A and 30% better in Case #B.

Note that the corresponding theoretical overhead of our solution is close to the simulated one for Case #A, but starts to diverge for Case #B. This is because first-order approximation is only accurate when the resilience parameters and pattern length are small compared to the MTBF, which is no longer true for Case #B. Specifically, we have:

• In Case #A, the optimal subset of levels is {2, 4}. The optimal pattern has length W opt = 1052s and consists of N 2 = 8 level-2 checkpoints followed by N 4 = 1 level-4 checkpoint, meaning that we have a level-2 checkpoint every 131.5s of computation. So a level-2 checkpoint is saved every 141.5s and a level-4 checkpoint is saved every 1222s. On the other hand, the combined MTBF for errors at levels 1 and 2 (handled by level-2 checkpoints) is 864s and the combined MTBF for errors at levels 3 and 4 (handled by level-4 checkpoints) is 6171s. Hence, we have 141.5 864 = 0.164 and 1222 6171 = 0.198, which are reasonably small, making our solution accurate.

• In Case #B, the optimal subset of levels is {1, 4}, and the optimal pattern has W opt = 223s, N 1 = 5 and N 4 = 1. Thus, we have a level-1 checkpoint every 44.6s of computation. So a level-1 checkpoint is saved every 45.6s and a level-4 checkpoint is saved every 298s.

The MTBF for errors at level 1 is 864s and the combined MTBF for errors at levels 2, 3 and 4 (handled by level-4 checkpoints) is 360s. Thus, we have 44.6 864 = 0.052, which is fine, but 298 360 = 0.828, which is too high and essentially makes the first-order solution inaccurate. Despite the difference between the theoretical and simulated overheads under Case #B, the proximity of our solution to Di et al.'s optimal numerical solution nevertheless shows the usefulness of first-order approximation for determining the optimal multi-level checkpointing patterns.

Using set of parameters (D)

Finally, we run simulations for eight cases, whose parameters are presented in Table 5. These parameters have been used by Di et al. [START_REF] Di | Toward an optimal online checkpoint solution under a two-level HPC checkpoint model[END_REF] to evaluate their two-level checkpointing model, and as such, each case consists of only two checkpointing levels. In their work, the authors proposed an optimal solution by solving complex mathematical equations using numerical method. Again, for each case, we compare the simulated overheads obtained with the three different approaches.

In this set of parameters, the MTBF has a large variation, ranging from more than 1 hour (Case 1) to less than 4 minutes (Case 8). Similarly, the checkpointing costs vary from 10s (Cases 4 and 5) to 300s (Case 8). Note that Cases 7 and 8 have both very short MTBFs and very high checkpointing costs, resulting in a lot of errors and recoveries. In particular, the checkpointing cost at level 2 in Case 8 (300s) is larger than the MTBF at level 1 (216s).

Results: Figure 6 presents the simulation results for the eight cases. First, we observe that the optimal algorithm by Di et al. only yields a slightly better simulated overhead compared to our simple first-order approximation solution (by less than 2% in Cases 1 to 6). However, our solution always improves significantly over Young/Daly's formula, from 2% (Case 1) up to 100% (Case 6). Due to their short MTBFs, Cases 7 and 8 stand out and incur much higher overheads compared to the first six cases (thus their results are presented in a separate plot). Still, considering Case 8, we are able to improve over Young/Daly's solution by as much as 2500% (in absolute value of the overhead), and we are off the optimal simulated overhead by only 300%. In addition, Figure 6 shows the theoretical overheads obtained both with our formula and the solution provided by Di et al. in [START_REF] Di | Toward an optimal online checkpoint solution under a two-level HPC checkpoint model[END_REF]. As expected, our first-order approximation remains accurate when the MTBF is large, as in Cases 1, 2 and 4. However, it becomes less accurate with shorter MTBFs and higher error rates, especially in Cases 7 and 8 (which do not represent healthy HPC platforms).

Summary of results

From the simulation results, we conclude that first-order approximation remains a valuable performance model for evaluating checkpointing solutions in HPC systems (as long as the error rates stay reasonably low). We have demonstrated, through an extensive set of simulations with a wide range of parameters, the usefulness of multi-level checkpointing (over using only one level of checkpoints) with significantly reduced overheads. The results also corroborate the analytical study by showing the benefit of selecting an optimal subset of levels among all the levels available. Overall, our approach achieves the optimal or near-optimal performance in almost all cases, except when the MTBF is too small, in which case even the optimal solution yields an unacceptably high overhead (e.g., Case 8 of Table 5).

Conclusion

This paper has studied multi-level checkpointing protocols, where different levels of checkpoints can be set; lower levels deal with frequent errors that can be recovered at low cost (for instance with a memory copy), while higher levels allow us to recover from all errors, such as node failures (for instance with a copy in stable storage). We consider a general scenario with k levels of faults, and we provide explicit formulas to characterize the optimal checkpointing pattern, up to first-order approximation. The overhead turns out to be of the order of k =1 √ 2λ C , which elegantly extends Young/Daly's classical formula. The first-order approximation to the optimal k-level checkpointing pattern uses rational numbers of checkpoints, and we prove that all segments should have equal lengths. We corroborate the theoretical study by an extensive set of simulations, demonstrating that greedily rounding the rational values leads to an overhead very close to the lower bound. Furthermore, we provide a dynamic programming algorithm to determine those levels that should be selected, and the simulations confirm the optimality of the subset of levels returned by the dynamic programming algorithm.

The problem of finding a first-order optimal pattern with an integer number of segments to minimize the overhead remains open. It may well be the case that such an integer pattern is not periodic at each level and uses different-length segments. However, the good news is that the rounding of the rational solution provided in this paper seems quite efficient in practice.

Figure 1 :

 1 Figure 1: Checkpointing patterns (highlighted using red bars) with (a) k = 3, (b) k = 1, and (c) k = 2 levels.

Figure 2 :

 2 Figure 2: Example of a 4-level pattern. Here, we let c = C1|C2| • • • |C denote the succession of checkpoints from level 1 to level .

Figure 4 :

 4 Figure4: Simulated and (corresponding) theoretical overheads for all possible subsets of levels with the best and worst roundings for each subset using set of parameters (B).

Figure 5 :

 5 Figure 5: Performance comparison of the three different approaches using two cases from Di et al. [7].

Figure 6 :

 6 Figure 6: Performance comparison of the three different approaches using 8 cases from Di et al. [8].

Table 1 :

 1 Sets of parameters (A) and (B) used as inputs for simulations.

	Set	From	Level	1	2	3	4
	(A)	Moody et al. [14] MTBF (s) 5.00e6 5.56e5 2.50e6 C (s) 0.5 4.5 1051	--
	(B)	Balaprakash et al. [1] MTBF (s) 3.60e4 7.20e4 1.44e5 7.20e5 C (s) 10 30 50 150

Table 2 :

 2 Simulation results using set of parameters (A).Levels N 1 N 2 N 3 W opt (s) Sim. Ov. Th. Ov. Abs. Diff. Rel. Diff. Th. L.B.

	{3} --1 2.96e4	7.74e-2 7.11e-2	0.63%	8.14%	7.11e-2
	{1,3}	14 -1 3.09e4 13 -1 3.09e4	7.40e-2 6.85e-2 7.39e-2 6.85e-2	0.55% 0.54%	7.43% 7.31%	6.85e-2
	{2,3}	-35 1 7.27e4 3.44e-2 3.33e-2 0.11% -34 1 7.25e4 3.46e-2 3.33e-2 0.13%	3.20% 3.33e-2 3.76%
	{1,2,3}	33 33 1 7.27e4 32 32 1 7.24e4	3.46e-2 3.35e-2 3.45e-2 3.35e-2	0.11% 0.10%	3.18% 2.90%	3.35e-2

  Simulated and (corresponding) theoretical overheads for all possible subsets of levels with the best and worst roundings for each subset using set of parameters (A).

		0.12			Sim. Overhead (Young/Daly) Sim. Overhead (Di et al.)
		0.10			Theoretical L.B. (Optimal Subset) Sim. Overhead (Best Rounding)
					Corresp. Theoretical Overhead
					Sim. Overhead (Worst Rounding)
		0.08			Corresp. Theoretical Overhead
	Overhead	0.06
		0.04
		0.02
		0.00	{3}	{1,3}	{2,3}	{1,2,3}
				Subsets of Levels
	0.00 0.02 0.04 0.06 0.08 0.10 0.14 0.16 0.18 Figure 3: {4} Overhead 0.12	{1,4} {2,4} {3,4} {1,2,4} {1,3,4} {2,3,4} {1,2,3,4} Sim. Overhead (Young/Daly) Sim. Overhead (Di et al.) Theoretical L.B. (Optimal Subset) Sim. Overhead (Best Rounding) Corresp. Theoretical Overhead Sim. Overhead (Worst Rounding) Corresp. Theoretical Overhead
				Subsets of Levels

Table 3 :

 3 Simulation results using set of parameters (B).Levels N 1 N 2 N 3 N 4 W opt (s) Sim. Ov. Th. Ov. Abs. Diff. Rel. Diff. Th. L.B.

	{4} ---1 2.45e3 1.43e-1 1.22e-1	1.9%	13.3%	1.22e-1
	{1,4}	5 --1 3.79e3 1.18e-1 1.05e-1 4 --1 3.61e3 1.18e-1 1.05e-1	1.3% 1.3%	11.0% 11.0%	1.05e-1
	{2,4} -5 -1 6.00e3 1.11e-1 1.00e-1	1.1%	9.9%	1.00e-1
	{3,4}	--11 1 1.55e4 9.96e-2 9.02e-2 --10 1 1.44e4 9.91e-2 9.01e-2	0.94% 0.90%	9.44% 9.08%	9.01e-2
		9 3 -1 6.41e3 1.11e-1 1.03e-1	0.8%	7.2%
	{1,2,4}	6 2 -1 5.21e3 1.13e-1 1.04e-1 6 3 -1 5.84e3 1.11e-1 1.03e-1	0.9% 0.8%	8.0% 7.2%	1.02e-1
		4 2 -1 4.74e3 1.17e-1 1.05e-1	1.2%	10.3%
		21 -7 1 1.58e4 9.72e-2 8.99e-2	0.73%	7.51%
	{1,3,4}	18 -6 1 1.40e4 9.82e-2 8.98e-2 0.84% 14 -7 1 1.04e4 9.68e-2 9.01e-2 0.67%	8.55% 6.92%	8.96e-2
		12 -6 1 1.26e4 9.85e-2 9.04e-2	0.81%	8.22%
		-16 4 1 1.70e4 1.07e-1 9.75e-2	0.95%	8.9%
	{2,3,4}	-12 3 1 1.36e4 1.04e-1 9.73e-2 -12 4 1 1.47e4 1.05e-1 9.68e-2	0.67% 0.82%	6.4% 7.8%	9.68e-2
		-9 3 1 1.17e4 1.05e-1 9.75e-2	0.75%	7.1%
		24 8 4 1 1.66e4 1.09e-1 1.00e-1	0.9%	8.2%
		18 6 3 1 1.32e4 1.08e-1 9.99e-2	0.81%	7.5%
		12 4 4 1 1.15e4 1.11e-1 1.03e-1	0.8%	7.2%
	{1,2,3,4}	9 3 3 1 9.17e3 1.14e-1 1.05e-1 16 8 4 1 1.51e4 1.08e-1 9.95e-2	0.9% 0.85%	7.9% 7.9%	9.92e-2
		12 6 3 1 1.20e4 1.09e-1 1.00e-1	0.9%	8.3%
		8 4 4 1 1.05e4 1.16e-1 1.05e-1	1.1%	9.5%
		6 3 3 1 8.33e3 1.19e-1 1.08e-1	1.1%	9.2%

Table 4 :

 4 Set of parameters (C) used as input for simulations.

				Set (C), from Di et al. [7]
				Level	1	2	3	4
				C (s)	8	10 80	90
			Case #A	R (s)	8	10 80	90
				MTBF (s) 2160 1440 8640 21600
				C (s)	1	20 60	70
			Case #B	R (s)	1	10 30	35
				MTBF (s) 864 864 1080 1440
		2.0
			Young/Daly's formula (Simulated)
			Our first-order approx. (Simulated)
			Our first-order approx. (Theoretical)
			Di et al.'s algorithm (Simulated)
		1.5
	Overhead	1.0
		0.5
		0.0	Case #A		Case #B

Table 5 :

 5 Set of parameters (D) used as input for simulations.

			Set (D), from Di et al. [8]
		Level	1	2	Level	1	2
	Case 1	C (s) MTBF (s) 3600 21600 20 50 Case 5	C (s) MTBF (s) 432 2160 10 40
	Case 2	C (s) MTBF (s) 1728 8640 20 50 Case 6	C (s) MTBF (s) 432 2160 100 20
	Case 3	C (s) MTBF (s) 864 4320 20 100 Case 7	C (s) MTBF (s) 288 1440 40 200
	Case 4	C (s) MTBF (s) 864 4320 10 40 Case 8	C (s) MTBF (s) 216 1440 50 300

We use the terms fault, failure and error indifferently.RR n°8930

RR n°8930

In practice, failures from different checkpointing levels can exhibit potential correlation[START_REF] Heien | Modeling and tolerating heterogeneous failures in large parallel systems[END_REF][START_REF] Di | Optimization of multi-level checkpoint model for large scale HPC applications[END_REF]. Consideration of correlated failures is beyond the scope of this paper.

The original optimal solution of[START_REF] Di | Toward an optimal online checkpoint solution under a two-level HPC checkpoint model[END_REF] considers faults in checkpointing but not during recoveries. We adapt its solution to exclude faults in checkpointing so to be consistent with the model in this paper for a fair comparison. The results reported herein are based on this modified solution.RR n°8930

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgment

We would like to thank the anonymous reviewers for their valuable comments, which have greatly improved the quality of this paper. This research was funded in part by the European project SCoRPiO, by the LABEX MILYON (ANR-10-LABX-0070) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French

National Research Agency (ANR), by the PIA ELCI project, and by the ANR Rescue project. Yves Robert is with Institut Universitaire de France.