Nisrine Jrad

Marco Congedo

Identification of spatial and temporal features of EEG

Keywords: Brain computer interface (BCI), Electroencephalography (EEG), Global Field Power (GFP), Support Vector Machine (SVM), Event related Potentials (ERP), spatial and temporal feature extraction

Identification of spatial and temporal features of EEG

Nisrine Jrad, Marco

Introduction

Brain Computer Interface (BCI) is a tool for direct communication between human and external devices [START_REF] Dornhege | Towards Brain-Computer interfacing[END_REF]. Many BCIs make use of Electroencephalography (EEG) signals to categorize two or more classes and associate them to simple computer commands. EEG is a popular technology for BCI purposes due to its low-cost, portability and non-invasiveness. Classification of EEG is not an easy task because this latter can only measure blurred cortical activities due to the diffusion of the skull and the skin [START_REF] Nunez | Electric Field of the Brain[END_REF]. Furthermore, recordings are highly contaminated by noise of various sources. Interestingly, EEG signals often reveal various spatial and temporal characteristics. Thus, it is important to characterize both spatial and temporal dynamics of EEG data to provide reliable BCI control.

Usually spatial spatial decomposition is performed to extract spatial features [START_REF] Blankertz | Single-trial analysis and classification of ERP components -a tutorial[END_REF][START_REF] Müller-Gerking | Designing optimal spatial filters for single-trial EEG classification in a movement task[END_REF][START_REF] Gouy-Pailler | Nonstationary brain source separation for multiclass motor imagery[END_REF][START_REF] Rivet | xDAWN Algorithm to Enhance Evoked Potentials: Application to Brain Computer Interface[END_REF]. Signal-processing criteria like the Signal-to-Noise Ratio and ratio of class variances [START_REF] Blankertz | Single-trial analysis and classification of ERP components -a tutorial[END_REF][START_REF] Rivet | xDAWN Algorithm to Enhance Evoked Potentials: Application to Brain Computer Interface[END_REF][START_REF] Hoffmann | Spatial filters for the classification of event-related potentials[END_REF] have been often employed for Evoked Response Potential (ERP). These methods are based on the assumption of the instantaneous and approximately linear relation between the amplitudes of the generating cerebral electrophysiological current sources and the amplitude of the observed scalp potential field. The idea is to find a linear transformation of the data optimizing the extraction of the relevant EEG spatial feature and the noise suppression. The performance of such filters mainly depends on the accuracy of spatial covariance estimations and is jeopardized by the non-stationary nature of the noise. Although a relation might be found between the objective functions of these aforementioned criteria and class separability, this relation has, to the best of our knowledge, never been addressed explicitly. Hence, spatial filter can be suboptimal from the classification point of view. Another limitation of these methods is that spatial decomposition is performed in a specific time window. However, ERPs reflect several temporal components, thus, spatial decomposition should be performed for each interesting interval occurring in the pre-fixed window.

Recently, some algorithms have been proposed to study where the discriminative information lies into the spatio-temporal plane. They visualize a matrix of separability measures into the spatio-temporal plane of the experimental conditions. The matrix is obtained by computing a separability index for each pair of spatial electrode measurement and time sample. Several measures of separability have been used, for instance the signed-r 2 -values [START_REF] Blankertz | Single-trial analysis and classification of ERP components -a tutorial[END_REF], Fisher score and Student's t-statistic [START_REF] Müller | Machine learning techniques for brain-computer interfaces[END_REF], or the area under the ROC curve [START_REF] Green | Signal detection theory and psychophysics[END_REF]. Separability matrix should be sought as to automatically determine intervals with fairly constant spatial patterns and high separability values. This proves difficult and heuristic techniques are often employed to approximate interval borders. In addition, the three first aforementioned measures rely on the assumption that class distributions are Gaussian, which is seldom verified.

To overcome all these drawbacks, we develop a spatio-temporal data driven decomposition technique. A two-stage feature extraction technique is proposed. First, a time feature extraction is performed based on Global Field Power (GFP) [START_REF] Lehmann | Reference-free identification of components of checkerboard-evoked multichannel potential fields[END_REF], defined for each time sample as the sum of the square potential across electrodes. Second, a spatially weighted SVM (sw-SVM) is used to learn for each time interval a sparse spatial filter optimizing directly the classification performance [START_REF] Jrad | sw-SVM : sensor weighting Support Vector Machines for EEG-based Brain-Computer Interfaces[END_REF]. Finally, the ensemble of sw-SVMs obtained on selected temporal features are combined using a weighted average, to get a robust decision function. The remainder of this paper is organized as follows. The proposed method is introduced in Section 2. Section 3 accounts for experiment and data set description. Section 4 shows and discusses experimental results. Finally, Section 5 holds our conclusions.

Method

Background brain activities, irrelevant to BCI tasks, continuously generate EEG signals that can be recorded almost anywhere over the scalp. These signals interfere with the EEG signals triggered by stimuli. Interestingly, post-stimulus signals are generated in particular regions of the brain at a given interval of time. Hence, they present specific spatial and temporal characteristics. This section describes a method for analyzing ERPs spatial and temporal characteristics considering the sequence :

• a temporal feature extraction,

• a spatial feature extraction embedded in a classification scheme,

• an ensemble of classifiers learning technique.

Problem description

In the following, we consider BCI applications with two classes of action providing a training set of labeled trials from which a decision function is learned. The decision function should correctly classify unlabeled trials. Let us consider an EEG post-stimulus trial recorded over S electrodes in a short time period of T samples as a matrix Xp ∈ R S×T . Hence, the entire available set of data can be denoted {( X1 , y 1 ), ..., ( Xp , y p ), ..., ( XP , y P )} with y p ∈ {-1, 1} the class labels. Our task consists on finding the spatio-temporal features that maximize discrimination between two classes.

Temporal features

To select temporal intervals in the ERP where discriminative peaks appear, Global Field Power (GFP) [START_REF] Lehmann | Reference-free identification of components of checkerboard-evoked multichannel potential fields[END_REF] is computed on the difference of the grand averages of the two class post-stimulus trials as following:

(GF P 2 ) t = 1 S S s=1     p|yp=1 ( Xp ) s,t - p|yp=-1 ( Xp ) s,t     2 , (1)
(GF P 2 ) t is the square of Global Field Power at a given time sample t and ( Xp ) s,t is the EEG recording of trial p at electrode s and p th post-stimuls time sample t.

Pronounced deflections with large peaks, denoting big dissimilarities between the two activities, are associated with large GFP values. Windows involving significant temporal features are chosen as intervals where GFP is high relative to the background EEG activity.

To select significant windows we require a statistical threshold for the observed GFP of the difference grand average trials in the two classes. Such threshold is estimated with a resampling method as the 95 th percentile (5% type I error rate) of the appropriate empirical null distribution. For Q and R observed single trials in classes labeled 1 and -1, respectively, we resample Q and R trials with random onset from the entire EEG recording. We compute the difference of the grand average of the Q and R random trials and retain the maximum value of GFP. The procedure is repeated 1000 times and the sought threshold is the 95 th percentile of such max-GPF null distribution after 10% trimming. The trimming makes the estimated GFP more robust with respect to outliers given by eye blinks and other large-amplitude artefacts. Taking the max-GPF at each resampling ensures that the nominal type I error rate is preserved regardless the number of windows that will be declared significant.

Noteworthily, contiguous samples with high GFP coincide with stable deflection configurations where spatial characteristics of the field remains unchanged [START_REF] Lehmann | Reference-free identification of components of checkerboard-evoked multichannel potential fields[END_REF]. Since within each selected time window the spatial pattern is fairly constant, average across time is calculated. Averaging over time rules out aberrant values, reduces signal variability and attenuates noise. Besides, it reduces dramatically time dimensionality to I where I is the number of significant time features.

Spatial features and classifier : sw-SVM

Temporal filtering provides us with X p ∈ R S×I trials. Each column vector (x p ) i ∈ R S reflects a spatial characteristic at a temporal feature i ∈ {1, ..., I}. Hereafter, x p will refer to (x p ) i for convenience. Hence, I spatial filters are learned over the different time components.

In this work, spatial filtering is learned jointly with a classifier in the theoretical framework of SVM. This method was proposed in [START_REF] Jrad | sw-SVM : sensor weighting Support Vector Machines for EEG-based Brain-Computer Interfaces[END_REF] and called sw-SVM for spatially weighted SVM. It has the advantage of learning a spatial filter so as to improve separability of classes whilst reducing classification errors. It involves spatial feature weights in the primal SVM optimization problem and tunes these weights as hyper-parameters of SVM. We denote by

d = [d 1 , d 2 , ..., d S ] T ∈ R S the
where w ∈ R d×1 is the normal vector, b ∈ R is an offset, ξ p are slack variables that ensure a solution in case data are not linearly-separable, and C is the regularization parameter that controls the trade-off between a low training error and a large margin. By setting to zero the derivatives of the partial associated Lagrangian according to the primal variables w, b and ξ p the optimization problem of the dual formulation can be written as :

min D J(D) subject to S s=1 d 2 s = 1. (3)
with

J(D) =        max α 1 T α -1 2 α T Y T X T D T DXY α subject to y T α = 0 and 0 ≤ α p ≤ C ∀p ∈ {1, . . . , P },
with α is the vector of Lagrangian multipliers, X = {x 1 , ..., x P }, y T = {y 1 , ..., y P } and Y = Diag(y). The optimization problem, as stated in (3), is inherently non-convex since the boundary of the unit ball given by { S s=1 d 2 s = 1} is not a convex set. To remove this non-convexity, we perform a variable substitution D = D T D. We arrive at the following optimization problem which is convex in D:

min D J( D) subject to S s=1 ds = 1. (4)
This is clearly an instance of the Multiple Kernel Learning (MKL) problem proposed in [START_REF] Rakotomamonjy | SimpleMKL[END_REF] where one homogeneous degree one polynomial kernel is used over electrode s samples and ds is its corresponding positive mixing coefficient. Authors of [START_REF] Rakotomamonjy | SimpleMKL[END_REF] prove that the search for the optimal D is convex, yielding fast convergence toward the optimal conditional solution. Hence, the optimization problem can be solved efficiently using a gradient descent as in SimpleMKL [START_REF] Rakotomamonjy | SimpleMKL[END_REF]. For effect of the ℓ 1 constraints over ds , the sought spatial filters will be sparse. Linear sw-SVM can be extended to a non-linear sw-SVM by replacing inner products with a suitable kernel.

To solve this sw-SVM optimization problem, we proceed by alternating the search for α and D. For D fixed, the problem is reduced to a ℓ 1 soft margin SVM with the only difference being that x p is replaced by Dx p in the inequality constraint. Primal and dual objective functions of SVM are convex, and their solution is obtained by any of the available SVM solver [START_REF] Canu | SVM and Kernel Methods Matlab Toolbox[END_REF]. For the so obtained α, the minimization problem over D is smooth and convex [START_REF] Lemaréchal | Practical aspects of the moreau-yosida regularization: Theoretical preliminaries[END_REF]. Hence, we can use a reduced gradient method which converges for such functions [START_REF] Luenberger | Linear and Nonlinear Programming[END_REF]. Once the gradient of J( D) is computed, D is updated by using a descent direction ensuring that the equality constraint and the non-negativity constraints on { ds } are satisfied. These two steps are iterated until a stopping criterion is reached. The stopping criterion we chose is based on a norm variation of the sensor weights.

Ensemble of sw-SVM classifiers

As seen above, a way to reduce EEG variability is to perform signal averaging across time. Another way to reduce this influence, from a classification point of view, is to use an ensemble of classifiers. This technique has proven efficient in [START_REF] Grandvalet | Bagging equalizes influence[END_REF][START_REF] Rakotomamonjy | BCI Competition III : Dataset II -Ensemble of SVMs for BCI P300 speller[END_REF] where a set of classifiers were learned on different subsets of a given training set. Classifier outputs were then averaged. In this work, we propose to learn an ensemble of sw-SVMs, each one on a given temporal feature and to process a weighted average on sw-SVM outputs. Weights are used to keep significant classifiers and discard insignificant ones.

Weights are set as the product of two functions growing proportionally with the accuracies of the two (positive and negative) classes, evaluated on a validation set. If T P , F P , T N and F N hold for True Positive, False Positive, True Negative and False Negative of a given sw-SVM classifier, respectively, then Positive Predective Value (P P V ) and Negative Predective Value (NP V ) can be defined as :

P P V = T P T P + F P , NP V = T N F N + T N ,
and weights associated to the classifier are such as:

0 if P P V < 0.5 or NP V < 0.5 tan(2(P P V -0.5)) × tan(2(NP V -0.5)) otherwise (5)
This weighting strategy is ideal for unbalanced data sets since it privileges classifiers that jointly present good accuracies for both classes.

Experiment

The proposed method was evaluated on a visual feedback Error Related Potentials (ErrP) [START_REF] Miltner | Event-related brain potentials following incorrect feedback in a time-estimation task: Evidence for a generic neural system for error detection[END_REF] experiment. Eight healthy volunteers (including three women) participated in this experiment. All subjects were BCI-naive at the time of the experiment. Subjects had to retain the position of an ensemble of two to nine digits. The digits were displayed as a sequence in square boxes and evenly distributed along a circle. When the sequence disappeared, a target digit was shown and subjects were asked to click on the box where it previously appeared. A visual feedback indicates whether the answer was wrong or correct.

The experiment involved two sessions that lasted together approximately half an hour. Each session consisted of six blocks of six trials, for a total of 6 × 6 × 2 = 72 trials.

The temporal order of each trial, illustrated in Figure 1, is detailed next. The score, initially zero, was displayed for 3000 ms followed by a fixation cross, which was in turn displayed for 3000 ms. Then the memorization sequence started with variable duration depending on the number of digits the subject had to memorize. When it ended the subject was asked to click on the box where the target digit had appeared. Once the subject had answered the interface was paused for 1500 ms and then turned the clicked box into green upon a correct answer or into red upon an erroneous answer. This feedback lasted for 2500 ms. The 1500 ms preceding the feedback was introduced to avoid any contamination of ErrP by beta rebound motor phenomena linked to mouse clicking [START_REF] Jurkiewicz | Post-movement beta rebound is generated in motor cortex: Evidence from neuromagnetic recordings[END_REF]. The subject was then asked to report if the feedback (correct/error) matched his expectation by a mouse click ("yes"/"no"). Following his answer a random break of 1000 ms to 1500 ms preceded the beginning of the new trial. The number of digits was adapted with an algorithm tuned to allow about 20% errors for all subjects. The mean error rate (standard deviation) was equal to 17.87(±4.64)% of the trials. Recordings of the EEG were made using 31 sensors from the extended 10/20 system. The Mitsar 202 DC EEG amplifier was used for signal acquisition. We used one ground positioned on the forehead and two earlobe references. The acquisition software links the two references digitally and computes one common reference for all cephalic sensors. During acquisition, EEG was band-pass filtered in the range 0.1 -70 Hz and digitized at 500 Hz.

Experimental results

Preprocessing

Raw EEG potentials were first re-referenced to the common average by subtracting from each sensor the average potential (over the 31 sensors) for each time sample. Many studies report two peaks, Ne and Pe, as the main components of Error Related Potential components [START_REF] Miltner | Event-related brain potentials following incorrect feedback in a time-estimation task: Evidence for a generic neural system for error detection[END_REF]. Ne shows up about 250 ms after the response as a sharp negative peak and Pe shows up about 300 to 500 ms after the response as a broader positive peak. According to these knowledge, only a window of 1000 ms posterior to the stimulus has been considered for each trial, which results in 500 samples per sensor. A 1 -10 Hz 4 th order Butterworth filter was applied as error related potentials are known to be a relatively slow cortical potential. No artifact rejection algorithm was applied and all trials were kept for analysis.

Cross Validation

Single trial classification of ErrPs was assessed using a 5-Cross Validation technique. Because of the limited number of trials, temporal and spatial extraction were performed on four subsets. sw-SVM was trained with a given regularization parameter C (selected from a set of pre-defined values of C) and weighted using [START_REF] Gouy-Pailler | Nonstationary brain source separation for multiclass motor imagery[END_REF]. Performances were computed on the remaining subset. For each value of C, this process was repeated five times for a given subject and averaged. The highest average accuracy was reported.

Results

Figure 2 shows the average of the difference error-minus-correct for channel FCz of subject 7. It also reports GFP computed on the difference average. Five components are to be noted. A negative deflection can be seen around 240ms after the feedback and a second positive component occurs about 350ms. Three more peaks are also detected; a negative deflection around 500ms, a less pronounced negative deflection around 700ms and a small positive deflection around 800ms. Scalp potentials topographies associated with the five extracted temporal features are also shown on Figure 2. The first negative peak seems to be occipital whereas the second positive peak covers a rather fronto-central area. The third peak covers a parieto-central area, the fourth peak covers the whole right hemisphere and the last one is more central.

Figure 2 shows accuracies for error and correct classes for each sw-SVM and their corresponding weights (normalized between 0 and 1). Only sw-SVMs learned on the most pronounced peaks (second and third) show good accuracies in both classes and are thus retained.

An important question is whether time interval selection, found by GFP, are consistent across different partitions of the data and across subjects. Figure 3 shows, for each subject, and each of the five partitions, temporal intervals (in white) selected on the ErrP data set. Because of the very small number of trials used in each partition, some inter-partition differences can be noted in these data, but overall, the procedure appears robust and meaningfull. Latencies, thus selected time intervals, are different from subject to subject, which is not surprising. However, for almost all subjects, an important activity is noted between 400 and 600ms. These findings confirm those of [START_REF] Miltner | Event-related brain potentials following incorrect feedback in a time-estimation task: Evidence for a generic neural system for error detection[END_REF] where a negative deflection, following an incorrect visual feedback of a time-production task, peaked at 330ms with a duration of 260ms. This witnesses in favor of the effectiveness and the consistency of the proposed temporal feature extraction.

Concerning classification results, figure 4 shows the 5-Cross Validation performance provided by a classical SVM approach where all electrodes were used, the sw-SVM where only one spatial filter was used on the whole trial duration and the proposed method where spatio-temporal features were extracted. The proposed method proved constantly superior to SVM and sw-SVM. A paired Wilcoxon signed-rank test was evaluated to compare the proposed method to SVM and sw-SVM and p-values of 0.0078 and 0.0391 were obtained. We conclude that, inclusion of temporal features along with learning an ensemble of classifiers, provide with superior performance.

Conclusion

Spatio-temporal feature identification was addressed. An analysis of Global Field Power highlighted time periods of interest where effects are likely to be the most robust yielding to a data-driven temporal feature extraction. For each temporal feature, a spatial filter was learned jointly with a classifier in the SVM theoretical framework. Spatial filters were learned to optimize classification performance. A weighted averaging on the so obtained ensemble of classifiers yielded to a robust final decision function. Experimental results on Error-related Potentials illustrate the efficiency of the method from a physiological and a machine learning points of view. Further research may extract all relevant aspects of brain post-stimulus dynamics recorded in EEG (spatio-temporal-frequential).

  spatial filter and D a matrix with d on the diagonal. Matrix D is learned by solving the sw-SVM optimization problem: min w,b,ξ,D 1 2 w 2 + C P p=1 ξ p subject to y p ( w, Dx p + b) ≥ 1 -ξ p and ξ p ≥ 0 ∀p ∈ {1, . . . , P } and S s=1 d 2 s = 1 ∀s ∈ {1, . . . , S},

Figure 1 :

 1 Figure 1: Temporal diagram of one ErrP trial.

Figure 2 :

 2 Figure2: Top : GFP computed on the difference of the grand average error-minus-correct for 1s trials, selected intervals and topographies associated. Middle : the difference computed on electrode FCz. Bottom : accuracies for error (blue bar) and correct (green bar) classes and sw-SVM associated weights (red bar, normalized between 0 and 1).

Figure 3 :

 3 Figure 3: Selected time intervals are shown in white pixels, for each of the 8 subjects and 5 partitions. Each matrix refers to a subject where columns hold time-course and rows hold partitions.

Figure 4 :

 4 Figure 4: From left to right: performances of classical SVM, sw-SVM and the proposed method for the 8 subjects. Mean (std) accuracies across the 8 subjects are of 70.71(10.77), 80.71(6.61) and 87.80(3.63) respectively.

Acknowledgment

This work has been supported through the project OpenViBE2 of the ANR (National Research Agency), France. Authors would like to acknowledge Sandra Rousseau from GIPSA-Lab Grenoble INP for providing the data.


O Beaumont
email: olivier.beaumont@labri.fr

E M Daoudi
email: 2mdaoudi@sciences.univ-oujda.ac.ma

N Maillard
email: 3nmaillard@inf.ufrgs.br

P Manneback
email: pierre.manneback@fpms.ac.be

J.-L Roch

Tradeoff to minimize extra-computations and stopping criterion tests for parallel iterative schemes

.

Introduction

In general, iterative schemes consist in re-iterating some computation until global convergence is reached (i.e. until the stopping test criterion is fulfilled). Standard algorithm performs the stopping test after each iteration, so that each step of the algorithm can be decomposed into two phases: i. computation of the "body" of the iteration and ii. control of the stopping test. Let n iter denote the exact number of iterations obtained by the standard algorithm. The execution time T standard of the standard algorithm is given by T standard = n iter (t comp + t synch ), where t comp denotes the execution time of one iteration and t synch denotes the cost of evaluating a stopping test.

On a parallel architecture, the evaluation of the stopping criterion test is often achieved using global synchronization. Unfortunately, for a parallel implementation on a distributed memory architecture (typically a cluster of processors or a computational grid), it is well known that such a global synchronization is extremely expensive since it requires a global waiting time that lower performances. Indeed, each processor must wait for the contribution of the latest processor in order to continue its work. In order to mitigate the problem of global synchronizations, different approaches have been proposed and studied in the literature.

A first approach consists in gathering all global synchronizations of one iteration in order to perform them together. This technique is used in [START_REF] Baserman | Preconditioned CG methods for sparse matrices on massively parallel machines[END_REF] for Conjugate Gradient method, but it does not reduce the number of stopping criterion tests and consequently the number of global synchronizations for a parallel implementation.

A second approach, which has been widely studied in the literature, consists in desynchronizing communications and iterations, resulting in an asynchronous algorithm [START_REF] Bertsekas | Parallel and distributed computation, numerical methods[END_REF][START_REF] Baz | An efficient termination method for asynchronous iterative algorithms on message passing architectures -LAAS Report[END_REF]. Unfortunately, the semantic of the sequential algorithm is not conserved and consequently, some properties may be lost. For instance, in this context, the iteration scheme is not preserved, so that ensuring convergence becomes a critical issue (for both the convergence criterion and the number of iterations). Moreover the communication costs may be high for these methods [START_REF] Bertsekas | Parallel and distributed computation, numerical methods[END_REF][START_REF] Baz | An efficient termination method for asynchronous iterative algorithms on message passing architectures -LAAS Report[END_REF]. A third approach, called k-steps in the sequel, is often used in practice. It consists in performing the stopping tests after each group of k iterations (i.e. stopping tests are evaluated at iterations k, 2k, . . . , q.k, . . . ), while the stopping test is not fulfilled. This method performs #I = n iter +k iterations and #T = n iter /k tests. However, since n iter is unknown a priori, the crucial problem is to choose k in order to achieve a good trade-off between #I and #T .

The problem of reducing global synchronization cost without modifying the semantic of the standard sequential algorithm and consequently the behavior of the convergence constitutes the main goal of this paper. This problem is also considered in the field of checkpoint/restart of parallel applications: a global synchronization is sometimes perfomed between all processes before each checkpoint in order to ensure a consistent global state [START_REF] Elnozahy | A survey of rollbackrecovery protocols in message-passing systems[END_REF]. This technique is mainly used to track bugs [START_REF] Kalaiselvi | A survey of checkpointing algorithms for parallel and distributed computers[END_REF]; an error may occur on a processor and causes the crash of the application at timestep n iter . Then, the application is restarted from the last global checkpoint till it crashes again. In order to find the exact time when the bug occured, new checkpoints are computed. A major concern is then to find a strategy that minimizes the number of checkpoint, each one being related to a global synchronization. In the framework of checkpoint restart, simple amortized techniques are used, consisting in doubling the time interval between two consecutive checkpoints [START_REF] Kalaiselvi | A survey of checkpointing algorithms for parallel and distributed computers[END_REF].

In what follows, we concentrate on a general iterative methods, but we restrict to the case where the only global synchronization is the one required by the control of the stopping test. Moreover, we assume that iterations converge to a fixed point, if the stopping test condition is fulfilled at the end of iteration n, it will also be fulfilled at the end of any following iteration.

In this work, we propose a new control technique, original to the best of our knowledge, which delivers the result of iteration n iter after a small number #T = log 1+o (1) n iter of stopping criterion tests while requiring the computation of order #I = n iter + o(n iter ) iterations only. The advantage of our technique lies in the fact that the structure of the initial algorithm is not modified, and consequently the convergence is ensured. On the other hand, the outputs of our algorithm are identical to the outputs of the standard algorithm. Our method is based on an amortized technique inspired from Floyd's and Brent's algorithm to detect periodicity in a sequence [START_REF] Knuth | The Art of Computer Programming[END_REF][START_REF] Cormen | Introduction to Algorithms -Second Edition[END_REF]. It consists in computing two numbers n 1 and n 1 such that n 1 ≤ n iter ≤ n 1 , and then determining the exact value of n iter . Like in the k-steps method, the method we propose performs a limited number of stopping tests only, but tests are no longer performed at regular steps anymore, but rather at steps ρ f (0) , ρ f (1) , • • • , ρ f (i) = n 1 , where 1 < ρ < 2 and f satisfies ∀i, f (i) ≤ i. The choice of f is based on a tradeoff between #I and #T . In this paper, we provide asymptotic results for different choices of f : namely f (i) = i α , 0 < α ≤ 1; f (i) = i log i and f (i) = i log i . We prove that n iter can be determined with very little more than #T = log n iter tests -which is a lower bound for #T -while performing an asymptotic optimal number of iterations #I = n iter + o(n iter ).

In Section 2, we present in detail the k-step method. Sections 3 and 4 are devoted to the presentation of ρ-amortized and generalized ρ-amortized methods which constitute our principal contribution. We close the paper by a conclusion.

Synchronous algorithm and the k-step method

In k-step method, stopping tests are no longer performed after each iteration (like in the standard algorithm), but rather after each group of k successive iterations, for instance at iterations k, 2k, • • • , qk, • • • where q ∈ N * . Let q denote the first index such that the stopping test is fulfilled. Then, this strategy requires qk ≥ n iter iterations and q = n iter k synchronizations. It is clear that the choice of k depends on various criteria: not only on the method itself, but also on the input data. In [START_REF] Bernard (pierre-Eric | Denis) -Large Scale Simulation of Parallel Molecular Dynamics[END_REF], P.E. Bernard et al. use such a strategy to control the convergence and the load balance of a dynamic molecular simulation. In this case, the execution stops as soon as the stopping test is fulfilled, i.e. after qk ≥ n iter iterations. Consequently, the output values of the simulation are different from the results of the standard algorithm, which provides the output values at the end of iteration n iter .

Nevertheless, this method can be modified in order to provide the same outputs as the standard algorithm. In this case, the modified method will be called k-step method. More precisely, the k-step method can be decomposed into two phases.

• During the first phase, the stopping tests are performed after each group of k successive iterations (just as described above).

• During the second phase, the exact iteration n iter is determined by applying the principle of the standard algorithm (stopping tests are then performed after each iteration) but starting from the iteration (q -1)k + 1, for which the context of the iteration has been saved. Since

n iter ∈ [(q -1)k + 1, qk], this phase requires (n iter -(q -1)k) extra iterations and (n iter -(q -1)k) synchronizations.
The two phases lead to a total of

qk + (n iter -(q -1)k) = n iter + k, iterations and q + (n iter -(q -1)k) = q + r, synchronizations
where 1 ≤ r < k, so that at most q + k -1 synchronizations are required.

Therefore, the total cost of the k-step algorithm is given by

T k-step ≤ (n iter + k).t comp + n iter k + n iter mod k .t synch .
It is clear that the choice of k is completely dependent on the value of n iter which is unknown. If k is large, a few number of synchronizations is needed in the first phase but during the second phase, this number becomes large (of order of k). In a dual way, a small value of k leads to a large number of synchronizations during the first phase (of order n iter k ) but during the second phase the number of synchronizations will be reduced. Compared to the standard algorithm, the improvement of the k-step method is given by:

T standard -T k-step ≥ -k.t comp + n iter k .(k -1) + 1 .t synch Therefore, if (the integer) k is chosen so that k 2 k-1 ≤ n iter t synch
tcomp , the k-step method is more efficient (in terms of number of synchronization steps) than the standard one.

ρ-amortized control

This technique is based on an amortized technique inspired from Floyd's and Brent's algorithm to detect periodicity in a sequence [START_REF] Knuth | The Art of Computer Programming[END_REF]. It consists in two phases. The first phase consists in computing two numbers n 1 and n 1 such that n 1 ≤ n iter ≤ n 1 and the second phase consists in determining the exact value of n iter .

Like the k-step technique, this method performs, during the first phase, a limited number of stopping tests. However tests are no longer performed at regular steps, but rather at steps ρ

0 , ρ 1 , • • • , ρ i , • • • . Let ρ k 1 = n 1 denote the first index k 1 such that the stopping test is fulfilled. Then, n 1 ρ = ρ k 1 -1 < n iter ≤ ρ k 1 = n 1 .
Therefore, this phase requires ρ k 1 iterations and k 1 tests. During the second phase, we apply a recursive technique or a dichotomic search for determining the exact value n iter .

Recursive ρ-amortized control

The iteration n iter , is obtained by applying recursively the same process as for the first phase, but starting from the iteration n 1 ρ + 1 = ρ k 1 -1 + 1, whose context has been saved. Stopping tests are performed at iterations

n 1 ρ + 1, n 1 ρ + ρ 1 , • • • , n 1 ρ + ρ k 2 , where k 2 (k 2 ≤ log(n 1 /ρ) ), is the smallest index such that the stopping test is fulfilled. Then, n 1 ρ < n 1 ρ + ρ k 2 -1 < n iter ≤ n 1 ρ + ρ k 2 ≤ n 1 .
Let n 2 = ρ k 2 , then we re-apply the same process, starting from the iteration n 1 ρ + n 2 ρ + 1. Thus, we recursively determine a sequence of indices k l , with l ≤ log ρ n iter , such that

l-1 i=1 ρ k i -1 + ρ k l -1 < n iter ≤ l-1 i=1 ρ k i -1 + ρ k l .
Therefore, by construction of sequence of indices k 1 , . . . , k l ,

n iter = log ρ n iter i=1 ρ k i -1 .
• The total number S of synchronizations is given by S

= k 1 + k 2 + • • • + k l . Since k i ≤ log ρ n iter , S = 1 i=l k i ≤ log ρ (n iter ) 2 .
• The total number of iterations N (in the worst case) is given by

N = l i=1 ρ k i = ρ. l i=1 ρ k i -1 ≤ ρ.n iter .
Finally, the total executing time T ρ-amortized of the ρ-amortized control technique is given by

T ρ-amortized ≤ ρ.n iter .t comp + log ρ (n iter ) 2 .t synch .
Therefore, the recursive ρ-amortized method allows to decrease exponentially the number of global synchronizations, while increasing execution time of order ρn iter .

In what follows, we provide choices for ρ such that the proposed algorithm will be more efficient than the standard one. From previous equations,

T ρ-amorti T standard ≤ ρ t cal t cal + t synch + log 2 ρ (n iter ) n iter t synch t cal + t synch .
Thus, if log 2 ρ (n iter ) < n iter and t synch > t comp , then, by choosing ρ < 1 + t synch t cal , the ρ-amortized method is more efficient than the standard method. Moreover, the performance of the ρ-amortized method is also better than the performance of the k-step algorithm.

Dichotomic ρ-amortized control

The value n iter , that lies in the search interval [ρ k 1 -1 , ρ k 1 ], is obtained using dichotomic technique which consists in log(ρ k 1 -ρ k 1 -1 ) steps and proceeds in the following manner.

First step:

compute all iterations

ρ k 1 -1 + 1, ρ k 1 -1 + 2, • • • , mid, where mid is the middle of the search interval [ρ k 1 -1 + 1, ρ k 1 ],
perform the stopping test only for iteration mid. Indeed, in this case, the iteration n iter lies in one of the two following reduced search intervals:

[ρ k 1 -1 + 1, mid] or [mid + 1, ρ k 1 ]).
For the remaining steps, we apply recursively the same processes, using the new computed search interval until the length of the search interval becomes 1.

Since, at each step we divide the search interval by 2, then

• the number of stopping tests is log(ρ k 1 -ρ k 1 -1 ) < log ρ k 1 = k 1 log ρ,
• the number of iterations is bounded by 2(ρ k 1 -ρ k 1 -1 ).

From phases 1 and 2, we deduce that

• the overall number #I of performed iterations satisfies

#I ≤ ρ k 1 + 2(ρ k 1 -ρ k 1 -1 ) = (3ρ -2).ρ k 1 -1 < (3ρ -2)n iter ,
• the overall number #T of performed synchronizations is bounded by

#T < k 1 + (k 1 ) log ρ = (log ρ + 1)(k 1 )
Since ρ k 1 -1 < n iter ≤ ρ k 1 by construction, then the number S of synchronizations is bounded by: S < (log ρ + 1)(log n iter + 1).

4 Generalized ρ-amortized

Introduction and preliminary results

The scheme we propose in this section is based on the ρ-amortized method. The main difference is that tests are no longer performed at steps

ρ 0 , ρ 1 , • • • , ρ i , • • • but rather at steps ρ f (0) , ρ f (1) , • • • , ρ f (i) , • • • where 1 < ρ < 2 and f satisfies ∀i, f (i) ≤ i.
Our aim is to determine f so that we can find the exact values of n iter

• at minimal synchronization cost (in terms of stopping tests), i.e. as close as possible to θ(log n iter )

• at minimal iteration cost (in term of number of iterations), i.e. as close as possible to n iter .

The proposed scheme is decomposed into two phases:

• during the first phase, the stopping tests are performed after iterations ρ f (0) ,

ρ f (1) , • • • , ρ f (i) , • • • .
We stop the iterations at the first index k such that the stopping test is first fulfilled. Then,

ρ f (k) < n iter < ρ f (k+1) .
At the end of this step, S c = (k + 1) tests have been performed and at most ρ f (k+1) iterations have been computed, so at most ρ f (k+1) -ρ f (k) extra iterations have been computed.

• The second phase is devoted to determine n iter , starting from iteration ρ f (k) + 1. In order to determine the exact value of n iter , we perform a dichotomic search of n iter between ρ f (k) and ρ f (k+1) . At each step of the dichotomic search we need to (re-)compute all the iterates between the lower bound and the middle of the search interval, but we test the stopping criterion for the middle value only. Thus, the overall dichotomic search induces at most

R d = 2(ρ f (k+1) -ρ f (k) )
extra computations and

S d = log 2 ρ f (k+1) -log ρ f (k) < log 2 n iter tests since ρ < 2 and ∀i, f (i) ≤ i.
Finally,

• the overall number R of extra iteration computations is bounded by 3(ρ f (k+1)ρ f (k) ) and therefore R = θ(ρ f (k+1) -ρ f (k) ). In all the sequel, upper bounds are provided on R that are close to n iter .

• The overall number of stopping tests S is S = S c + S d ≤ k + 1 + log 2 n iter tests; in the sequel, in order to bound S we study the number S c of tests performed during the first phase, before the dictotomic search.

In what follows, we provide asymptotic results for R and S c associated to different choices of f , namely

f (k) = k α , 0 < α < 1, f (k) = k log k and f (k) = k log * k .
We prove that it is possible to determine n iter with very little more than log n iter tests, while performing of order o(n iter ) extra iteration computations. Proof: By construction,

ρ k α < n iter ≤ ρ (k+1) α ⇒ k α < log ρ n iter ≤ (k + 1) α ⇒ k < (log ρ n iter ) 1 α ≤ k + 1,
what achieves the proof of Lemma 1.

Lemma 2

The number R of extra computations is of order o(n iter ).

Proof:

ρ f (k+1) -ρ f (k) = ρ (k+1) α -ρ k α = ρ k α ρ (k+1) α -k α -1 = ρ k α log(ρ)αk α-1 + o(k α-1 ) .
By construction, ρ k α < n iter and log(ρ)αk α-1 = o(1) (since α < 1). Therefore,

ρ f (k+1) -ρ f (k) = o(n iter ), what achieves the proof of Lemma 2.
Therefore, if f (k) = k α , the number of stopping criterion tests is poly-logarithmic in n iter , whereas the overall number of computed iterations is of order (n iter + o(n iter )).

f (i) = i

Lemma 4

The number R of extra computations is of order o(n iter ).

Proof:

ρ f (k+1) -ρ f (k) = ρ k+1 log(k+1) -ρ k log k = ρ k log k ρ k+1 log(k+1) -k log k -1 , and since k + 1 log(k + 1) - k log k < 1 log k ρ f (k+1) -ρ f (k) < ρ k log k ρ 1 log k -1
Now, since 1 < ρ < 2 and log k > 1, when n iter → ∞ we have:

ρ 1 log k -1 < log ρ log k + log ρ log k 2 < ρ -1 log k < 1 log k . Moreover, ρ f (k) = ρ k log k < n iter by construction; thus, ρ k+1 log(k+1) -ρ k log k < n iter log k = o(n iter ), which achieves the proof of Lemma 4. Therefore, if f (k) = k
log k , the number of stopping tests is asymptotically close to log n iter log log n iter , whereas the overall number of computed iterations is (n iter + o(n iter )) ∼ n iter .

f

(i) = i log * i
Recall that log * is defined by log * k = min{i ≥ 0, log (i) k ≤ 1}.

Lemma 5

The number S c of tests during the first phase is lower than S c ≤ 1 + log n iter (1 + log * log ρ n iter ).

Proof: By construction,

k log * k < log ρ n iter ≤ k + 1 log * (k + 1) < k log * k (1 + 1 k ). Moreover, k < exp k log * (k) ⇒ log * k < log * k log * (k) + 1 and k > k log * (k) 1 + 1 k ⇒ log * k > log * k log * k 1 + 1 k .
Therefore, log * k -1 < log * (log ρ n iter ) < log * k and, since

k < log * k log ρ n iter < k + 1, k < log ρ n iter 1 + log * (log ρ n iter ) < k + 1.
This achieves the proof Lemma 5.

Lemma 6

The number R of extra computations is of order o(n iter ).

ρ f (k+1) -ρ f (k) = ρ k+1 log * (k+1) -ρ k log * (k) = ρ k log * (k) ρ k+1 log * (k+1) -k log * (k) -1 < ρ k log * (k) ρ 1 log * k -1 .
Now, since 1 < ρ < 2 and log * k > 1, we have:

ρ 1 log * k -1 < log ρ log * k + log ρ log * k 2 < ρ -1 log * k < 1 log * k . Moreover, ρ f (k) = ρ k log * k < n iter ; thus, ρ k+1 log * (k+1) -ρ k log * k < n iter log * k = o(n iter )
, which achieves the proof of Lemma 6.

Therefore, if f (k) = k log * k , the number of stopping tests is asymptotically bounded by log n iter log * log n iter ), whereas the overall number of computed iterations is (n iter + o(n iter )). It is therefore possible to achieve a number of tests close to the lower bound log n iter while performing an asymptotocally optimal number n iter of iteration computations.

Summary of theoretical results

Table 1 summarizes the theoretical complexity results we have obtained in terms of the number of iterations and number of synchronizations, for different choices of the function f where 1 < ρ < 2; n iter denotes the number of iterations required by the standard algorithm which is unknown a priori.

This table shows that n iter can be determined with very little more than log n iter synchronizations, which is a lower bound of the total number of tests, while performing an asymptotic optimal number of iterations (n iter + o(n iter )).

Functions Iterations

Synchronizations Synchronizations S c (first phase) S (both phases)

f (i) = i α , 1 < α < 1 n iter + o(n iter ) ∼ (log ρ n iter ) 1 α ) ∼ log 2 n iter f (i) = i log i n iter + o(n iter ) ∼ log n iter log log n iter ∼ log 2 n iter f (i) = i log * i n iter + o(n iter ) ∼ log n iter log * log n iter ∼ log 2 n iter
Table 1: Asymptotic cost of the three algorithms when n iter increases. The overall number of synchronizations S is S ≤ S c + log 2 n iter ∼ log 2 n iter .

Experiments on a Synthetic Benchmark

When it comes to implementing our methods, a fundamental difficulty, when compared to theoretical proof, lies in the variations of the duration of each process. Even if we may assume a good load-balancing of computations between processors at each step, in practice the time of each iteration may vary due to external effects. Indeed, this phenomenon particularly appears on architectures with a large number of processors. To analyze the impact of a decrease in the number of global synchronizations, we consider for the sake of simplicity that communications occur only for the evaluation of the stopping criterion. For 1 ≤ i ≤ p, let X i be the random variable denoting the duration of the computations performed by processor P i at each iteration step. Let X t i denote the duration of iteration t on processor P i . We assume that all X t i are independent and follow the same distribution law with expectation E(X). If the stopping criterion is evaluated after each step, then the expectation of the duration of each step, from a given step τ , is

E p max i=1 X τ i .
But, if the synchronization is performed only after K steps only, then the expectation of the amortized duration of one step is given by

1 K E p max i=1 τ +K t=τ +1 X t i = E p max i=1 1 K τ +K t=τ +1 X t i ;
which is always smaller than E(max p i=1 X t i ), from Jensen's inequality. Thus, increasing the number K of iterations between two synchronizations will always be better than synchronizing the processes at each step.

Moreover, in previous ρ-amortized algorithms, K increases and becomes very large when n iter increases. Also, using the strong law of large numbers,

lim K→∞ 1 K τ +K t=τ +1 X t i = E(X).
Thus, for large values of K, the duration of each step tends to be equal to the mean value, and the effects of variations due to external effects such as hardware tend to diminish. This property is a major advantage for ρ-amortized methods in cases where the duration of a local iteration on a single processor may vary.

To illustrate this, the following experiment has been performed. A MPI program is performing n iter iterations. At each iteration, each process makes a Lapack dgesv call to solve a N × N system of linear equations. After the resolution, the process i = 1, . . . , p exchanges the result vector with its neighbors i -1, i + 1 (with exception of the processes 0 and p -1, of course, which only send their vectors to the processes 1 and p -2, respectively). After the communication phase, a global synchronization is done through a MPI Allgather, a MPI collective communication that broadcasts to all processes the fusion of pieces of data collected on all other ones. Such an iteration is meant to be representative of a real numerical code, for instance for domain decomposition computations, where a local computation such as a LU factorization is performed, followed by an exchange of the domain intersections between the processes, and eventually by a test of the convergence of the norm of the global solution.

Two implementations of the global synchronization have been made. The first one, called "Global Allgather", broadcasts to all the processes the total solution vector (thus, the volume of the total communication is Θ(pN )). The other implementation first computes locally the norm of each part of the solution vector, and then makes the broadcast of the partial norms (thus with a cost of Θ(p)). We call it "Reduced Allgather". This second solution is more realistic, yet is only valid when the norm of the global vector may be obtained as a linear combination of the partial results.

Based on this MPI template code, it is simple to test and time different schemes for the iterative control of the convergence. We present here the measures made on the standard control algorithm, on the k-step one, and on the generalized ρamortized method, obtained with f (i) = i/ log i (with the value ρ set to 1.95, and k = 10). This includes, when it is necessary, the dichotomic steps. All the necessary operations have been programmed, as in a realistic implementation (copies of the solution for the restarts, norm computation, etc. . . ). The dichotomic search has been programmed such as to perform the worst case number of iterations.

The testing platform is the INRIA's I-Cluster2. It is constituted of 20 nodes Itanium-2, all chosen between the 100 available in order to be connected by a unique switch Fast Ethernet. Notice that this choice, obviously, minimizes the time required for a global synchronization.

We have design a limit experimentation where the kstep method with a small k can be better than the ρ-amortized ones: the number of processor (20) is small and the number n iter = 1000 of iterations is small too. On this platform, the size of the matrix for the Lapack call has been fixed to N = 300, and the convergence value n iter = 1000 is imposed. Typically, such a program runs in about 3 minutes on 20 I-Cluster2 nodes. Each version of the control algorithm has been run 10 times, in order to obtain statistically significant results.

For each run, the total number of iterations R has been measured, as well as the number of synchronization and the time that they lasted. Of course, the timing results highly depend on the hardware and software implementation. Table 2 presents the results obtained on each algorithm, for each one of the two implemented Allgathers.

Of course, the number of iterations and synchronizations is equal with both Allgathers. The standard deviations are not presented in the table, because they are insignificant. These measures show, without surprise, the number of synchronizations that are spared by the k-step and i/ log i algorithms. Due to the very low value for k, the k-step algorithm appears as specially interesting. Regarding timings, as far as 38 seconds are gained, on a total of some 180 (i.e. 20%), due to a reduced use of synchronizations, in the case of the Global Allgather. With an optimized Reduced Allgather, the gain is less (some 2 seconds), yet visible. These measures, even for a template application running on few nodes, may be seen as promising. The mere increase in the number of nodes or, even simpler, the same runs, but on 20 nodes that would not be interconnected by a complete network, would simply amplify the gain that is already measured with this limited benchmark.

Conclusion

Minimizing the number of global synchronizations in a distributed computation enables to increase efficiency of synchronous algorithms. In this paper, we focus on iterative methods where a global stopping criterion is evaluated at each step. We analyzed various strategies in order to decrease the number of synchronizations while ensuring that the computed output result is the same than the one delivered by the sequential method. We propose a new control technique named ρ-amortized, original to the best of our knowledge, where the global test is performed only after iteration ρ f (i) with 1 < ρ < 2 (f being a non-decreasing function). Analyzing various trade-offs for ρ and f , we prove that, when the total number of iterations n iter performed by the sequential method is unknown, the whole number of global synchronizations may be reduced close to the lower theoretical bound (log n iter ) on the overall number of tests, while performing an asymptotically optimal number of iterations (n iter + o(n iter )).

Experiments have been run on a synthetic benchmark, meant as a template for scientific computing based on iterative methods. Running on a small number of nodes, using a complete interconnecting network, they have shown a clear gain in time spent doing synchronization. This implementation, as well as theoretical considerations, show that the proposed methods are robust and may be implemented

4. 2 f 1 Lemma 1

 211 (i) = i α , where 0 < α < The number S c of tests during the first phase is lower thanS c ≤ 1 + (log ρ n iter ) 1 α tests.

log i Lemma 3

 3 The number S c of tests during the first phase is lower than S c ≤ log n iter . log(log n iter )) 1 + O 1 log log log n iter . Proof: By construction, ρ k log k < n iter ≤ ρ k+1 log(k+1

  ) . log ρ n iter log k ⇒ k < log ρ n iter . log(log ρ n iter . log k) ⇒ k < log ρ n iter . log log ρ n iter + log ρ n iter . log log k Since k = O(log n iter ), we deduce that k < log ρ n iter log log ρ n iter . 1 + O 1 log log log n iter ∼ n iter →∞ log ρ n iter log log ρ n iter which achieves the proof of Lemma 3.

	Therefore,
	k log k	< log ρ n iter ≤	k + 1 log(k + 1)	⇒ k <

Table 2 :

 2 Number of iterations, of synchronizations, and associated cumulated time for each, using the Global and the Reduced Allgather.

	Global Allgather
		Standard Algorithm k-step algorithm i/ log i algorithm
	n iter	1000	1000	100
	R	1000	1047	1078
	S	1000	28	72
	S c	1000	?	?
	S d	0	?	?
	T iter (mean)	185.7 s	181.4 s	188.2 s
	T synch (mean)	40.2 s	1.0 s	2.1 s
	T iter /n iter (mean)	185.7 ms	181.4 ms	188.2 ms
	T synch /n iter (mean)	40.2 ms	1.0 ms	2.1 ms
	Reduced Allgather
		Standard Algorithm k-step algorithm i/ log i algorithm
	T synch (mean)	2.8 s	0.9 s	2.8 s

efficiently.

This work is supported by the "Comité Mixte Franco-Marocain -Action Intégrée" MA/01/19.