Patrick Loiseau
email: patrick.loiseau@eurecom.fr

Claire Médigue
email: claire.medigue@inria.fr

Paulo Gonçalves
email: paulo.goncalves@ens-lyon.fr

Najmeddine Attia
email: najmeddine.attia@inria.fr

Stéphane Seuret
email: seuret@u-pec.fr

François Cottin
email: francois.cottin@inserm.fr

Denis Chemla
email: denis.chemla@bct.aphp.fr

Michel Sorine
email: michel.sorine@inria.fr

Julien Barral
email: barral@math.univ-paris13.fr

Large deviations estimates for the multiscale analysis of heart rate variability

Keywords: Heart rate variability, multiscale analysis, multifractal analysis, large deviations

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

The characteristics of the inter-beat-times (RR) signal, and in particular of heart rate variability (HRV), constitute an important fingerprint of the heart state. They reflect complex mechanisms, controlled at short term by the autonomic nervous system or at longer term by circulatory hormones, for instance. Since the pioneering works [START_REF] Womack | The analysis of respiratory sinus arrhythmia using spectral analysis and digital filtering[END_REF][START_REF] Sayers | Analysis of heart rate variability[END_REF], they have been extensively studied in the last decades with various mathematical methods (Fourier analysis, time-frequency analysis, complex demodulation, etc.). At the light of this long series of works, it is now unanimously accepted that spectral analysis of the RR signals is insightful to characterize HRV at the different scales of interest for the cardiologists, and to detect different heart pathologies. Multiscale approaches are intrinsically different from classical spectral analyses as these latter strongly rely on data stationarity and focus on a very limited (two, sometimes three) range of time scales of interest through the so-called (very-)low and high frequency bandwidths [START_REF] Akselrod | Power spectrum analysis of heart rate fluctuation: a quantitative probe of beat-to-beat cardiovascular control[END_REF][START_REF] Malliani | Cardiovascular neural regulation explored in the frequency domain[END_REF][START_REF] Van Laar | Spectral analysis of fetal heart rate variability for fetal surveillance: review of the literature[END_REF][START_REF] Mainardi | On the quantification of heart rate variability spectral parameters using time frequency and time-varying methods[END_REF]. Alternately, different techniques accounting for multiple time-scales have recently been applied to heart rate variability: a wavelet method is used in [START_REF] Ivanov | Scaling behaviour of heartbeat intervals obtained by wavelet-based time-series analysis[END_REF] to characterize the nonstationary behavior of the RR signals, a multi-scale probabilistic analysis of RR fluctuations is proposed in [START_REF] Kiyono | Multiscale probability density function analysis: Non-gaussian and scale-invariant fluctuations of healthy human heart rate[END_REF], while it is the intrinsic complexity of heart rate data that is addressed using entropy measures in [START_REF] Ferrario | Comparison of entropy-based regularity estimators: Application to the fetal heart rate signal for the identification of fetal distress[END_REF], linear or non-linear models in [START_REF] Baselli | Cardiovascular variability signals: towards the identification of a closed-loop model of the neural control mechanisms[END_REF][START_REF] Jarisch | Statistical modeling of fetal heart rate variability[END_REF][START_REF] Signorini | Linear and nonlinear parameters for the analysis of fetal heart rate signal from cardiotocographic recordings[END_REF] or stochastic point processes in [START_REF] Barbieri | A pointprocess model of human heartbeat intervals: new definitions of heart rate and heart rate variability[END_REF]. Also, a collection of studies elaborate on the fractal dimension (or equivalently on the local regularity) of RR signals to identify a variety of physiological conditions [START_REF] Hopkins | A comparative study of fetal rate variability analysis techniques[END_REF][START_REF] Nakamura | Autonomic control of heart rate during physical exercise and fractal dimension of heart rate variability[END_REF][START_REF] Yamamoto | On the fractal nature of heart rate variability in humans: effects of data length and β-adrenergic blockade[END_REF][START_REF] Nunes Amaral | Behavioral-independent features of complex heartbeat dynamics[END_REF].

Multifractal analysis formalized the intuitive idea that the local regularity could itself vary erratically along time and significantly deviate from its baseline (see [START_REF] Lopes | Fractal and multifractal analysis: A review[END_REF] for a medically oriented introduction to this theory and a vast list of biomedical applications of it). For instance, interesting works demonstrated the relevance of multifractal-derived indicators to diagnose possible cardiovascular stresses or to assess the autonomic nervous system development of foetuses [START_REF] Nakamura | Local holder exponent analysis of heart rate variability in preterm infants[END_REF][START_REF] Nakamura | Multifractal analysis of heart rate variability[END_REF].

Nonetheless and despite all this anthology of promising results, still today, the complexity of heart beat rate remains elusive. The present contribution lies in the vein of multifractal analysis and leverages recent advances in large deviations theory to improve the precision of state-of-the-art estimators of large deviations spectra from heart beat rate time series. In particular, multifractal analysis is used here to measure the scaling or non-scaling properties of RR signals, and leads to remarkable discrimination between several behaviors, namely congestive heart failure, atrial fibrillation and physiological behavior.

Multifractal analysis is a multi-scale analysis which represents the structure of the variability of irregular (fractal ) signals such as the RR signal, through scaling functions and spectra which describe the statistics of the signal's roughness. In some signals, called monofractal, the roughness is described by a single exponent; whereas in the general case of multifractals, there is more heterogeneity and the signal exhibits a full interval of roughness exponents. The question of mono versus multifractality of RR signals in different types of subjects have been largely debated, see e.g. [START_REF] Kiyono | Critical scale invariance in a healthy human heart rate[END_REF][START_REF] Ching | Multifractality and scale invariance in human heartbeat dynamics[END_REF][START_REF] Makowiec | Multifractal estimates of monofractality in RR-heart series in power spectrum ranges[END_REF]. The notion of monofractality, however, is mostly theoretical and a range of roughness exponents is always observed in practice. Therefore, we do not discuss this question here and we focus on the possible shapes of multifractal spectra and on the reliability of their estimates.

Two main lines of thought coexist in multifractal analysis. The most popular in HRV analysis is based on the Legendre transform of the scaling function, termed Legendre spectrum. This approach always assumes scaling properties of the data. Used for instance in [START_REF] Ivanov | Multifractality in human heartbeat dynamics[END_REF] and [START_REF] Sassi | Multifractality and heart rate variability[END_REF], it revealed distinct multifractal behaviors for two sets of individuals: more linear scaling functions, i.e. narrower Legendre spectra, for subjects that suffer heart failure (or atrial fibrillation) than for healthy subjects. The other axis of multifractal analysis relies on the so-called large deviations spectrum, which in a nutshell, amounts to compute at each observation scale the logarithmic frequency of occurrence of the different values of the roughness exponent, and was first considered in [START_REF] Mandelbrot | Multiplications aléatoires itérées et distributions invariantes par moyenne pondérée aléatoire[END_REF], in the context of Mandelbrot cascades. Such an analysis is performed in [START_REF] Meyer | Self-affine fractal variability of human heartbeat interval dynamics in health and disease[END_REF] where the authors also analyze two groups of subjects (healthy and heart failure), using an original definition of the roughness exponent based on detrended fluctuation analysis (DFA, a method introduced and studied in [START_REF] Peng | Mosaic organization of dna nucleotides[END_REF][START_REF] Hu | Effect of trends on detrended fluctuation analysis[END_REF][START_REF] Chen | Effect of nonstationarities on detrended fluctuation analysis[END_REF]). They observe that heart failure is characterized by a non-concave large deviations spectrum whereas the spectrum of healthy subjects is strictly concave. Although it is an intricate difference to interpret, this study remarkably illustrates that, compared to the Legendre spectrum -which by definition is a convex-shaped function -the large deviations spectrum conveys additional information.

In practice though, to face the tricky problem of estimating large deviations spectra the authors of [START_REF] Meyer | Self-affine fractal variability of human heartbeat interval dynamics in health and disease[END_REF] used an ad-hoc method that implies an arbitrary choice of crucial and extremely sensitive parameters, and moreover assumes the implicit existence of a scale-invariance property over a reasonable scale range.

To overcome these limitations, [START_REF] Makowiec | Multifractal estimates of monofractality in RR-heart series in power spectrum ranges[END_REF] proposed to estimate Legendre spectra at different scales by Legendre transforming a sequence of scaling functions derived from a piecewise regression over small scale intervals. However, this procedure yields unstable estimates that moreover, do not easily relate to well-known characteristics of HRV.

In this work, we exploit the theory of large deviations to derive a systematic and theoretically sound procedure to estimate the large deviations spectra at each scale. Our estimation procedure presents a number of advantages as compared to previous methods: it is fully adaptive, it relies on strong theoretical justification and it provides one spectrum estimation for each scale. This last point is maybe the most precious, as it permits to question the presence of a scaling behavior that is implicitly assumed in all prior studies. Applied to RR time series, it will quantify, in the range of scales corresponding to short-term variability, different non-scaling behaviors for healthy subjects and subjects with atrial fibrillation, while subjects with heart failure turn out to exhibit scaling properties close to that of a synthetic signals obtained from standard multifractal statistically self-similar objects. This observation has important consequences, notably in the context of adequate models identification.

We also show that the proposed algorithm for large deviations spectrum estimation is sensitive to events in the signal, via the presence of nonconcavities, imperceptible to the Legendre spectrum. We illustrate this property in the context of RR signals via the presence/absence of extrasystoles.

Let us stress here that a major advantage of the proposed data treatment is that it does not rely on a stationarity assumption and therefore on the sensitive choice of a short stationary time window. Typically, we will consider here hours-long signals, which gives the proposed method a valuable robustness as compared, e.g., to methods based on Fourier analysis.

The paper is organized as follows. We describe in Section 2 the theoretical aspects of multifractal analysis and the algorithm that we use to estimate large deviations spectra. In Section 3, we demonstrate the sensitivity of the algorithm to interesting properties of the signal, first on the controlled example of the Brownian motion in multifractal time (Sec. 3.1) and then on RR signals presenting extrasystoles (Sec. 3.2). In paragrah 3 of the same Section, we study the scaling or non-scaling behavior of RR-signals. We conclude in Section 4.

Multifractal analysis of signals

We consider dyadic scales: for an integer n ≥ 1, the dyadic interval

I n (k), 1 ≤ k ≤ 2 n is defined as I n (k) = [(k -1)2 -n , k2 -n ] and its length is |I n (k)| = 2 -n .
The set of all dyadic intervals at scale n is denoted by G n .

Theory 2.1.1. Various notions of roughness grain exponents

A roughness exponent reflects the amplitude of variations of signal X within a dyadic interval. Several definitions can be used.

Oscillation-based grain exponent. The oscillation of X in the interval

I n (k) is defined as Osc X (I n (k)) = sup{|X(u) -X(v)| : u, v ∈ I n (k)},
and the corresponding oscillation grain exponent is defined as

α X,Osc (I n (k)) = log Osc X (I n (k)) log 2 -n . (1)
It is the unique exponent α = α(X, n, k) such that the oscillation behaves as a power law of the interval size: Osc X (I n (k)) = 2 -nα . This exponent is very natural to measure the signal's roughness at points where it is not smooth; at points where the signal is smooth, one can identify a smoothness exponent by using oscillations of higher orders.

Wavelet leader-based grain exponent. The (discrete) wavelet coefficients of the signal X are defined as w X (n, k) = ψ n,k , X , where ψ n,k (t) = 2 n ψ 0 (2 n tk), k ∈ Z, n ∈ Z, and •, • stand for the inner product in the L 2 ([0, 1]) signal space. The so-called wavelet ψ 0 ∈ L 2 ([0, 1]) is compactly supported (it can be defined using a finite dyadic filter [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF], and is characterized by its number of vanishing moments N ψ > 1, an integer such that t k ψ 0 (t) dt = 0, ∀k = 0, ..., N ψ -1. Provided some now well-known assumptions on ψ 0 , the wavelet coefficients w X (n, k) of a signal X entirely characterize X, which can be reconstructed as X = n∈Z,k∈Z 2 -n/2 w X (n, k)ψ n,k . Moreover, the decreasing rate of these coefficients, when n goes to infinity, characterizes both the global and the local regularity of the signal X, it is thus also natural to use the w X (n, k) as a measure of the roughness of a signal. Nevertheless, theory and practice show that wavelet coefficients are not easy to deal with (for instance, a coefficient can be very small, generating numerical artefacts). This leads to rather study the so-called wavelet-leaders [START_REF] Jaffard | Wavelet techniques in multifractal analysis, Fractal Geometry and Applications: A Jubilee of Benoît Mandelbrot[END_REF][START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF] (or to use the closely related wavelet transform modulus maxima method [START_REF] Muzy | Multifractal formalism for fractal signals: the structure fonction approach versus the wavelet transform modulus maxima method[END_REF]), which correspond to the thinned series:

L X (n, k) = sup{w X (n , k ) : I n (k ) ⊂ I n (k -1) ∪ I n (k) ∪ I n (k + 1)}.
As for the oscillation, the wavelet-leader grain exponent is then defined as:

α X,L (I n (k)) = log L X (n, k) log 2 -n . (2)
It is known that this exponent is asymptotically equivalent to the oscillationbased exponent when n goes to infinity for values smaller that 1, i.e. at "rough points".

Multifractal detrended fluctuation analysis. The two aforementioned grain exponents permitted to state precise theoretical results on typical multifractal objects. Yet, a more empirical approach, known as multifractal detrended fluctuation analysis (MDFA) turned out to be more popular in the context of RR signals analysis [START_REF] Meyer | Self-affine fractal variability of human heartbeat interval dynamics in health and disease[END_REF][START_REF] Makowiec | Multifractal estimates of monofractality in RR-heart series in power spectrum ranges[END_REF][START_REF] Stanley | Statistical physics and physiology: monofractal and multifractal approaches[END_REF][START_REF] Kantelhardt | Multifractal detrended fluctuation analysis of nonstationary time series[END_REF]. In this method, whose principle is close to the wavelet method, one considers a discrete signal X and the coefficients

D X (I n (k)) = Ñ 1 #I n (k) x∈In(k) X(x) -P (dp) In(k) p é 1/p
, where p ≥ 2 is an integer, P

In(k) is the polynomial trend of X of degree d p ≥ 1 in the interval I n (k), and here #I n (k) stands for the cardinality of the set of points of I n (k) at which X is defined. This yields another grain exponent:

α X,D (I n (k)) = log D X (n, k) log 2 -n . (3)
According to previous empirical studies, p = 2 and d p = 2 seem to be the most appropriate values to use in practice, thus we opt for this same choice.

As we shall see, the three definitions yield very similar results in practice; in particular the observed roughness exponents on RR signals are smaller than 1. Therefore, we will favor the most natural choice of the oscillationbased exponent.

Multifractal spectra

Multifractal analysis aims at characterizing the heterogeneity of a roughness exponent of a signal via a multifractal spectrum. This distribution-like approach strongly relates to large deviations theory, which allows quantifying the probability of rare events. Classically, there are three different definitions of multifractal spectra.

The Legendre spectrum. Consider the sequence of scaling functions defined on R associated to the grain exponent α X :

τ X,n (q) = - 1 n log 2 I∈Gn 2 -nqα X (I) . (4)
The Legendre spectrum is defined as the Legendre-Fenchel transform of the function τ X (q) = lim inf n→∞ τ X,n (q):

τ * X (α) = inf q∈R {αq -τ X (q)}. (5)
By construction, it is always a concave function.

The large deviations spectrum. Another possibility to describe the heterogeneity of roughness exponents is to count at each scale n, the number of dyadic intervals where the grain exponent is close (up to ) to a particular value, and to define the "logarithmic frequency":

f X,n (α, ) = log #{I ∈ G n : α -≤ α X (I) ≤ α + } log 2 n . ( 6
)
The following double limit function defines the large deviations spectrum:

f X (α) = lim →0 lim sup n→∞ f X,n (α, ), (7)
a sort of continuous histogram (in term of the logarithmic frequency) measuring the occurrence of each exponent α.

The Hausdorff spectrum. The term multifractal was given in reference to the third definition of the spectrum which corresponds to the Hausdorff dimension of the iso-α sets:

g X (α) = dim H {t ∈ [0, 1] : lim inf n→∞ α X (I n ( t2 n )) = α}. (8)
Based on geometric information, this theoretical spectrum is impossible to estimate in practice, and we will naturally discard it in the sequel.

Relationship between the different spectra. One always has g X (α) ≤ f X (α) ≤ τ * X (α) [START_REF] Riedi | Multifractal processes[END_REF][START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF]. Moreover, in a number of specific cases where there is some self-similarity, it is known (for the oscillation-based exponent and the wavelet-based exponent) that the three spectra coincide [START_REF] Rand | The singularity spectrum f (α) for cookie-cutters[END_REF][START_REF] Jaffard | Multifractal formalism for functions. i. results valid for all functions. ii self-similar functions[END_REF][START_REF] Barral | Multifractal analysis of complex random cascades[END_REF] and the multifractal formalism is then said to hold. Moreover, if the sequence of scaling functions converges on some interval Q to τ X (q), and the limit is differentiable, then for any α = τ X (q), q ∈ Q, we have:

f X (α) = lim →0 lim n→∞ f X,n (α, ) = τ * X (α). (9)
In words, the large deviations spectrum is obtained by using limit only (no lim sup is required) and equals the Legendre spectrum [START_REF] Dembo | Large Deviations Techniques and Applications[END_REF].

Let us stress that the collection of logarithmic frequencies f X,n (α, ) yields a refined picture of the distribution of grain exponents at each scale, sometimes disclosing information that remains hidden in the function τ * X,n . Typically, if the multifractal formalism does not hold, the relation f X ≤ τ * X reduces to a strict inequality, and it may occur that the estimation of f X exhibits a non concave spectrum inevitably different from that of τ * X , as it is the case in the numerical studies of [START_REF] Riedi | Multifractal Properties of TCP traffic: a numerical study[END_REF][START_REF] Meyer | Self-affine fractal variability of human heartbeat interval dynamics in health and disease[END_REF].

Our approach will only exploit the functions τ X,n (q), τ X,n (q) and τ X,n (q) and provide large deviations spectra estimates as the parametrized graphs (α n (q), f X,n (α n (q), n (α n (q))), where α n (q) = τ X,n (q) and n (α n (q))) is expressed suitably through τ X,n (q) as in [START_REF] Baselli | Cardiovascular variability signals: towards the identification of a closed-loop model of the neural control mechanisms[END_REF].

Practical estimation of multifractal spectra

When analyzing a discrete signal (corresponding to sampled data or to a genuine time series, like that of the RR intervals), the number of available scales is limited above by the finest resolution and limited below by the signal size. It is then only possible to compute the different quantities at finite scales n. To estimate the Legendre spectrum, one simply computes the quantities:

τ * X,n (α) = inf q∈R {αq -τ X,n (q)},
at the different scales of interest. On the other hand, estimation of the large deviations spectrum is more complex because of the double limit involved in definition [START_REF] Ivanov | Scaling behaviour of heartbeat intervals obtained by wavelet-based time-series analysis[END_REF]. A method has been proposed in [START_REF] Canus | Continuous large deviation multifractal spectrum: definition and estimation[END_REF][START_REF] Véhel | On various multifractal spectra[END_REF] (and used in [START_REF] Meyer | Self-affine fractal variability of human heartbeat interval dynamics in health and disease[END_REF]), that involves ad-hoc tuning of arbitrary parameters. Moreover, it only produces one smoothed spectrum corresponding to a regression procedure over a particular scale range. In the present work, we use a fully adaptive algorithm proposed in [START_REF] Barral | On the estimation of the large deviations spectrum[END_REF], to estimate the large deviations spectrum at each scale n. The idea ruling this algorithm, grounded on large deviations theory and on geometric properties of multifractal objects, is to parameterize in expression [START_REF] Mainardi | On the quantification of heart rate variability spectral parameters using time frequency and time-varying methods[END_REF], the bin width to the analyzing scale n and to the grain exponent value α n = τ X,n (q) as follows:

n (α n ) = Ã log log(n) τ X,n (q) n log 2 , (10)
the practical computation of τ X,n (q) (which exhausts the range of possible α at each scale) and τ X,n (q) being directly based on formulas obtained from the successive differentiation of (4). Doing this for values of q taken over a thin enough grid provides the parametrized spectrum (α n (q), f X,n (α n (q), n (α n (q))) announced at the end of the previous section. Then, it is shown in [START_REF] Barral | On the estimation of the large deviations spectrum[END_REF], that for well-known multifractal signals, this method does converge towards the actual large deviations spectrum f X (α), which makes this adaptive algorithm a theoretically sound method to estimate large deviations spectra, and therefore a good candidate to investigate multifractal properties of RR signals. It is important to notice here that at a finite scale, the inequality f X,n (α n , n (α n )) ≤ τ * X,n (α n ) may not be realized for all α's (only f X (α) ≤ τ * X (α) is true in general). Normalization of the signal. In contrast to most known methods, this algorithm uses roughness exponents computed at a scale n, and not from a regression over a range of values. This is an interesting property that produces a sequence of spectra at the different scales. However, it necessitates to properly and systematically normalize the signal so that the range of estimated exponents defined in (1) does not drift with scale. Here, we chose a multiplicative factor such that after normalization, the signal gets as close as possible to a self-affine fractal. More precisely, we start estimating the fractal dimension D g of its graph by linearly regressing the quantity n + log 2 I∈Gn Osc X (I) = n(1-τ X,n (1)) against the scale n, so that n(1-τ X,n (1)) behaves like A + nD g . Then, for an arbitrary fixed scale n 0 belonging to this regression range, we divide the signal X by the unique constant C yielding n 0 (1 -τ X/C,n 0 (1)) = n 0 D g (i.e. zeroing the intercept A). Numerically, this normalization procedure is almost insensitive to the chosen scale n 0 , and it has the advantage to be data-driven, necessitating no a priori knowledge on the signal. Naturally though, the closer the quantity I∈Gn Osc X (I) is to a power law behavior (i.e., n(1 -τ X,n (1)) ∼ A + nD g ), the less sensitive is the normalization with respect to n 0 .

Other normalizations are possible that would only impact the resulting spectra up to a slight horizontal translation and would not affect the results of the next section on (non-)scaling and non-concavity.

The advantage of the proposed method is that it permits, in practice, to estimate one large deviations spectrum for each scale n, and thus to question in an interesting manner the scale-invariance property through the superimposition of the different spectra. Moreover, as we shall see, it allows nonconcave spectrum estimates that are reminiscent of certain signal specificities (which remain undisclosed in a Legendre spectrum).

Results and discussion

A toy example: the Brownian motion in multifractal time

Before discussing multifractal analysis of RR signals and its interpretation, we illustrate the different aspects of the estimation's method on a controlled multifractal signal: the Brownian motion in multifractal time. The ordinary Brownian motion B(t) (Fig. 1-(a)) is (up to a multiplicative factor) the only centered Gaussian process with independent and stationary increments. Theoretically it is a monofractal signal, with the single possible roughness exponent H = 0.5. To synthesize a controlled multifractal process, we form the compound process X(t) = B(M (t)), where M is a dyadic multiplicative cascade with equally probable weights W ∈ {0.2, 0.8}. A realization of the resulting process, called Brownian motion in multifractal time, is displayed in Fig. 1-(c). The theoretical multifractal spectrum of X reads f X (α) = f M (2α) almost surely, where f M is the multifractal spectrum of the cascade M . As the multifractal formalism holds here, f M readily derives from the Legendre transform of τ M (q) = -log 2 (0.2 q + 0.8 q ) (see [START_REF] Riedi | Multifractal processes[END_REF] for more details). Numerically, we generate the signals at a sufficiently fine resolution J = 18, such that the generated trace over the 2 18 + 1 points reasonably approaches the continuous limit process. This imposes the maximum analyzing scale n = J . Fig. 2 illustrates the multifractal analysis (corresponding to the oscillationbased roughness exponent): at scale n, the oscillation is computed in each interval of size 2 -n (there are 2 n consecutive closed dyadic intervals of length 2 J-n + 1 points) from which the corresponding roughness exponent α is deduced. Fig. 2 shows that although the oscillations' amplitude naturally decreases with n, the range of the roughness exponents stabilizes across the scales, a characteristic expression of scale-invariance property.

Using the algorithm described in the previous section, we can then estimate for each scale n the large deviations spectrum of (7) that describes the statistical repartition of the roughness exponents. resulting spectra for scales n = 13 to n = 8. The fair superimposition of the estimated spectra at all scales with the theoretical curve, is the evidence of the underlying scaling properties of B(M (t)), at the same time that it validates our choice for n (α n ) in Eq. ( 10) and the proposed normalization procedure. Fig. 3-(b) and 3-(c) show that the estimation algorithm performs equally well with other roughness exponents: namely with the wavelet-leader-based and with the MDFA exponents. Therefore, and without loss of generality, we choose for the rest of our study to consider oscillation-based exponents only, as they are the most natural to interpret. Let us finally mention that here, and in all the following experiments, we use values of q between -100 and +100 to compute the spectra. Whereas most authors restrict to values of q between -5 and +5, this allows us to fully explore the spectrum over its entire support. We now illustrate the sensitivity of the estimation algorithm to various signal transformations. We present in Fig. 4 the effect of signal quantization on the spectrum estimate. As the quantization step-size increases, the left part of the spectrum remains sensibly unchanged while its right part develops a stepping effect due to the rarefaction of possible values for the roughness exponent. As the minimum non-zero oscillation that can be observed at each scale n is constant and equal to the quantization step, the maximum estimated roughness exponent increases at the coarser scales. Fig. 5 illustrates the effect of brutal variations on the spectrum estimates. At first (see Fig. 5-(a)), we add to the signal three local events similar in shape to extrasystoles in cardiac signals. An extrasystole is a premature cardiac excitation, followed by a compensatory pause, more or less marked, and results in two successive beats, a short RR followed by a longer one. Fig. 5-(b) provides a first illustration one of the main asset of our estimation algorithm: it faithfully reproduces non-concavely shaped spectra. In the present case, the non-concavity is due to a few very large oscillations creating locally very small exponents (around 0.1 here) which significantly depart from the rest of the exponents. Then, by spectrum continuity, these exponents create a non-concave variation that cannot be observed with the Legendre spectrum. The rest of the spectrum remains unaltered. To illustrate further such a phenomenon that occurs in RR signals of subjects with heart failure, we add more events of extrasystole kind, with two different amplitude in Fig. 5-(c). Then, the corresponding spectrum in Fig. 5-(d) presents two clear concavity changes, the first one caused by the biggest "extrasystoles", the second one by the smallest ones.

RR physionet signals study

We now turn to the analysis of RR signals and show how the proposed estimation algorithm reveals so far unobserved properties of these signals.

In this section, most analyzed signal are of length 2 13 + 1 points and we focus on the scales 7 to 12 that are related to the short term variability of the signal. This window size (J = 13) has been chosen for it allows to observe weak non-stationarity of our signals: the neighboring windows show almost the same characteristics. Tab. 1 describes the essential characteristics f (α) of the considered scales and their (assumed) relationship with the different physiological factors. Contrarily to the Brownian motion in multifractal time studied in the previous section, RR signals contain strong periodic components ("physiological oscillations") at the scales considered in our analysis (see Tab. 1). Therefore, it is not clear that detrending is not needed here before performing an oscillation-based multifractal analysis. To check that such detrending is indeed not necessary, we performed the same analysis with the waveletleader-based exponent (which is blind to trends) instead of the oscillationbased exponent, and verified that all the results presented in this section remain similar (for both non-concavity and non-scaling). In the rest of this section, we thus present only the results obtained with the oscillation-based exponent.

Database

We used signals from the Physionet public database2 . They consist of a set of 15 heart beat (RR-interval) times series, 5 in health, 5 in congestive heart failure, and 5 in atrial fibrillation, recorded at 125Hz. Each time series is about 24 hours long.

Concavity changes due to extra-systoles

Non-concavities in RR signals have been discussed in [START_REF] Meyer | Self-affine fractal variability of human heartbeat interval dynamics in health and disease[END_REF], where they are claimed to be a sign of heart failure. The authors suggest that the concavity change might be due to the lumping of two sub-spectra, respectively related to the parasympathetic and sympathetic activities. Here, we emphasize that such concavity changes can be mainly caused only by the presence of extrasystoles.

We consider in Fig. 6 the case of a subject with congestive heart failure. The corresponding RR time series exhibits several extrasystoles arising from two different locations in the heart and with approximately two different coupling intervals (i.e. whose amplitudes take on approximately two different values). These two types of events produce in the estimated large deviations spectrum (red curve) two inflection points whose location on the α and f (α) axes relates to the extrasystole's amplitude and to the frequency of their occurrence, respectively (this behavior motivated the operations done in Fig. 5). Removing the extrasystoles in the original signal smooths out these non concavities.

Scaling behavior study

In this section, we discuss the fundamental assumption that the analyzed RR signals possess a scale invariance structure. It is this assumption that bears all estimation methods that rely on a systematic linear regression of the logarithmic structure function defined in (4). In the opposite, as our estimator yields a large deviations spectrum at each analyzing scale, it permits to investigate this scale invariance property. In particular, the scale-dependent spectra should superimpose in presence of scale invariance. Fig. 7 presents the analysis of a RR signal recorded on a control subject. The non superimposition of the large deviations spectra estimated for different scales casts doubt on the scale invariance hypothesis, even though the logarithmic structure functions n • τ X,n (q) do seem to behave linearly for q's close to zero. Indeed, the procedure used to detect the linear behavior of n • τ X,n (q) is applied to each individual q in some range (e.g., q ∈ [-5, 5]), which does not make sure that this scaling is uniform in q, i.e., that the functions τ X,n possess scale invariance. In this sense, considering the parameters q individually only provides scale invariance in a weak sense. In contrast, our approach is, by design, sensitive to uniform scale invariance (all well-known multifractal models do have scale invariance in such a strong sense).

The non-scaling observed on Fig. 7 naturally leads to study how this phenomenon is altered, or not, in presence of pathology, and to have some measurement of non-scaling/scaling. To fix ideas, Fig. 8 presents spectra estimated in healthy, atrial fibrillation and congestive heart failure, after removing extrasystoles (this removing is natural since the distribution of extrasystoles is so erratic that it unavoidably creates non-scaling in the left part of the spectrum). The global non-scaling behavior seems strengthened with atrial fibrillation, while with congestive heart failure, some scale invariance appears (at least in the left part corresponding to the biggest oscillations and thus less sensitive to quantization).

In order to better quantify these observations, we propose to calculate the area A(n) beneath each large deviations spectrum at each scale n, and to plot the evolution with the analyzing scale of the ratio A(n)/A(n max ) for n in a fixed interval [n min , n max ].

Fig. 9-(a) compares the evolution of this simple empirical criterion for the RR signals in the 3 different classes of subjects, and it also includes the same calculations for 5 simulations of Brownian motion in multifractal time, which is a pure statistically self-similar process. As expected, the calculated area remains almost constant with n for BMMT (green curve), and it is close to follow the same behavior in congestive heart failure (blue curve), while it significantly increases in healthy (red curve), a phenomenon considerably strengthened in atrial fibrillation (black curve). This discrepancy between scaling versus non-scaling behavior turns out to be better emphasized by restricting the calculation of this area criterion to the left part of the spectra, Fig. 9-(b). Moreover, in the left part of the spectra, our measurement of the (non-)scaling property is almost insensitive to the acquisition sampling frequency. For the right part of the spectra, Fig. 9-(c) shows that it is quite sensitive to the quantization of the signal.

It is also worth noticing that the proposed criterion is insensible to the chosen normalization procedure applied prior to estimating the large devia- scaling property by estimating via linear regression the slope of the graph A q (n)/A q (n max ) versus n. The closer the slope is to 0, the more we consider that the signal possesses scaling properties for the exponents in [-q, q], [0, q] and [-q, 0]. Fig. 10 illustrates this new scaling property measure. As for the simpler analysis of Fig. 9, computing surfaces only on the left part of the spectra (i.e., using the mapping A 0,q as in Fig. 10-(b)) yields the best distinction between different classes due to a lower sensitivity to quantization. Fig. 10 strengthens our former observations of Fig. 9 (when q = 100): the non-scaling noticed for healthy and atrial fibrillation subjects is persistent as q increases. Also, Fig. 10 confirms that the scaling property of subjects with congestive heart failure is the closest to the reference statistically self-similar process BMMT, even though a slight non-scaling remains visible for these subjects. We verified that the same results hold if we perform the same calculations with the Legendre spectra τ * n instead of large deviations ones. From a physiological viewpoint, we may interpret the absence of scaling behavior in RR signals recorded from control subjects, as the ability of heart rate to vary over a wide range of amplitude in a very short time, but to stabilize around mean values over larger periods. That is coherent with a large deviations spectrum whose roughness support narrows as the scale gets coarser. This is particularly true for the left side of the spectra, which corresponds to relatively large heart rate oscillations over small time scales. This phenomenon is amplified in atrial fibrillation which exhibit a richer HRV. Conversely, scale invariant spectra, as in the case of congestive heart failure, may be the sign of a deficient heart with uniform variability at all scales. Only an adaptive non-parameterized estimator providing a sequence of scale-dependent multifractal spectra that do not presuppose scale invariance, is able to reveal such subtle non-scaling behaviors.

More generally though, we believe that large deviations spectra, along with the efficient estimation procedure we propose, is a right tool to determine whether scale invariance properties hold true or not for a sampled finite size time series.

Conclusion

In this paper, we revisited the multifractal analysis of heart inter-beat intervals from the angle of large deviations theory. We proposed an adaptive non-parametric algorithm to estimate the large deviations spectrum of a (a) Surface under the whole spectrum Surface Figure 10: Reinforced scaling/non-scaling criterion for RR signals and Brownian motion in multifractal time: "slope" of the mapping n ∈ [n min , n max ] → A q (n)/A q (n max ), as a function of q ∈ (0, 100]. The scales interval for the computation of the slope is [n min , n max ] = [START_REF] Ferrario | Comparison of entropy-based regularity estimators: Application to the fetal heart rate signal for the identification of fetal distress[END_REF][START_REF] Signorini | Linear and nonlinear parameters for the analysis of fetal heart rate signal from cardiotocographic recordings[END_REF]. The surface A q (n)/A q (n max ) corresponds to: (a) A q (n) = A -q,q (n): surface under the spectrum restricted to [τ n (q), τ n (-q)] -(b) A q (n) = A 0,q (n): surface under the spectrum restricted to [τ n (q), τ n (0)] -(c) A q (n) = A -q,0 (n): surface under the spectrum restricted to [τ n (0), τ n (-q)]. The different curves correspond to the mean and the error bar to the std, computed on the 5 signals of each classe (BMMT, healthy, AF, CHF).

series at each scale. Combined with our original estimator, the advantage of the large deviations spectrum is twofold. Firstly, by essence and compared to Legendre estimates, it allows exhibiting richer scaling structures, such as non concavities reminiscent of a depletion or of an overexpression of roughness exponents within a given interval. Secondly, because it does not assume any a priori scale invariance in the data, it permits to a posteriori validate the effectiveness of such property. We firstly demonstrated the accuracy of the proposed estimation procedure and its sensitivity to non-concavities on toy examples based on the well-known Brownian motion in multifractal time. Then we showed that in RR signals, such non-concavities may relate to physiological factors such as the presence of extrasystoles. We finally provided a way to measure the scale invariance of a signal, which applied to RR signals shows that non-scaling seems to be the norm for control subjects, and that the reacher is the HRV, the lesser the signal presents scaling.

Regarding analytical processes to model HRV time series, so far multifractal analyses implicitly considered that scale invariance was an inherent property of RR signals. Our preliminary study seems to challenges this common thought and refocuses on scale invariance property as a possible discriminant indicator between population groups according to diverse physiological features or experimental conditions.

While the proposed method to analyze the HRV characteristics shows very promising results, it appears necessary to complement our preliminary study by a thorough clinical analysis of larger groups of subjects to assess the significance of our results and understand their limits.

Finally, it would be interesting to apply the deviation theory used in this paper to the Multifractal Detrended Cross Correlation Analysis considered in [START_REF] Zhou | Multifractal detrended cross-correlation analysis for two nonstationary signals[END_REF] after the introduction of the Detrended Cross Correlation Analysis in [START_REF] Podobnik | Detrended cross-correlation analysis: A new method for analyzing two non-stationary time series[END_REF][START_REF] Podobnik | Crosscorrelations between volume change and price change[END_REF], in order to analyse the fine scale invariance properties of partition functions build on the correlations of two non-stationary series.

Figure 1 :

 1 Figure 1: Brownian motion in multifractal time: (a) ordinary Brownian motion B -(b) multifractal measure M obtained by a random dyadic cascade with weights 0.2 and 0.8 of equal probability -(c) Brownian motion B in multifractal time M . The signals are generated at scale J = 18 (they have 2 18 points).

Fig. 3 -Figure 2 :

 32 Figure 2: Oscillation variability, and its transformation into roughness exponent's variability for the Brownian motion in multifractal time at scales n = 12 and n = 8: (a) oscillations -(b) oscillation-based roughness exponents. The signal is generated at scale 18.

Figure 3 :

 3 Figure 3: Large deviations spectra for the Brownian motion in multifractal time: (a) with the oscillation-based roughness exponent -(b) with the MDFA-based roughness exponent -(c) with the wavelet-based roughness exponent. The signal is generated at scale 18 and estimations are presented at scales n = 8 (black) to n = 13 (red). The blue curve is the theoretical spectrum.

Figure 4 :

 4 Figure 4: Large deviations spectra estimates of the signal of Fig. 1-(c) with increments re-quantized on 100 values (scales n = 8 (black) to n = 13 (red)). The blue curve is the theoretical spectrum.

Figure 5 :

 5 Figure 5: Sensitivity of the estimation's algorithm to brutal variations: (a) signal of Fig. 1-(a) with 3 added "extrasystoles" -(b) corresponding large deviations spectrum estimation (red) and Legendre spectrum estimation (black) at scale n = 10. The blue curve is the theoretical spectrum (of the signal without peaks). -(c) signal of Fig. 1-(a) with 30 added "extrasystoles" of two distinct size (5 big and 25 small) -(d) corresponding large deviations spectrum estimation (red), and Legendre spectrum estimation (black) at scale n = 10. The blue curve is the theoretical spectrum (of the signal without peaks).

Figure 6 :

 6 Figure 6: Example of non-concave spectrum on the RR signal of a subject with heart failure (2 13 beats): (a) the RR signal containing about ten series of few consecutive extrasystoles -(b) large deviations spectra estimations with the extrasystoles (red) and Legendre spectrum estimation (black) -(c) large deviations spectra estimations without the extrasystoles (magenta) and Legendre spectrum estimation (black). The spectra are represented at scale n = 11. The normalization used to obtain the spectra of (b) and (c) are derived from the signal with and without extrasystoles respectively.

Figure 7 :

 7 Figure 7: Analysis of the scale invariance hypothesis on the RR signal of a control subject (2 13 beats): (a) the RR signal -(b) large deviations spectra estimations at scales n = 12 to n = 7 (red to black), Legendre spectrum estimated with regression in the range n = 8 to n = 11 (blue) -(c) functions n • τ X,n (q) vs scale n.

Figure 9 :

 9 Figure 9: Scaling/non-scaling criterion for RR signals and Brownian motion in multifractal time: surface beneath the spectrum vs scale n. The surface is computed: (a, whole spectrum) on [α min (n), α max (n)] -(b, left part of the spectrum) on [α min (n), τ n (0)] and (c, right part of the spectrum) on [τ n (0), α max (n)]. The different curves correspond to the mean and the error bar to the std, computed on the 5 signals of each classe (BMMT, healthy, AF, CHF).

  RR controls stat. on RR AF stat. on RR CHF stat. on BMMTs q (b) over [0, q] (c) over [-q, 0]

Table 1 :

 1 Scales of interest. Numbers of dyadic intervals are given for the resolution J = 13, which is the standard value generally used in the following analyses. For RR signals, this value corresponds to an approximate duration of 2h15min (at 60 bpm).

	relative	nb. of RR	nb dyad.	physiol.
	scale	in dyad. int.	int.	oscillation
	J -1	3	4096
	J -2	5	2048	respiration
	J -3	9	1024	baroreflex
	J -4	17	512
	J -5 J -6	33 65	256 128	chemoreflex

http://www.physionet.org/challenge/chaos/

$ This work was supported by the French national research agency (ANR) project DMASC, coordinated by J. Barral.

tions spectrum. The previous discussion of the non-scaling/scaling leads to seriously question the assumption often made that RR-signals possess the scale invariance property. At least, one has to consider this notion as relative, and our surface computations provide a way to compare the respective scale invariance of two signals. However, we considered the whole spectra (and their left and right sides) estimated for q ∈ [-100, 100] (q ∈ [0, 100] and q ∈ [-100, 0] respectively), whereas in most related works scaling is assumed in smaller q's intervals (typically q ∈ [-5, 5]). Consequently, to be more convincing, our calculation should provide the same comparison of surfaces for spectra corresponding to the ranges [-q, q], [0, q] and [-q, 0] for any q > 0, and show their evolution as q increases. Specifically, the scaling property is now measured as follows. For each q, we define the mapping A -q,q : n → A -q,q (n). To a scale n in an interval [n min , n max ], it associates the surface A -q,q (n) under the scale-n spectrum restricted to [τ n (q), τ n (-q)]. We define similarly the mappings A 0,q and A -q,0 corresponding to the spectrum restricted to [τ n (q), τ n (0)] and [τ n (0), τ n (-q)] respectively. (These mappings are given for q = 100 in Fig. 9.) From each mapping A q ∈ {A -q,q , A 0,q , A -q,0 }, we then obtain a measure of the


Florian Lemarié
email: florian@atmos.ucla.edu

Laurent Debreu
email: laurent.debreu@imag.fr

Eric Blayo
email: eric.blayo@imag.fr

Optimal control of the convergence rate of Schwarz Waveform Relaxation algorithms

In this study we present a non-overlapping Schwarz waveform relaxation method applied to the one dimensional unsteady diffusion equation. We derive efficient interface conditions using an optimal control approach once the problem is discretized. Those conditions are compared to the usual optimized conditions derived at the PDE level by solving a min-max problem. The performance of the proposed methodology is illustrated by numerical experiments.

Introduction

Schwarz-like domain decomposition methods are very popular in mathematics, computational sciences, and engineering notably for the implementation of coupling strategies. This type of method, originally introduced for stationary problems, can be extended to evolution problems by adapting the waveform relaxation algorithms to provide the so-called Schwarz waveform relaxation method [START_REF] Gander | Optimal convergence for overlapping and non-overlapping Schwarz waveform relaxation[END_REF][START_REF] Gander | Methodes de relaxation d'ondes pour l'equation de la chaleur en dimension 1[END_REF]. The idea behind this method is to separate the spatial domain, over which the time-evolution problem is defined, into subdomains. The resulting time-dependent problems are then solved separately on each subdomains. An iterative process with an exchange of boundary conditions at the interface between the subdomains is then applied to achieve the convergence to the solution of the original problem. To accelerate the convergence speed of the iterative process, it is possible to derive efficient interface conditions by solving an optimization problem related to the convergence rate of the method [e.g. ; 1, 5].

In this study, we specifically address the optimization problem arising from the use of Robin type transmission conditions in the framework of a non-overlapping Schwarz waveform relaxation. For this type of problem, the existing work has been achieved mainly at the PDE level, giving rise to the optimized Schwarz waveform relaxation algorithm [START_REF] Bennequin | Optimized Schwarz waveform relaxation methods for convection reaction diffusion problems[END_REF][START_REF] Gander | Optimal convergence for overlapping and non-overlapping Schwarz waveform relaxation[END_REF][START_REF] Gander | Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF]. The objective here is to use the optimal control theory paradigm [START_REF] Lions | Contrôle optimal des systèmes gouvernés par des équations aux dérivées partielles[END_REF] to find parameters optimized at the discrete level, and thus to systematically make a comparison with the parameters determined at the PDE level. This paper is organized as follows : in section 2 we briefly recall the basics of optimized Schwarz methods in the framework of a time evolution problem. Section 3 is dedicated to the determination of the optimal control problem that we intend to address. Finally, in section 4 we apply our approach to a diffusion problem.

L u = f , in Ω × [0, T ], (1)
Bu = g, on ∂Ω × [0, T ], (2)
where L and B are two partial differential operators, and f the forcing. This problem is complemented by an initial condition

u(x, 0)=u 0 (x), x ∈ Ω . (3)
We consider a splitting of the domain Ω into two non-overlapping domains Ω 1 and Ω 2 communicating through their common interface Γ . The operator L introduced previously is split into two operators L j restricted to Ω j ( j = 1, 2). By noting F 1 , F 2 , G 1 and G 2 the operators defining the interface conditions, the alternating form of the Schwarz waveform relaxation algorithm reads

       L 1 u k 1 = f 1 , in Ω 1 × [0, T ], u k 1 (x, 0)=u o (x), x ∈ Ω 1 , B 1 u k 1 (x,t)=g 1 , in [0, T ] × ∂Ω 1 , F 1 u k 1 (0,t)=F 2 u k-1 2 (0,t), in Γ × [0, T ],        L 2 u k 2 = f 2 , in Ω 2 × [0, T ], u k 2 (x, 0)=u o (x), x ∈ Ω 2 , B 2 u k 2 (x,t)=g 2 , in [0, T ] × ∂Ω 2 , G 2 u k 2 (0,t)=G 1 u k 1 (0,t), in Γ × [0, T ], (4)
where k = 1, 2,... is the iteration number, and the initial guess u 0 2 (0,t) must be given. The operators F j and G j must be chosen to impose the desired consistency of the solution on the interface Γ . We consider here the one-dimensional diffusion equation with constant (possibly discontinuous) diffusion coefficients κ j (κ j > 0, j = 1, 2). We define

L j = ∂ t -κ j ∂ 2 x , Ω 1 =(-L 1 , 0), Ω 2 =(0, L 2 ) (L 1 , L 2 ∈ R + )
, and Γ = {x = 0}. In this context, we require the equality of the subproblems solutions and of their normal fluxes on the interface Γ ,

u 1 (0,t)=u 2 (0,t), κ 1 ∂ x u 1 (0,t)=κ 2 ∂ x u 2 (0,t), t ∈ [0, T ]. (5)
To obtain such a consistency we use mixed boundary conditions of Robin type

F j = -κ j ∂ x + p 1 , G j = κ j ∂ x + p 2 , ( j = 1, 2),
where p 1 and p 2 are two parameters that can be optimally chosen to improve the convergence speed of the Schwarz method. Algorithm (4) with two-sided Robin conditions (i.e. for p 1 = p 2 ) is well-posed for any choice of p 1 and p 2 such that p 1 + p 2 > 0. This result can be shown using a priori energy estimates, as described in [START_REF] Gander | Methodes de relaxation d'ondes pour l'equation de la chaleur en dimension 1[END_REF].

Optimization of the convergence factor

To demonstrate the convergence of algorithm (4) a classical approach [e.g. 6] is to define the error e k j between the exact solution u and the iterates u k j . A Fourier analysis enables the transformation of the original PDEs into ODEs that can be solved analytically. The analytical solution on each subdomain is then used to define a convergence factor ρ of the corresponding Schwarz algorithm. For a diffusion problem, defined on subdomains of infinite size (i.e. assuming L 1 , L 2 → ∞), we get

ρ(p 1 , p 2 , ω)= (p 2 - √ iωκ 2 ) (p 2 + √ iωκ 1 ) (p 1 - √ iωκ 1 ) (p 1 + √ iωκ 2 ) , (6)
where p 1 and p 2 are two degrees of freedom which can be tuned to accelerate the convergence speed. In ( 6), i = √ -1, and ω ∈ R is the angular frequency arising from a Fourier transform in time on e k j . A general approach to choose the Robin parameters p 1 and p 2 is to solve a minimax problem [START_REF] Gander | Optimal convergence for overlapping and non-overlapping Schwarz waveform relaxation[END_REF] min

p 1 ,p 2 ∈R max ω∈[ω min ,ω max ] ρ(p 1 , p 2 , ω) . ( 7
)
Because we work in practice on a discrete problem the frequencies allowed by the temporal grid range from ω min = π/T to ω max = π/∆t, where ∆t is the time step of the temporal discretization. For the diffusion problem under consideration here, the analytical solution of the optimization problem [START_REF] Lemarié | Optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients[END_REF] has been derived in [START_REF] Lemarié | Toward an optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients, part 1 : the constant coefficients case[END_REF] in a general two-sided case (i.e. with p 1 = p 2 ) with discontinuous coefficients

κ 1 = κ 2 .
For the sake of simplicity, we consider in the present study the continuous case (κ 1 = κ 2 = κ) and we recall the result found in [START_REF] Lemarié | Toward an optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients, part 1 : the constant coefficients case[END_REF] in this case.

Theorem 1. Under the assumption κ 1 = κ 2 = κ, the optimal parameters p 1 and p 2 of the minmax problem [START_REF] Lemarié | Optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients[END_REF] are given by

p 1 = α √ 2κ 4 8 + v 2 -v , p 2 = α √ 2κ 4 8 + v 2 + v , where α =(ω min ω max ) 1/4 , β = α -1 ( √ ω min + √ ω max ) and v =    2 β -1 if β ≥ 1 + √ 5, 2β 2 -12 if √ 6 ≤ β < 1 + √ 5, 0 if 2 < β < √ 6.
It is worth mentioning that even if the diffusion coefficients are continuous the twosided case provides a faster convergence than the one-sided case studied in [START_REF] Gander | Methodes de relaxation d'ondes pour l'equation de la chaleur en dimension 1[END_REF] (Fig. 1).

General remarks :

• The usual methodology to optimize the convergence at the continuous level comes with a few assumptions that may lead to inaccuracies once the problem is discretized. For example, as discussed in [START_REF] Lemarié | Optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients[END_REF] (Sec. 5), the infinite domain assumption used to determine the convergence factor (6) may lead to appreciable differences in the optimized parameters compared to an approach taking the finiteness of the subdomains into account. We numerically found that the infinite domain assumption is valid as long as the dimensionless Fourier number Fo = κ j /(L 2 j ω) (with L j the size of subdomain Ω j ) of the problem does not exceed a critical value Fo c = 0.02.

• The optimization problem [START_REF] Lemarié | Optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients[END_REF] aims at minimizing the maximum value of ρ(p 1 , p 2 , ω) over the entire interval [ω min , ω max ]. This provides a very robust method general enough to deal with the worst case scenario when all the temporal frequencies are present in the error. An even more efficient way to proceed would be to adjust the values of p 1 and p 2 at each iteration so that those parameters are efficiently chosen to "fight" the remaining frequencies in the error. Convergence factor optimized at the PDE level in the one-sided case (black line) [START_REF] Gander | Methodes de relaxation d'ondes pour l'equation de la chaleur en dimension 1[END_REF] and in the two-sided case (dashed black line) [START_REF] Lemarié | Toward an optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients, part 1 : the constant coefficients case[END_REF], for κ = 10 -2 ms -1 , ∆t = 10 s, and T = 2 13 ∆t.

Optimal control of the Robin parameters

To investigate the robustness of the optimized parameters once the problem is discretized, the use of the optimal control theory appears as a natural choice. We aim at controlling the Robin parameter in order to get the best possible convergence speed in the sense of a given cost function J . Moreover, following the approach of [START_REF] Gander | A non-overlapping optimized Schwarz method which converges with arbitrarily weak dependence on h[END_REF] and the previous discussion, we consider the possibility to use different parameters p j for different steps of the iterative process. It is easy to check that by choosing different parameters at each iteration we still converge to the solution of the global problem. A first way to choose the parameters is to look, at each iteration k, for p k 1 and p k 2 minimizing the error at the interface. In this case the cost function that we intend to minimize at each iteration would be

J (p k 1 , p k 2 )= w 2 T 0 u k 1 (0,t) -u k 2 (0,t) 2 dt + w 2 T 0 κ 1 ∂ x u k 1 (0,t) -κ 2 ∂ x u k 2 (0,t) 2 dt. (8)
The constants w and w must be chosen to balance both terms, depending on the characteristics of the problem (see Sec. 4). The cost function ( 8) is designed in agreement with the consistency (5) we want to impose at the interface between subdomains. J provides a measure of the "inconsistency" of the solution at each iteration k, and is, thus, directly related to the order of magnitude of the errors e k j of the algorithm (as shown in Fig. 2). An other strategy could be to minimize the error at a given iteration K. The cost function would thus be

J (p k 1 , p k 2 ) k=1,K = w 2 T 0 u K 1 (0,t) -u K 2 (0,t) 2 dt + w 2 T 0 κ 1 ∂ x u K 1 (0,t) -κ 2 ∂ x u K 2 (0,t) 2 dt, (9)
leading to an optimization on 2K parameters. This latter approach is particularly interesting when we intend to obtain the best possible approximation of the exact solution after a number of iterations set in advance. We propose here to lead our study with this kind of approach with K = 5. The optimal control approach does not per se reduce the computational cost of the algorithm because many evaluations of the cost function are required during the minimization process (see algorithm 1). We use this approach as a tool to improve our understanding of the behavior of the Robin parameters in order to find new directions to further accelerate the convergence speed when Robin-type interface conditions are used. We denote by p ,num 1 and p ,num 2 the parameters found numerically by solving the optimal control problem. Those parameters correspond to two vectors of size K. Similarly we will denote by p ,ana 1 and p ,ana 2 the parameters found analytically (cf Theorem 1). We used Matlab for the computation (algorithm 1). Note that the well-posedness of the coupling problem (4) is not sufficient to ensure a well-posed optimal control problem. Some additional requirements on the convexity and regularity of the cost function are necessary. We do not provide here such a proof, however we empirically checked that the same solution of the optimal problem is obtained for a wide range of parameter values for the initial guess.

Algorithm 1 Optimal control %== Robin parameters found analytically : p1ana, p2ana %== Solution of the optimal control problem : p1opt, p2opt %== Initial guess ==% x 0 (1:2:2*K-1)=p1ana; x 0 (2:2:2*K )=p2ana; %== Solve the optimal control problem ==% %== the CalcJ function proceeds to K iterations of the %== Schwarz algorithm using 2K Robin parameters, %== and computes the associated cost function [START_REF] Lions | Contrôle optimal des systèmes gouvernés par des équations aux dérivées partielles[END_REF] x = fminsearch( @CalcJ, x 0 ); %== Retrieve the optimized parameters p1opt(1:K)=x(1:2:2*K-1); p2opt(1:K)=x(2:2:2*K );

Numerical experiments

We discretized problem (4) using a backward Euler scheme in time and a second order scheme defined on a staggered grid in space (see [START_REF] Lemarié | Toward an optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients, part 1 : the constant coefficients case[END_REF] for more details). We decompose the domain Ω into two non-overlapping subdomains Ω 1 =[-H, 0] and Ω 2 =[0, H] with H = 500 m. The diffusion coefficient is κ = 10 -2 m 2 s -1 and the total simulation time is T = 2 13 ∆t with ∆t = 10 s. The parameter values lead to a dimensionless Fourier number smaller than 0.02 so that the infinite domain assumption is valid. We simulate directly the error equations, i.e. f 1 = f 2 = 0 in (4) and u 0 (x)=0. We start the iteration with a random initial guess u 0 2 (0,t)(t ∈ [0, T ]) so that it contains a wide range of the temporal frequencies that can be resolved by the computational grid. This is done to allow a fair comparison as the parameters optimized at the PDE level are optimized assuming that the full range [ω min , ω max ] is present in the error. We first perform the Optimized non-overlapping Schwarz Method (referred as to OSM case) using p ,ana 1 and p ,ana 2 and then using an optimal control of the Robin parameters with K = 5 (referred as to OptCon case). We first check that the minimization of cost function J consistently implies the reduction of the errors e j ∞ of the associated algorithm (Fig. 2). For our experiments, we chose w = 1 and w = H/κ in [START_REF] Lions | Contrôle optimal des systèmes gouvernés par des équations aux dérivées partielles[END_REF]. We notice that in the OptCon case the convergence speed is significantly improved compared to the OSM case. Indeed, 9 iterations of the OSM are required to obtain the same accuracy than the OptCon case after only 5 iterations. In order to have more insight on the way the parameters p ,num 1 and p ,num 2 evolve throughout the iterations we plot, in Fig. 3, the corresponding convergence factor (6) at each iteration. It is striking to realize that the optimal convergence is obtained through a combination of 2-point (equivalent to the one-sided case) and 3point (equivalent to the two-sided case) equioscillations sometimes shifted along the ω-axis to adapt to the temporal frequencies still present in the error. The first two iterations aim at working mainly on the high-frequency components while the last three iterations are optimized to work on the low-frequency component. The adap-tivity of the Robin parameters from one iteration to the other brings more flexibility to the method enabling more scale selectivity.

Conclusion

Due to its simplicity, the use of Robin-type transmission conditions is very attractive when one wants to couple unsteady problems defined on non-overlapping subdomains. Once the Robin parameters are properly chosen one can achieve a fast convergence [START_REF] Gander | Optimal convergence for overlapping and non-overlapping Schwarz waveform relaxation[END_REF]. In the present study we showed that there is still room for improvement in the design of the Robin conditions. If the Robin parameters are adjusted from one iteration to the other we showed, thanks to an optimal control approach, that we can significantly improve the convergence speed. It is important to emphasize that the optimal control paradigm proposed in this study is general enough to be used with any type of PDE and an arbitrary number of subdomains.

  Fig.1. Convergence factor optimized at the PDE level in the one-sided case (black line)[START_REF] Gander | Methodes de relaxation d'ondes pour l'equation de la chaleur en dimension 1[END_REF] and in the two-sided case (dashed black line)[START_REF] Lemarié | Toward an optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients, part 1 : the constant coefficients case[END_REF], for κ = 10 -2 ms -1 , ∆t = 10 s, and T = 2 13 ∆t.

Fig. 2 .

 2 Fig. 2. Evolution of the L ∞ -norm of the error (left) and of the cost function J (right) with respect to the iterates k in the OSM and OptCon cases.

Fig. 3 .

 3 Fig. 3. Sequence of convergence factors ρ(ω) resulting from the optimal control of the Robin parameters determined to get the best possible convergence after K = 5 iterations.

2

  Optimization of the convergence at the PDE level 2.1 Model problem and Optimized Schwarz Methods Let us consider Ω a bounded open set of R. The model problem is to find u such that u satisfies over a time period [0, T ]

Acknowledgement. This research was partially supported by the ANR project COMMA (COupling in Multi-physics and multi-scale problems: Models and Algorithms) and by the INRIA project-team MOISE (Modelling, Observation and Identification for Environmental Sciences). We are thankful to Héloïse Pelen (ENS Lyon) for her contribution during her masters internship.