Oliver Orejola
email: oorejola@tulane.edu

Gustavo Didier
email: gdidier@tulane.edu

Patrice Abry
email: patrice.abry@ens-lyon.fr

Herwig Wendt
email: herwig.wendt@irit.fr

Hurst multimodality detection based on large wavelet random matrices

Keywords: self-similarity, operator fractional Brownian motion, wavelets, random matrices, high dimensions

In the modern world, systems are routinely monitored by multiple sensors, generating "Big Data" in the form of a large collection of time series. In this paper, we put forward a statistical methodology for detecting multimodality in the distribution of Hurst exponents in high-dimensional fractal systems. The methodology relies on the analysis of the distribution of the log-eigenvalues of large wavelet random matrices. Depending on the presence of a single or many Hurst exponents, we show that the wavelet empirical log-spectral distribution displays one or many modes, respectively, in the threefold limit as dimension, sample size and scale go to infinity. This allows for the construction of a unimodality test for the Hurst exponent distribution. Monte Carlo simulations show that the proposed methodology attains satisfactory power for realistic sample sizes.

I. INTRODUCTION

Context: scale invariance. Scale invariance has been observed in signals coming from a wide range and variety of contexts in physics and engineering, see, e.g., [START_REF] Sornette | Critical Phenomena in Natural Sciences: Chaos, Fractals, Selforganization and Disorder: Concepts and Tools[END_REF], [START_REF] Abry | Wavelet analysis of long-range-dependent traffic[END_REF]. A signal X is called scale invariant, or fractal, when its temporal dynamics lack a characteristic scale. The focus of the analysis is thus on identifying the scaling exponents that relate the continuum of scales [START_REF] Flandrin | Wavelet analysis and synthesis of fractional Brownian motion[END_REF]. A cornerstone model of scale invariance is selfsimilarity: X is called self-similar when its finite-dimensional distributions (f.d.d.) are invariant under suitable time scaling, i.e., {X(t)} t∈R f.d.d.

= {a H X(t/a)} t∈R , a > 0, where the scaling exponent 0 < H < 1 is called the Hurst parameter. The most prominent example is fractional Brownian motion (fBm), the only Gaussian, self-similar process with stationary increments [START_REF] Pipiras | Long-Range Dependence and Self-Similarity[END_REF]. Estimation of the Hurst parameter H plays a key role in signal processing tasks such as characterization, diagnosis, classification and detection. The so-named wavelet transform provides the analytical basis for well-established estimation methodologies for H [START_REF] Veitch | A wavelet-based joint estimator of the parameters of long-range dependence[END_REF]. Challenge: high-dimensional time series. The modeling of self-similarity in applications has remained so far based on the univariate fBm model. Yet, in the modern world of "Big Data," a plethora of sensors monitor natural and artificial systems, generating large data sets in the form of several joint G.D.'s long term visits to ENS de Lyon were supported by the school, the CNRS and the Simons Foundation collaboration grant #714014. H.W. supported by ANR-18-CE45-0007 MUTATION, France. time series. This can be seen in many fields of application. In neuroscience, for example, the number of macroscopic brain activity time series ranges from hundreds (MEG data) to several tens of thousands (fMRI data) [START_REF] Ciuciu | Scale-free and multifractal time dynamics of fMRI signals during rest and task[END_REF]. Likewise, in climate studies, dealing with large numbers of measured components has become standard [START_REF] Isotta | The climate of daily precipitation in the Alps: development and analysis of a high-resolution grid dataset from pan-Alpine raingauge data[END_REF]. For such high-dimensional data, a multitude of scaling laws -i.e., Hurst multimodality -implies different large scale behavior of the system along possibly non-canonical coordinate axes. Ignoring Hurst multimodality in data may lead to arbitrarily large estimation biases due to the so-named dominance and amplitude effects (see [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF]). Related work: self-similarity in high dimensions. A recently proposed multivariate extension of fBm, operator fractional Brownian motion (ofBm), [START_REF] Mason | Sample path properties of operator-selfsimilar Gaussian random fields[END_REF]- [START_REF] Didier | Exponents, symmetry groups and classification of operator fractional Brownian motions[END_REF] frames the modeling of scale invariance in a high-dimensional time series setting [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF]. The so-named wavelet eigenanalysis methodology [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF] was shown to lead to efficient and robust estimation of Hurst exponents in both multivariate (fixed dimensions) and highdimensional (the dimension p(n) of the measurements grows as a function of the time series size n) settings [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF], [START_REF] Abry | Detecting and estimating multivariate self-similar sources in high-dimensional mixtures[END_REF], [START_REF] Boniece | On multivariate non-Gaussian scale invariance: fractional Lévy processes and wavelet estimation[END_REF]. Use of this model in applications immediately leads to a first critical question: how many different scaling laws exist in the possibly very large number of time series? In the statistical signal processing literature, the problem of identifying the number and properties of sources in multivariate or highdimensional noisy signals has been studied for decades [START_REF] Comon | Handbook of Blind Source Separation: Independent Component Analysis and Applications[END_REF]- [START_REF] Da Costa | Comparison of model order selection techniques for high-resolution parameter estimation algorithms[END_REF]. Examples of the proposed techniques include principal component analysis, factor analysis and sparse graphical Gaussian models [START_REF] Johnson | Applied Multivariate Statistical Analysis[END_REF]. Nevertheless, there has been a paucity of estimation methodologies for both high-dimensional and scale invariant signals; see a contrario [START_REF] Zhang | Identifying cointegration by eigenanalysis[END_REF] or [START_REF] Lucas | Counting the number of different scaling exponents in multivariate scale-free dynamics: Clustering by bootstrap in the wavelet domain[END_REF], proposing a bootstrap-based method for counting the number of distinct Hurst exponents based on wavelet eigenanalysis, yet in a multivariate context with modest dimension. A key related difficulty is the study of random matrices under dependence, as they emerge in wavelet eigenanalysis, which is still a very active area of research [START_REF] Liu | On the Marčenko-Pastur law for linear time series[END_REF], [START_REF] Erdős | Random matrices with slow correlation decay[END_REF]. This is even more so in regards to the presence of fractional memory [START_REF] Merlevède | On the empirical spectral distribution for matrices with long memory and independent rows[END_REF]. Goal, contributions and outline. The goal of the present work is to study, in the context of high-dimensional asymptotic limits (as the sample size, the number of components and the scale go to infinity at a fixed rate), the properties of the statistical distribution of the multivariate eigenwavelet-based estimation of the vector of Hurst exponents, both theoretically and practically. Notably, aiming to test the equality of all Hurst exponents, the uni-vs. multi-modality of these distributions is carefully addressed. To that end, the definitions and properties of ofBm and the high-dimensional model are summarized in Sections II-A and II-B, while multivariate eigenwavelet estimation procedures are briefly sketched in Section II-C. Section III describes the core contribution of this work. It explains, first, how large random matrix theory leads to the modification of the high-dimensional limit formulation by incorporating the analysis of scales and establishes, second, the high-dimensional asymptotic behavior of the distribution of the estimator of the vector of Hurst exponents, thus permitting the discussion of uni-vs. multimodality. Furthermore, these theoretical asymptotic limits are tested in practice from extensive Monte Carlo simulations using synthetic ofBm sample paths, see Section IV. First, the uni-or multimodality of the distributions of estimators are shown to be practically observed with data of moderate and realistic sizes (see Section IV-B). Second, elaborating on the Hartigans' dip test statistic [START_REF] Hartigan | The dip test of unimodality[END_REF], we devise a procedure used to test unimodality in the distributions of the estimated Hurst exponents, and assess its practical power and relevance (see Section IV-B). Reported results indicate that the proposed test has satisfactory performance and can be readily applied to real-world high-dimensional data, including for sample sizes typical in neuroscience (cf. [START_REF] La Rocca | Revisiting functional connectivity for infraslow scale-free brain dynamics using complex wavelets[END_REF]).

II. SELF-SIMILARITY ANALYSIS AND MODELING A. Operator fractional Brownian motion

Operator fractional Brownian motion (ofBm) is a canonical model for multidimensional scale invariant structures in realworld data. We briefly recall its definition and some properties (see [START_REF] Didier | Integral representations and properties of operator fractional Brownian motions[END_REF] for the general definition and properties of ofBm).

Let B H,Σ (t) = B H1 (t), . . . , B Hp (t) t∈R denote a collection of p possibly correlated fBm components defined by their individual self-similarity exponents H = (H 1 , . . . , H p ), 0 < H 1 ≤ . . . ≤ H p < 1. Let Σ be a pointwise covariance matrix with entries (Σ) ℓ,ℓ ′ = σ ℓ σ ℓ ′ ρ ℓ,ℓ ′ , where σ 2 ℓ are the variances of the components and ρ ℓ,ℓ ′ their (pairwise) correlation coefficients. We define ofBm as the stochastic process B P,H,Σ (t) := P B H,Σ (t), where P is a real-valued, p×p invertible matrix that mixes the components (changes the scaling coordinates) of B H,Σ (t). OfBm consists of a multivariate Gaussian self-similar process with stationary increments. Moreover, it satisfies the (operator) self-similarity relation

{B P,H,Σ (t)} t∈R f.d.d. = {a H B P,H,Σ (t/a))} t∈R , (1)
∀a > 0. In (1), the matrix (Hurst) exponent is given by H = P diag(H)P -1 , and a H := = stands for the equality of finite-dimensional distributions.

B. High-dimensional model

Let π(dH) be a discrete distribution of Hurst exponents with (ordered) support {H 1 , . . . , H m }, m ∈ N. Given a vector H ∈ (0, 1) p of i.i.d. samples from π(dH), the process

Y (t) := B P,H,Σ (2)
as defined in Section II-A is, conditionally on H, a pvariate ofBm with Hurst matrix H. Given a time series {Y (t)} t=1,...,n , we further assume

p = p(n) → ∞ as n → ∞,
under which ( 2) is a high-dimensional model and has p = p(n) → ∞ Hurst exponents (on models of the general form (2) under weak dependence, see, for instance, [START_REF] Chang | Principal component analysis for secondorder stationary vector time series[END_REF]).

C. Scaling in the wavelet domain

Multivariate wavelet transform. Let ψ be a mother wavelet, i.e., a real-valued function such that R ψ 2 (t)dt = 1. For all k, j ∈ Z, the multivariate DWT of

{Y (t)} t∈R is de- fined as D Y (2 j , k) := (D Y1 (2 j , k), . . . , D Yp (2 j , k)), where D Y ℓ (2 j , k) := ⟨2 -j/2 ψ(2 -j t -k)|Y ℓ (t)⟩ ∈ R for ℓ ∈ {1, . . . , p}.
For a detailed introduction to wavelet transforms, see [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. It can be shown that the wavelet coefficients {D Y (2 j , k)} k∈Z of p-variate ofBm Y = B P,H,Σ satisfy, for every fixed octave j, the operator self-similarity relation [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF]

{D Y (2 j , k)} k∈Z f.d.d. = {2 j(H+ 1 2 I) D Y (1, k)} k∈Z . (3)
Wavelet random matrices and high-dimensional eigenanalysis. Given any p-variate process {Y (t)} t∈R (in particular, model ( 2)), the sample wavelet spectrum (variance) at octave j = j 1 , . . . , j 2 is given by the p × p wavelet random matrices

S Y (2 j ) = 1 n j nj k=1 D Y (2 j , k)D Y (2 j , k) * (4)
where n is the time series (sample) size and n j ≃ n/2 j is the number of wavelet coefficients available at scale 2 j . It has been shown and discussed [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF] that, in general, estimation based on the entrywise multiscale behavior of S Y (2 j ) is arbitrarily biased and effectively meaningless. Instead, it was proposed in [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF] that the estimation of the multivariate Hurst exponents should be based on the eigenvalues λ 1 (2 j ), . . . , λ p (2 j ) of the random matrix S Y (2 j ) as in (4). Notably, it was shown that the Hurst exponents H ℓ can be efficiently estimated by means of the weighted wavelet log-eigenvalues

H ℓ = j2 j=j1 w j log 2 λ ℓ (2 j ) 2 - 1 2 , ℓ = 1, . . . , p, (5)
where j w j = 0 and j jw j = 1. In fact, (5) has good statistical performance in both Gaussian and non-Gaussian frameworks [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF], [START_REF] Abry | Detecting and estimating multivariate self-similar sources in high-dimensional mixtures[END_REF], [START_REF] Boniece | On multivariate non-Gaussian scale invariance: fractional Lévy processes and wavelet estimation[END_REF].

III. HIGH-DIMENSIONAL LIMITS

A. High dimensions: Wavelet log-eigenvalue distribution in the three-way limit As indicated in the Introduction, the present work addresses high-dimensional limits. In contradistinction to the classical setting where the sample size goes to infinity for fixed number of components, i.e., n → +∞, p = p 0 , here we consider high-dimensional settings, where both the sample size and the number of components go jointly to infinity, i.e., n, p → +∞. In the context of sample covariance matrices, this is more sharply characterized by the two-way limit lim n→∞ p/n = c ∈ [0, ∞) [START_REF] Bai | Spectral Analysis of Large Dimensional Random Matrices[END_REF].

In the context of multivariate Hurst exponent estimation using wavelet-based linear regressions, as in ( 5), when both n, p → +∞, thorough analysis of scale invariant systems further implies that the range of scales where linear regressions are performed should also grow towards infinity: (j 1 , j 2 ) → +∞. Technically, in such settings, the high-dimensional behavior analyzed by means of large random matrices actually leads us to consider the three-way limit [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF], [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF] n, p, j → +∞,

p n/2 j → c ∈ [0, ∞), (6)
a non-trivial contribution of the present work, both theoretically and numerically.

B. High dimensions: Multimodality of Hurst exponents

Mathematical results in [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF], [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF] imply that, in the threeway limit (6) and under (2), the eigenvalues of the sample wavelet spectrum behave as

log 2 λ ℓ (2 j ) = C ℓ + j • (2H ℓ + 1) + o P (1), (7)
ℓ = 1, . . . , p, for some constant C ℓ , where o P (1) vanishes in probability. In fact, eigenvalue scaling relations of the type ( 7) are shown to hold for high-dimensional signal-plus-noise instances, broader than the noise-free model (2) [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF], [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF]. Further, assume measurements given by the highdimensional model [START_REF] Abry | Wavelet analysis of long-range-dependent traffic[END_REF]. In light of [START_REF] Isotta | The climate of daily precipitation in the Alps: development and analysis of a high-resolution grid dataset from pan-Alpine raingauge data[END_REF], we obtain the key result of this paper. Fundamental property of weighted log-eigenvalues of wavelet random matrices: the distribution of the p = p(n) estimates H ℓ resulting from weighted log-eigenvalues of wavelet random matrices as defined in (5) satisfies, for small ε > 0 and any q = 1, . . . , m, as n → +∞,

# ℓ = 1, ..., p(n) : H ℓ ∈ (H q -ε, H q + ε) p(n) P → π(H q ). (8)
Also, and importantly, we note that the convergence (8) of the wavelet empirical (log)spectral distribution to π(dH) takes place on the condition that n j = n/2 j is large enough w.r.t. p = p(n) so as to prevent significant bias stemming from eigenvalue repulsion (see [START_REF] Tao | Topics in Random Matrix Theory[END_REF]- [START_REF] Abry | Bootstrap for testing the equality of selfsimilarity exponents across multivariate time series[END_REF]), or for the random matrices (4) to have deficient rank (which yields some negatively infinite log-eigenvalues).

This key theoretical result is further inspected practically in Section IV-B and is illustrated in Fig. 1.

IV. PRACTICAL PERFORMANCE ASSESSMENT A. Monte Carlo simulation setting

To assess the practical relevance of the theoretical results stated in Section III, we make use of Monte Carlo simulations based on 1000 of independent realizations of p-variate measurements B P,H,I = P B I,H,I , where P is a randomly chosen orthogonal matrix and {B I,H,I (t)} t∈R made up of independent univariate fBms with n = 2 14 (such sample sizes are realistic, for example, in the context of the analysis of infraslow brain activity [START_REF] La Rocca | Revisiting functional connectivity for infraslow scale-free brain dynamics using complex wavelets[END_REF]). The univariate fBms are generated using the R package somebm (see [START_REF] Huang | somebm: some Brownian motions simulation functions[END_REF]). The pdimensional vector H is obtained by drawing p i.i.d. samples from Unif(H 1 , . . . , H m ), for each realization independently. Here, we limit ourselves to the consideration of m = 1, 2, 4, and 6. The wavelet transformation is generated by means of Mallat's algorithm based on a Daubechies filter with N ψ = 2. In our simulations, we set the dimension to p = 2 6 . The range of regression scales is chosen to be j 1 = 2 to j 2 = 5. Hence, p/n j < 1 for j = j 1 , . . . , j 2 , implying that all the wavelet random matrices used in the regression procedure have full rank.

B. Multimodality in practice

To illustrate the practical relevance of the high-dimensional asymptotic limits in ( 7) and ( 8), Fig. 1 provides an example of the high-dimensional behavior of the wavelet logeigenvalue distribution. In the plots, a unimodal (left) or a fourmodal (right) distribution strikingly emerges in the wavelet eigenspectrum depending on the presence of one or four modes, respectively, in the support of π(dH). To put this remarkable result in context, consider the traditional random matrix theoretic framework of sample covariance matrices. In this case, as p/n → c ∈ [0, ∞) and under unimodality, we observe a combination (a free convolution; see [START_REF] Anderson | An Introduction to Random Matrices[END_REF]) between a Marcenko-Pastur law and the spectral density of the univariate stochastic process whose independent copies form the sample covariance matrix. In particular, the Hurst exponent only appears implicitly in the limit eigenspectrum distribution, which is not available in closed form [START_REF] Merlevède | On the empirical spectral distribution for matrices with long memory and independent rows[END_REF].

C. Unimodality testing

To go beyond visual evidence, we are now interested in practically testing multimodality vs. unimodality in the distribution of Hurst exponents. More precisely, this means testing the null hypothesis H 0 : the distribution π(dH) of Hurst exponents is unimodal, i.e., that the support of π(dH) contains a single, unspecified value in (0, 1). We put forward a test based on the wavelet empirical (log)spectral distribution, i.e., the distribution of the estimates H ℓ . Indeed, in light of ( 8), this distribution is expected to be unimodal under H 0 , and multimodal otherwise. Unimodality test. With this in mind, we propose a unimodality testing methodology that builds upon the original idea behind Hartigan's dip test for multimodality (see [START_REF] Hartigan | The dip test of unimodality[END_REF]). Let

F (x) = 1 p p ℓ=1 1 { H ℓ ≤x} , x ∈ R,
be the empirical distribution function of the weighted wavelet log-eigenvalues H ℓ , ℓ = 1, . . . , p,

as in [START_REF] Veitch | A wavelet-based joint estimator of the parameters of long-range dependence[END_REF]. As in [START_REF] Hartigan | The dip test of unimodality[END_REF], let U be the class of all unimodal distribution functions, and define the statistic

d = inf G∈U sup x∈R | F (x) -G(x)|, (10)
that is, d is computed based on picking the unimodal distribution that provides the best fit to the empirical distribution function (see also [START_REF] Hartigan | The dip test of unimodality[END_REF], Section 4, on the details of how to compute [START_REF] Mason | Sample path properties of operator-selfsimilar Gaussian random fields[END_REF] in practice). The proposed test procedure is defined as reject

H 0 when d > d α , (11)
where the size of the test α ∈ [0, 1] is defined as [START_REF] Lehmann | Testing Statistical Hypotheses[END_REF] α = sup H0∈(0,1)

P d > d α supp π(dH) = H 0 . ( 12
)
The statistic d as in [START_REF] Mason | Sample path properties of operator-selfsimilar Gaussian random fields[END_REF] was computed using the R package diptest (see [START_REF] Maechler | diptest: Hartigan's Dip Test Statistic for Unimodality -Corrected[END_REF]); test significance was set to α = 0.05. Test statistics practical threshold. In regard to the computation of the rejection threshold d α , note that the observations are assumed independent in the original dip test [START_REF] Hartigan | The dip test of unimodality[END_REF], a condition that does not hold for the weighted wavelet logeigenvalues [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF]. For this reason, we simulate a large number of independent p-variate ofBms with Hurst matrix exponent H = diag(H 0 , . . . , H 0 ), H 0 ∈ (0, 1), calculate the corresponding test statistics [START_REF] Mason | Sample path properties of operator-selfsimilar Gaussian random fields[END_REF] and use [START_REF] Didier | Exponents, symmetry groups and classification of operator fractional Brownian motions[END_REF] to compute the rejection threshold d α for fixed α.

When computing the test statistic d ≡ d H0 under the null hypothesis H 0 , the sampling distribution of d H0 is found to be independent of H 0 (see Figure 2). This result is non-trivial, but it is in line with and corroborates previous studies of the properties of wavelet-based estimators for scale-free systems (e.g., [START_REF] Veitch | A wavelet-based joint estimator of the parameters of long-range dependence[END_REF]). Test performance. Results reported in Fig. 3 (leftmost point) show that, for different settings, the proposed unimodality test [START_REF] Didier | Integral representations and properties of operator fractional Brownian motions[END_REF] has the prescribed size α = 0.05 under H 0 . To study the power of the test, we first considered instances of the alternative hypothesis where π(dH) is supported on the set {H 1 , H 2 } of two distinct Hurst exponents. Namely, we compute the proportion of rejections for the pairs of values H 1 = 0.6 -∆, H 2 = 0.6 + ∆, with ∆ ∈ [0.0, 0.10]. Figure 3 shows that the test attains power greater than 0.80 provided that H 2 -H 1 ≥ 0.125.

Next, we repeated the analysis for instances of the alternative hypothesis where π(dH) is supported on a set of 4 or 6 distinct Hurst parameters, respectively, reported in Figure 3 (Red) and (Blue). In both cases, the parameters H 1 , ..., H m are initialized at 0.6. From there, the difference between adjacent Hurst parameters is increased so as to keep the parameters equidistant and centered around 0.6. As expected, the power of the test was shown to increase as a function of the distance between the Hurst exponents in both cases. Moreover, comparative inspection of the distinct lines in Fig. 3 reveals that, starting at the same value 0.6 under the null hypothesis, the speed of convergence of the power of the test to 1 is greater when the number of distinct Hurst parameters is larger.

Overall, these results confirm that the proposed test procedure is operational, has satisfactory performance and can be readily applied in the study of self-similarity in highdimensions.

V. CONCLUSIONS AND PERSPECTIVES

A statistical methodology for detecting multimodality in the distribution of Hurst exponents in high-dimensional fractal systems has been devised. It relies on the analysis of the distribution of the log-eigenvalues of large wavelet random matrices in the threefold limit as dimension, sample size and scale go to infinity. Depending on the presence of a single or several modes in the distribution of Hurst exponents, the wavelet empirical (log)spectral distribution also displays one or several modes. A practical unimodality test is further proposed for the Hurst exponent distribution. Future work includes the construction of methodology for the identification of the distribution π(dH) once unimodality has been rejected (cf. [START_REF] Lucas | Counting the number of different scaling exponents in multivariate scale-free dynamics: Clustering by bootstrap in the wavelet domain[END_REF] in the multivariate context). Studying high-dimensional limits is not a purely mathematical issue; instead, it is absolutely crucial for practical use in applications. Indeed, multiscale analysis implies by its own nature that the effective sample size is not large compared to the number of components at coarse scales. Real data modeling also calls for the investigation of the more general case where the components of X(t) in (2) are correlated. This requires specific mathematical efforts to be made in the future.

+∞

  k=0 log k (a)H k /k!, where f.d.d.

Fig. 1 .

 1 Fig. 1. Distribution of weighted wavelet log-eigenvalues Hq in the threeway limit limn→∞ p 2 j /n. (Left) p-variate ofBm with Hurst matrix H = P diag(0.5, . . . , 0.5)P -1 ∈ R p 2 . (Right) p-variate ofBm with Hurst matrix H = P diag(H 1 , . . . , Hp)P -1 ∈ R p 2 , where H 1 , . . . , Hp are sampled from π(dH) supported on {0.2, 0.4, 0.6, 0.8}. Red curves are numerical best fits. A unimodal (left) or a four-modal (right) distribution emerges in the wavelet log-eigenspectrum, as opposed to the traditional Marcenko-Pasturbased laws in the eigenspectrum of sample covariance matrices.

Fig. 2 .

 2 Fig. 2. Quantile-quantile plots for d H 0 . The plots compare the quantiles of the distributions of d H 0 , for different pairs of H 0 : (Left) H 0 = 0.8 vs. H 0 = 0.2 ; (Right) H 0 = 0.8 vs. H 0 = 0.5 ; and show that the (non-Gaussian) distribution of d H 0 is independent of H 0 .

Fig. 3 .

 3 Fig. 3. Test powers. Proportions of rejections over 1000 Monte Carlo runs for several alternative hypothesis: (Black) supp π(dH) = {H 1 , H 2 } (Red) supp π(dH) = {H 1 , H 2 , H 3 , H 4 } (Blue) supp π(dH) = {H 1 , ..., H 6 }.


Charles-Gérard Lucas

Patrice Abry

Herwig Wendt
email: herwig.wendt@irit.fr

Gustavo Didier
email: gdidier@tulane.edu

Drowsiness detection from polysomnographic data using multivariate selfsimilarity and eigen-wavelet analysis

Because drowsiness is a major cause in vehicle accidents, its automated detection is critical. Scale-free temporal dynamics is known to be typical of physiological and body rhythms. The present work quantifies the benefits of applying a recent and original multivariate selfsimilarity analysis to several modalities of polysomnographic measurements (heart rate, blood pressure, electroencephalogram and respiration), from the MIT-BIH Polysomnographic Database, to better classify drowsiness-related sleep stages.

Clinical relevance-This study shows that probing jointly temporal dynamics amongst polysomnographic measurements, with a proposed original multivariate multiscale approach, yields a gain of above 5% in the Area-under-Curve quantifying drowsiness-related sleep stage classification performance compared to univariate analysis.

I. INTRODUCTION

Context: Physiology, selfsimilarity and drowsiness. Physiological and body temporal rhythms are well-described by arhythmic or scale-free dynamics. This is notably the case of infraslow brain activity [START_REF] He | Scale-free properties of the functional magnetic resonance imaging signal during rest and task[END_REF], [START_REF] Ciuciu | Interplay between functional connectivity and scale-free dynamics in intrinsic fMRI networks[END_REF], [START_REF] La Rocca | Self-similarity and multifractality in human brain activity: A waveletbased analysis of scale-free brain dynamics[END_REF], heart rate variability [START_REF] Yamamoto | Coarse-graining spectral analysis: new method for studying heart rate variability[END_REF], [START_REF] Ivanov | Multifractality in human heartbeat dynamics[END_REF], [START_REF] Amaral | Modeling heart rate variability by stochastic feedback[END_REF], [START_REF] Doret | Fractal analysis and hurst parameter for intrapartum fetal heart rate variability analysis: a versatile alternative to frequency bands and lf/hf ratio[END_REF], [START_REF] Nakamura | Multiscale analysis of intensive longitudinal biomedical signals and its clinical applications[END_REF], human gait variability [START_REF] Benabdelkader | Gait recognition using image self-similarity[END_REF], [START_REF] Decker | Complexity and human gait[END_REF] or sleepstage dynamics [START_REF] Leon | A bayesian framework for bivariate multifractal analysis[END_REF]. Such scale-free dynamics are often well-modeled by selfsimilarity and quantified by the corresponding so-called selfsimilarity or Hurst exponent [START_REF] Pipiras | Long-Range Dependence and Self-Similarity[END_REF]. However, physiological data very often consist of several time series recorded jointly to study one same biological mechanisms or a pathology. This is for instance the case for polysomnographic measurements, of interest here, aiming to characterize sleep stages and to quantify sleep quality, from several non-invasive modalities related to cardiovascular, respiratory or macroscopic brain activities. However, most often, selfsimilarity is assessed in individual time series, resulting in a collection of univariate analyses. This is mostly due to the fact that, up to a recent past, selfsimilarity was defined only in univariate settings and modeled by the univariate fractional Brownian motion. However, recently a proper definition of multivariate selfsimilarity was proposed [START_REF] Didier | Integral representations and properties of operator fractional Brownian motions[END_REF] together with multivariate eigen-wavelet-based analysis [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF]. The present work thus aims to show the potential of multivariate selfsimilarity in physiological and body rhythm analyses through the example of drowsiness detection from polysomnographic data.

Related work: Drowsiness detection. Drowsiness is usually defined as an intermediate state between wake and sleep [START_REF] Yu | Driver drowsiness detection using condition-adaptive representation learning framework[END_REF]. Drowsiness is documented to play a major role in vehicle accidents. Therefore, drowsiness detection from portable non-invasive biomedical measurements constitutes a significant societal stake. It has often been performed from electromyogram data, ECG data or EEG data, making use of non-linear statistical signal processing tools such as sample entropy, selfsimilarity, multifractality, cf., e.g., [START_REF] Yu | Driver drowsiness detection using condition-adaptive representation learning framework[END_REF], [START_REF] Wang | Real-time ecg-based detection of fatigue driving using sample entropy[END_REF], [START_REF] Cahyadi | Muscle fatigue detections during arm movement using emg signal[END_REF], [START_REF] Ahn | Exploring neuro-physiological correlates of drivers' mental fatigue caused by sleep deprivation using simultaneous eeg, ecg, and fnirs data[END_REF], yet leaving essentially unexplored the use of multivariate scale-free dynamics (see a contrario [START_REF] Leon | A bayesian framework for bivariate multifractal analysis[END_REF]). Goals, contributions and outline. The present work aims to quantify the benefits of using the recently devised multivariate selfsimilarity-based eigen-wavelet analysis [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF] to detect drowsiness in polysomnographic data. To that end, the estimation of a set of M selfsimilarity exponents from Mvariate data using a multivariate wavelet eigen-value based analysis is decribed in Section II. These tools are applied to M = 4-variate polysomnographic data, obtained from the MIT-BIH Polysomnographic Database, well-documented in [START_REF] Ichimaru | Development of the polysomnographic database on cd-rom[END_REF], [START_REF] Goldberger | Physiobank, physiotoolkit, and physionet: components of a new research resource for complex physiologic signals[END_REF] and described in Section III. Drowsiness detection performance are reported in Section IV.

II. MULTIVARIATE EIGEN-WAVELET ANALYSIS

Multiscale or Wavelet spectrum. The M-variate time series to be analyzed are denoted Y = {Y 1 (t), . . . ,Y M (t)} t∈R . In the next sections, Y will correspond to the collection of polysomnographic signals used to detect drowsiness.

The discrete wavelet transform (DWT) coefficients of the component Y m of Y , D Y m (2 j , k) are defined as inner products between Y m and dilated and translated templates ψ 0 (2 -j tk) of a well-chosen reference pattern ψ 0 called the mother wavelet: [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF].

D Y m (2 j , k) = 2 -j ψ 0 (2 -j t -k)|Y m (t) ,
To study scale-free dynamics, one usually forms the wavelet spectrum S(2 j ) defined, at a given scale 2 j , as a M × M matrix of wavelet coefficients intercorrelation functions [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF]:

S m,m (2 j ) 1 n j n j ∑ k=1 D Y m (2 j , k)D Y m (2 j , k) * , (1)
where n j is the number of wavelet coefficients at scale 2 j . In addition, for each pair of components (Y m ,Y m ), the wavelet coherence function C m,m (2 j ) = S m,m (2 j )/ S m,m (2 j )S m ,m (2 j ), measures a series of scaledependent correlation coefficients [START_REF] Whitcher | Wavelet analysis of covariance with application to atmospheric time series[END_REF], [START_REF] Wendt | Multivariate Hadamard self-similarity: testing fractal connectivity[END_REF]. Univariate selfsimilarity parameter estimation. The classical univariate analysis of Y consists in estimating one selfsimilarity parameter independently for each time series. In other words, it amounts to using only the diagonal entry of the matrix wavelet spectrum S(2 j ). It has been welldocumented that, for the increments of univariate selfsimilar signals such as fractional Gaussian noise, S m,m (2 j ) should (asymptotically) behave as a power law with respect to the scales 2 j , with scaling exponents driven by the selfsimilarity parameters as 2H True m -1 [START_REF] Flandrin | Wavelet analysis and synthesis of fractional Brownian motion[END_REF]. Thus, parameter estimation can be efficiently performed by means of linear regressions of log 2 (S m,m (2 j )) against log 2 (2 j ) = j [START_REF] Veitch | A wavelet-based joint estimator of the parameters of long-range dependence[END_REF]:

ĤU m = j 2 ∑ j= j 1 v j log 2 (S m,m (2 j )) 2- 1 2 , m = 1, . . . , M, (2)
with v j such that ∑ j v j = 1 and ∑ j v j = 0. Classical multivariate selfsimilarity parameter estimation. To account for cross-dependencies amongst component temporal dynamics, a natural approach is to apply the same procedure to each entry of the wavelet spectrum S m,m (2 j ), thus defining a set of cross-selfsimilarity exponents:

Ĥm,m = j 2 ∑ j= j 1 v j log 2 (S m,m (2 j )) 2- 1 2 , m, m = 1, . . . , M.
(3) Eigen wavelet selfsimilarity parameter estimation. However, it has been recently documented that multivariate selfsimilarity can be better conducted based on a change of paradigm: Instead of considering for each pair of components, m and m , the behavior of S m,m (2 j ) along scales and then comparing between pairs of components the estimated selfsimilarity parameter Ĥm,m , as above, one can first consider at a given scale 2 j all components together, by computing the eigenvalues λ m (2 j ) of the spectrum matrix S(2 j ), and then study the behavior of these eigen-values along scales 2 j . It has indeed been shown that for multivariate selfsimilar processes the eigenvalues λ m (2 j ) should (asymptotically) behave as a power law with respect to the scales 2 j , with scaling exponents related to the selfsimilarity parameters as 2H True m -1 [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF]. However, the practical computation of the eigenvalues of the set of matrices S(2 j ) from finite sample size multivariate data suffers from a bias, referred to as the repulsion bias of significant impact at coarse scales, where the number of available wavelet coefficients becomes limited compared to the number of components M. This recently lead to the construction of a robust estimation based on the eigenvalues, from a collection of estimated wavelet spectra obtained from non-overlapping windows of width w, using an equal number n j 2 of wavelet coefficients D Y (2 j , k) at each scale 2 j (with n j 2 the number of wavelet coefficients at coarsest scale 2 j 2 ) [START_REF] Lucas | Bootstrap for testing the equality of selfsimilarity exponents across multivariate time series[END_REF]:

S (w) m,m (2 j ) 1 n j 2 wn j 2 ∑ k=1+(w-1)n j 2 D Y m (2 j , k)D Y m (2 j , k) * . (4)
The log-eigenvalues λ (w)

1 (2 j ), . . . , λ (w)
M (2 j ) of the matrices S (w) (2 j ) are then averaged to yield λm (2 j )

2 j 2 -j ∑ 2 j-j 2 w=1 log 2 (λ (w)
m (2 j )). Eigenwavelet-based estimates ĤM m of the multivariate selfsimilarity parameters are finally obtained by linear regressions of the log-averaged eigenvalues λm (2 j ) against log 2 2 j = j:

ĤM m = j 2 ∑ j= j 1 v j λm (2 j ) 2 - 1 2 , m = 1, . . . , M. (5)
III. DATABASIS MIT-BIH Polysomnographic database. Data used here are those made available in the MIT-BIH Polysomnographic Database 1 , documented in [START_REF] Ichimaru | Development of the polysomnographic database on cd-rom[END_REF], [START_REF] Goldberger | Physiobank, physiotoolkit, and physionet: components of a new research resource for complex physiologic signals[END_REF]. Data consist of a collection of 4 to 7 physiological measurements (cardiovascular, respiratory or macroscopic brain activities, stroke volume, oxygen saturation, eye movements or chin muscle responses) recorded for 16 male subjects in Boston's Beth Israel Hospital Sleep Laboratory in the context of chronic obstructive sleep apnea syndrome; the recordings of the first two subjects are split into two consecutive portions, resulting in a total of 18 multivariate time series (sampling 250Hz, data lengths from 77 to 390 min). Sleep stage expert annotations are available for each 30-second long non overlapping window. Data for the present study. Because the focus of this work is on detecting drowsiness, defined as transitions between the awake and stage1 stages, the goal here is to perform a classification of these two stages. Use is thus made only of the time windows corresponding to such annotations.

Further, to investigate the benefits of a joint analysis based on multivariate selfsimilarity, one needs to analyze the different modalities of polysomnographic data in the same frequency range. Because respiratory and heart rate rhythms contain relevant information in time scales ranging from 1 second to 1 minute, only the infraslow brain activity is considered and data are filtered and resampled at 4Hz.

Finally, to ensure a consistent number of features for detection and classification, only those 4 modalities of the polysomnographic data are used that are available for all subjects: heart rate (HR), blood pressure (BP), electroencephalogram (EEG) and respiration (RESP).

IV. SLEEP STAGE CLASSIFICATION

A. Analysis/classification set-up 1) Physiological time series analysis: Analysis is performed in 2-min long sliding windows, with 75% of overlap between successive windows, each hence containing N w = 480 = 4 × 2 × 60 samples.

Classification will be restricted to sequences sharing the same annotations for at least 4 consecutive 30-second windows. In total, 1753 and 561 such windows are available for the awake and stage 1 states, respectively.

The DWT coefficients D Y m (2 j , k) are computed using the least asymmetric Daubechies3 wavelet [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. Linear regressions are performed across scales ranging from 2 j 1 = 2 1 to 2 j 2 = 2 4 corresponding to frequencies from 1/8 to 2 Hz, or equivalently to time scales from 1/2 to 8 seconds.

Figs. 1 and 2 plot, as an example, log-wavelet spectrum and log-eigenvalues, for one 2-min long window, indicating power law scaling thus selfsimilarity across the chosen analysis scales. Scaling exponents are estimated using Eqs. 3, 2 and 5 from such functions.

2) Multifeature classification: We make use of a standard Random Forest Classifier [START_REF] Breiman | Random forests[END_REF], a learning method consisting in training a large number of decision trees from data resampled with replacement. Each tree is trained from a subset of features randomly selected among the N f available features, to reduce correlation between trees. The size of the subset list is chosen here as N f as in [START_REF] Breiman | Random forests[END_REF], and decision is made by majority vote.

Random Forest Classifiers are performed with N trees ∈ {10, 25, 50} trees. A diagonal cost matrix is set with coefficients w 1 = W v N a /N W and w 2 = N s /N W , where N a is the number of windows related to the state "awake" (class 0), N s is the number of windows related to the state "stage 1" (class 1), N W = N a + N s is the total number of windows and W v ∈ [0.001, 6] to tune the classification false alarm rate. Performances are assessed by cross-validation over N MC = 100 repetitions, with 80% of the available windows randomly and independently selected for the training sets.

Receiver Operational Characteristic (ROC) curves are computed by varying W v during training. ROC areas-undercurve (AUC) are used as performance score.

B. Single-feature classification

As a baseline, classification is first performed for each of the four modalities independently, using as a single feature the univariate selfsimilarity parameter, ĤU m , estimated using Eq. 2. Classification hence simply consists in comparing the feature against a threshold, and ROC curves are computed by varying the classification threshold as shown in Fig. 3. AUC, reported in Table I

C. Multi-feature classification

The goal is now to quantify the benefits of using several polysomnographic modalities jointly in classifying drowsiness related sleep stages.

1) Features: Three different strategies to combine information stemming from the different modalities are tested.

The simplest way consists in concatenating the four univariate ĤU m into a feature vector of dimension 4. This will be referred to as features1. Such a joint modality classification does not account for information related to cross-temporal dynamics between modalities as it only makes use of features computed independently on each modality. However, as illustrated in Fig. 1, there exist non-negligible cross-temporal dynamics at all scales between modalities, thus prompting for the use of such information to improve classification.

To account for cross-temporal dynamics across modalities, a second classification will be performed by adding to the 4 univariate features ĤU m , the 6 = 4 × 3/2 cross-selfsimilarity parameters Ĥm,m (m = m) estimated from the nondiagonal entries of the wavelet spectrum (3), resulting in a feature vector of dimension 10 (referred to as features2).

In addition, an original contribution of this work is to promote the use of the eigenwavelet analysis as a new way to quantify cross-temporal dynamics. Therefore, a third classification will be performed by adding to the 4 univariate features ĤU m , the 4 eigenvalue based multivariate selfsimilarity parameters ĤM m estimated as in [START_REF] Ivanov | Multifractality in human heartbeat dynamics[END_REF], resulting in a feature vector of dimension 8. This will be referred to as features3.

2) Performance: To compare performance between these three different multimodal classification strategies, Fig. 4(right plots) displays ROC curves and Table II further reports the AUC. These results yield the following conclusions. Multi-feature classifications significantly outperform single-feature classification. For multi-feature classifications, the performance is robust to the choice of the number of trees in the random forest procedure. Amongst the multi-feature classifications, the one based on the classical measures of cross-temporal dynamics (features2) does not improve classification performance compared to the simplest concatenation of the univariate features (features1). To the contrary, the promoted eigen-wavelet approach to probing cross-temporal dynamics, combined with the univariate features (features3), improves significantly drowsiness related sleep-stage classification performance.

V. CONCLUSIONS

The present work has first quantified the benefits of sleep multimodal monitoring in phase sleep classification. Second, it has shown that compared to simply combining by concatenation univariate features, the proposed multivariate eigenwavelet approach strengthens sleep stage classification. It clearly quantifies that cross-temporal dynamics amongst modalities do carry information relevant to the sleep state evaluation. This is because the smallest eigenvalues are able to sense to low-intensity yet meaningful M-wise dependencies amongst components that classical pairwise crosscorrelations may miss.

1Fig. 1 :Fig. 2 :

 12 Fig. 1: Log-wavelet spectrum for one subject and one 2-min long window.

Fig. 3 :

 3 Fig. 3: ROC curves (average ± standard deviation) for single-feature classification.

Fig. 4 :

 4 Fig. 4: ROC curves (average ± standard deviation), for the multivariate classifications with three different numbers of trees.

  , show that Blood Pressure yields the best classification performance. However, single-feature performance remains quite low.

	feature ĤU 1	feature ĤU 2	feature ĤU 3	feature ĤU 4
	52.71	67.59	55.15	57.44

TABLE I :

 I AUC for univariate classifications.

TABLE II :

 II AUC (average ± standard deviation) for the multivariate classification.

Work supported by PhD Grant DGA/AID (no 01D20019023)