Gilles Blanchard

Jean-Baptiste Fermanian

Nonasymptotic one-and two-sample tests in high dimension with unknown covariance structure

Keywords: Signal detection, Two-sample test, Relevant hypotheses, Minmax testing separation distance, Effective dimensionality

.

Introduction

We consider the following fundamental signal detection problem: given an i.i.d. sample X = (X i ) 1≤i≤n from a square integrable distribution P X on R d (or possibly a separable Hilbert space, under some conditions which will be discussed later) with µ = E[X 1 ], test the hypothesis of "η-closeness to zero" of the mean:

(H 0 (η)) : µ ≤ η, against (H 1 (η, δ)) : µ > η + δ.

(1)

In fact, we consider the following more general two-sample mean closeness testing problem: for X = (X i ) 1≤i≤n and Y = (Y i ) 1≤i≤m two independent samples of i.i.d.

variables with distributions P X , P Y on R d with respective means µ and ν, test the hypothesis of η-closeness (or similarity) of the two means, (H 0 (η)) : µν ≤ η, against (H 1 (η, δ)) : µν > η + δ.

(2)

Observe that we can always formally subsume setting (1) into setting (2), by letting m go to infinity and/or assuming (if needed) that the covariance of Y 1 is zero. Therefore, in the contribution section we will concentrate mainly on setting (2).

The problem (1) (and numerous extensions thereof) has been a long-time subject of attention in mathematical statistics. For the zero mean test problem, i.e. η = 0, the celebrated works of [START_REF] Ingster | Minimax nonparametric detection of signals in white Gaussian noise[END_REF][START_REF] Ingster | Asymptotically minimax hypothesis testing for nonparametric alternatives I-II-III[END_REF] in the Gaussian white noise model are seminal. In the case η > 0, the problems (1)-( 2) are known as testing for precise hypotheses, relevant hypotheses or relevant differences [START_REF] Berger | Testing Precise Hypotheses[END_REF]; this setting has found applications in particular in biostatistics for bioequivalence testing (see e.g. [START_REF] Wellek | Testing statistical hypotheses of equivalence[END_REF]. (See next sections for a more detailed discussion of related literature.) In this work, we will consider the situation where the involved distributions are either Gaussian or of bounded norm (and hence sub-Gaussian), but with unknown covariance matrix acting as a nuisance parameter.

We are interested in finding bounds on the separation distance δ, i.e. a bound on the minimum value of δ such that there exists a test with both Type I and Type II error rates bounded by a "small" prescribed quantity. Our interest here is more on the constructive side, so that we will concentrate on feasible procedures that are in particular adaptive to the covariances of the involved distributions. A matching lower bound (for any fixed covariance structure) will be provided in the Gaussian setting. We emphasize that our focus is on finite sample (i.e. nonasymptotic) results, as will be discussed below.

Relation to white noise model in nonparametric statistics

In the isotropic Gaussian case (white noise) with known variance , and for η = 0, the signal detection problem (1) has been studied in much generality, in particular in the infinite-dimensional setting where R d is replaced by a separable Hilbert space. In this situation, due to the fact that the white noise model on an infinite-dimensional Hilbert space cannot be represented by a random variable taking values in that space, the canonical model which is considered instead is the Gaussian sequence model for the coordinates of each of the observations in an orthonormal basis (in fact the Gaussian sequence model with known variance is usually considered with a single observation of the sequence):

X (i) = µ (i) + σε (i) , i ∈ N >0 , (3)
where (ε (i) ) i≥1 is an i.i.d. standard normal sequence. This fundamental model in nonparametric statistics allows to represent in a clean way many functional spaces of interest for the signal µ through geometrical properties of its expansion coefficients (µ (i) ) i≥1 in a suitable basis. Since in that infinite-dimensional setting the alternative µ 2 > δ 2 is "too big" and gives rise to trivial separation rates, the usual focus is on considering restricted alternatives of the form µ ∈ F ; µ ≥ δ 2 , for a given nonparametric set F . Classical alternatives of interest include in particular ℓ 2 ellipsoids (corresponding to Hilbert norms of different strengths), ℓ p bodies, and Besov bodies. Interpreted in functional spaces, these alternatives correspond respectively to balls in Sobolev spaces (typically when considering Fourier basis coefficient expansions) or in Besov spaces (for suitable wavelet basis coefficient expansions).

The literature on these topics is profound and extensive, see e.g. [START_REF] Ingster | Nonparametric goodness-of-fit testing under Gaussian models[END_REF] for a comprehensive overview. The case of certain classes of ℓ 2 -ellipsoids appears to have been studied first by [START_REF] Ingster | Minimax nonparametric detection of signals in white Gaussian noise[END_REF] and [START_REF] Ermakov | Minimax detection of a signal in a Gaussian white noise[END_REF], then a remarkable series of works of [START_REF] Ingster | Asymptotically minimax hypothesis testing for nonparametric alternatives I-II-III[END_REF] and [START_REF] Ingster | Minimax detection of a signal for Besov bodies and balls[END_REF] established minimax testing rates for general ℓ 2 ellipsoids as well as other alternatives. V. Spokoiny's contribution is prominent in this body of literature, in particular for dealing with the case of Besov bodies [START_REF] Lepski | Minimax nonparametric hypothesis testing: the case of an inhomogeneous alternative[END_REF] as well as considering the problem of statistical adaptivity over a family of alternatives [START_REF] Spokoiny | Adaptive hypothesis testing using wavelets[END_REF].

This very limited overview of the topic of testing in the white noise model is meant to contrast with the setting considered here. On the one hand, we will not consider a particular form of alternative; on the other hand, we assume that the observations can truly be represented as elements in a possibly infinite-dimensional separable Hilbert space. Under the Gaussian assumption, this means that the covariance operator Σ of the noise process is assumed to have a finite trace, which also prevents the triviality problem mentioned above for the white Gaussian noise setting. If we represent the observation coordinates in a diagonalizing basis of Σ, our setting in the Gaussian setting amounts to the Gaussian sequence model (3) wherein the constant parameter σ is replaced by a square integrable sequence (σ (i) ) i≥1 . Note that formally normalizing the i-th observation coordinate by σ (i) would give rise again to model (3), however the separation distance would then be measured in the weak norm Σ 1 /2 µ .

Relation to "modern" and high-dimensional statistics

Since we only consider test separation distance without a specific alternative, the setup we consider can be considered as less elaborate, at least in the sense of asymptotic theory, than the settings with various non-parametric alternatives discussed above. On the other hand, our focus is specifically on the following points:

1. Finite-sample analysis; 2. Non-Gaussian data (we will only consider bounded data here); 3. Robustness to misspecification (here under the form of the relaxed composite null µ 2 ≤ η 2 , also called relevant hypothesis testing).

These features have been rightly identified by V. Spokoiny as the defining features of "modern" approach to statistics [START_REF] Spokoiny | Parametric estimation. Finite sample theory[END_REF][START_REF] Spokoiny | Basics of modern mathematical statistics[END_REF]. The problem of testing a null hypothesis defined as a neighborhood rather than an exact match has been tackled under different settings in the statistics literature, especially for the two-sample testing case. For example, motivated by bioequivalence testing between populations, [START_REF] Munk | Nonparametric validation of similar distributions and assessment of goodness of fit[END_REF] consider the problem of testing closeness of two real distributions as measured in Mallows distance, [START_REF] Dette | Nonparametric comparison of several regression functions: exact and asymptotic theory[END_REF] that of closeness in L 2 distance of two nonparametric (Hölder regular) regression functions; Dette, Kokot, and Aue (2020), the closeness in supremum norm distance of two mean functions in a Banach functional data setting; [START_REF] Dette | Testing relevant hypotheses in functional time series via self-normalization[END_REF], the closeness in L 2 distance of the functional mean of time series. In all cases, the underlying principle is to estimate the target distance -as will be also case in the present paper -and the data is not always assumed to be Gaussian, but the corresponding analysis based on Gaussian asymptotic theory. To estimate the quantiles of the test statistic, [START_REF] Dette | Nonparametric comparison of several regression functions: exact and asymptotic theory[END_REF] choose to estimate the variance, [START_REF] Dette | Testing relevant hypotheses in functional time series via self-normalization[END_REF] use a self-normalized procedure and give asymptotic bounds; [START_REF] Dette | Functional data analysis in the Banach space of continuous functions[END_REF] propose a bootstrap approach and obtain an asymptotic convergence of the test statistic. In the present paper our approach is a direct estimation of the variance with nonasymptotic guarantees.

Taking the above aspects into account in the theory, in particular nonasymptotic analysis, is motivated by a large number of high-dimensional applications, where it appears that relying on traditional asymptotic of Gaussian parametric or non-parametric theory can possibly be problematic if done without care. Finite sample theory allows to delineate more precisely in which situations traditional approximations still can be relied upon, and to study non-standard asymptotics, in particular when key parameters, such as dimensionality, can themselves depend on the sample size n. It is also of use when considering multiple testing scenarios, where multiplicity has to be taken into account precisely.

Another fruitful modern insight is that high-dimensional statistical models tend to blur the line between parametric and non-parametric point of views. Precise non-asymptotic results in a finite-dimensional setting, but where the role of key model parameters (in particular, dimensionality or effective dimensionality) is precisely analyzed, can provide key theoretical components for analyzing nonparametric settings. In the signal testing framework considered in the present paper, this way of thinking has in particular been pioneered by [START_REF] Baraud | Non-asymptotic minimax rates of testing in signal detection[END_REF], who obtained sharp non-asymptotic results for the problem (1) in the case η = 0, and for the finite-dimensional counterpart of the white noise model (3), i.e. the isotropic setting Σ = σ 2 I d in dimension d. Baraud further demonstrated that this result provided a valuable and versatile tool to analyze models of typical interest in high-dimensional statistics (such as sparse alternatives) as well as non-parametric alternatives (such as those mentioned in the previous section).

A key insight from Baraud's work is that the minimum separation distance in that setting is O(d 1 /4 σ/ √ n), in contrast with minimax estimation distance for µ, which is Θ(d 1 /2 σ/ √ n): the testing separation distance is smaller than the minimax estimation error by a factor d 1 /4 . Analyzing precisely the role of dimensionality (ambient or effective) in minimax testing separation rates and the difference with minimax estimation rates has been a subject of interest in recent literature in various settings, highlighting similar related phenomena. For instance, [START_REF] Lam-Weil | Local minimax rates for closeness testing of discrete distributions[END_REF] consider the problem of testing equality of two high-dimensional multinomial distributions and study the minimum ℓ 1 separation distance in a vicinity of a reference distribution π (which implicitly determines a notion of local effective dimensionality). Since this model has bounded data, our analysis could be applied in that setting, however it concerns separation in ℓ 2 distance (the separation in ℓ 1 distance exhibits considerably more involved behavior). [START_REF] Ostrovskii | Near-Optimal Model Discrimination with Non-Disclosure[END_REF] consider a different type of two-sample testing problem, in a regression context, where the goal is to determine which one of the two distributions has a given (known to the user) regression vector. They give a sharp bound on the minimum separation distance between the two regression vectors including the role of the dimension, also exhibiting a difference with estimation rates.

Coming back to our model, the results of [START_REF] Baraud | Non-asymptotic minimax rates of testing in signal detection[END_REF] provide a sharp answer, but only in the case η = 0 and for isotropic Gaussian (white noise) data with known variance. Still in the Gaussian isotropic case, the minimum separation rates for any value of η ≥ 0 were precisely characterized by [START_REF] Blanchard | Minimax Euclidean separation rates for testing convex hypotheses in R d[END_REF]. We also consider the Gaussian setting in the present work, but analyze the generalized situation where the covariance matrix Σ can be arbitrary (and unknown). In this situation, the role of the dimensionality d is played by proxy quantities depending on Σ, sometimes called effective dimensionality or effective rank. For the signal testing problem however, it turns out that the proxy dimensionalities for testing and estimation differ. Namely, for η = 0, we find that the minimax separation distance is O(d

1 /4 * Σ ∞ /n), where d * := Σ 2 2 / Σ 2 ∞ , while the minimax estimation distance for µ is Θ(d 1 /2 e Σ ∞ /n),
where d e := Σ 1 / Σ ∞ . (Notice that d * ≤ d e ≤ d in general, while these quantities are all equal in the isotropic setting.) Furthermore, we also study the estimation of key quantities Σ 1 /2 ∞ and Σ 2 determining the proxy dimensionality and the testing threshold 4 .

A crucial mathematical tool in high-dimensional statistics is to obtain sharp concentration inequalities for quadratic forms of random vectors. These are closely related to technical tools used in the present work. An important point in such inequalities is to quantify as precisely as possible up to which point quadratic forms of non-Gaussian vectors can mimic the Gaussian behavior (i.e. that of central and non-central weighted chi-squared statistics). This topic has received a good deal of attention in the recent years and V. Spokoiny also made substantial contributions to that area [START_REF] Spokoiny | Basics of modern mathematical statistics[END_REF][START_REF] Spokoiny | Sharp deviation bounds for quadratic forms[END_REF]. In the present work, we derive from scratch the needed concentration inequalities; we discuss in more detail the relation to V. Spokoiny's own work and to related literature in Section 2.4.

Relation to machine learning and kernel mean embeddings of distributions

An application setting which motivated us to consider in detail the case of bounded data is that of testing of the data distribution via kernel mean embedding (KME) methods, a principle which has garnered a lot of attention in the machine literature since the seminal paper of [START_REF] Smola | A Hilbert space embedding for distributions[END_REF]. It has been advocated in particular for two-sample [START_REF] Gretton | A kernel two-sample test[END_REF] and goodness-offit [START_REF] Chwialkowski | A kernel test of goodness of fit[END_REF] testing; see [START_REF] Muandet | Kernel mean embedding of distributions: a review and beyond[END_REF] for a recent overview.

We describe the KME principle briefly. Assume Z is a random variable with distribution P Z taking values in the measurable space Z, and that one has at hand a fixed mapping Φ : Z → H, where H is a separable Hilbert space. To this mapping is associated a reproducing kernel Hilbert space (rkHs) H ′ with kernel k(z, z ′ ) := Φ(z), Φ(z ′ ) .

Assuming the variable X = Φ(Z) is Bochner integrable5 (which is the case in particular when the mapping Φ is bounded), the kernel mean embedding of P Z is defined as Φ(P Z ) := E[Φ(Z)] ∈ H (using a rather natural overload of notation for Φ). The maximum mean discrepancy (MMD) between distributions P, Q in the domain of definition of Φ is defined as the semimetric

MMD k (P, Q) := Φ(P) -Φ(Q) .
Since MMD k (P, Q) > 0 implies P = Q, this principle can be used for simple goodness-of-fit testing (testing for P Z = P 0 for some known distribution P 0 , given an i.i.d. sample from P Z ) and two-sample testing (testing for P Z = P Z ′ , given two independent i.i.d. samples from P Z and P Z ′ ); in each case, the test statistic is a suitable estimator of MMD k (P Z , P 0 ), resp. MMD k (P Z , P Z ′ ) from the observed data. More generally one may want to test the relaxed null hypothesis MMD k (P, Q) ≤ η and analyze the power of the test in terms of the MMD separation itself. This is indeed a particular case of (1)-(2), when considering the Hilbert-valued variable X = Φ(Z) and, for two-sample testing, Y = Φ(Z ′ ).

A common situation is when Φ is bounded in norm by some constant L, or equivalently in terms of the kernel, sup z∈Z k(z, z) ≤ L 2 . This ensures in particular that Φ is defined on all distributions. Analyzing our original setting with norm-bounded but potentially infinite-dimensional data is therefore suited to this case. [START_REF] Gretton | A kernel two-sample test[END_REF] derive the asymptotic distribution of the (suitably renormalized) MMD test statistic, which is identical to the one we use below (once interpreted in the KME setting). Unsurprisingly, a Gaussian limiting behavior is identified. Our study analyzes this behavior from a non-asymptotic point of view; this can be particularly of interest for situation where the mapping Φ (or equivalently the associated kernel) is to depend on the sample size, or when performing a large number of such tests in parallel: in this case uniformly valid nonasymptotic bounds are a a valuable tool for further analysis. See [START_REF] Marienwald | High-dimensional multi-task averaging and application to kernel mean embedding[END_REF] for such a multiple test scenario in the context of so-called multiple task averaging. Multiple tests can also be aggregated to test a global hypothesis, see [START_REF] Fromont | Kernels based tests with non-asymptotic bootstrap approaches for two-sample problems[END_REF] in the context of two-sample testing based on the KME approach.

In our study, the power of the test is investigated for alternatives of the form (1)-( 2), which, interpreted in the KME setting, correspond to MMD k (P, Q) ≥ η + δ. The power of KME-based tests (in the goodness-of-fit case) was also investigated by [START_REF] Balasubramanian | On the optimality of kernelembedding based goodness-of-fit tests[END_REF], but for alternatives measured in a χ 2 distance separation, more precisely, of the form Q ∈ F ; χ 2 (P 0 , Q) ≥ δ , where F is a nonparametric set of distributions whose density with respect to P 0 is approximated at a given rate by functions in the rkHs H ′ associated to k, in the sense of interpolation with L 2 (P 0 ). This is close in spirit to nonparametric points of view discussed in Section 1.1, in the sense that χ 2 -separation alone is too weak to get nontrivial separation rates and one has to additionally consider intersection with nonparametric sets of interest. Again, because we choose to analyze alternatives measured in MMD k -separation itself, the results we obtain in this setting have a different nature.

Overview of contributions

The main contribution of this paper is to give upper bounds on the optimal (minmax) testing separation distance for problems (1) and (2) over classes of probability distributions with fixed covariance matrix Σ for sample X, as well as S for sample Y in the two-sample case. The covariance structures are considered as nuisance parameters and we investigate precisely how they influence the testing separation distance. Let P be a family of distributions for the two samples (we consider the Gaussian setting and the bounded setting), and P Σ,S the subsets of distributions of P with Cov[X 1 ] = Σ, and Cov[Y 1 ] = S (in the two-sample case). Consider the sets of distributions

H 0 (η, Σ, S) := {P ∈ P Σ,S |P satisfies H 0 (η)} , A δ (η, Σ, S) := {P ∈ P Σ,S |P satisfies H 1 (η, δ)} ,
then the optimal separation distance is, for α ∈ (0, 1):

δ * (α, Σ, S, η) = inf δ ≥ 0 ∃ test T : sup P∈H0 P(T = 1) + sup P∈A δ P(T = 0) ≤ α .
(4) In the Gaussian setting, we establish that δ * is upper bounded up to a constant factor via

δ * (α, Σ, S, η) σκ α max 1, min d 1 4 * , d 1 2 * σκ α η , (5)
(Theorem 7) where κ α := log(α), and, in the one-sample case, σ 2 := Σ op /n is a scalar variance factor and d * := Tr Σ 2 / Σ 2 op a notion of effective dimension. In the two-sample case, we obtain also (5), with σ 2 := M m,n op , and

d * := Tr M 2 m,n /σ 4
, where M m,n := (Σ/n+S/m) (Theorem 8). In the one-sample case, this result can be formulated equivalently in terms of sample complexity n * needed to detect at given error level α and separation distance δ for problem (1):

n * (α, Σ, S, η) Σ op κ α δ -1 max δ -1 , d 1 2 * (max(δ, η)) -1 . ( 6
)
This result is established first when assuming that Σ, S are known, then we show that it holds as well when they are unknown (under some mild assumptions on the sample size, see Corollary 14 for an explicit statement in the one-sample case and condition (27) there). Matching minimax lower bounds are given for one and two-sample problems in the Gaussian setting. In the bounded setting, we derive upper bounds only, which take the same flavor as ( 5) under some mild assumptions on the sample sizes.

Organization of the paper

We present in Section 2 our main results. In order to cover both the Gaussian and bounded settings under the same umbrella, we start in Section 2.1 by a generic result: assuming some suitable concentration for an estimate U of the squared signal norm µ 2 holds (Assumption 1), as well as for estimators of its quantiles (Assumption 2), for the problems ( 1) and ( 2) we propose in Theorem 3 sufficient conditions on δ such that we can control the Type I and Type II errors of a test T based on U . In the following sections, the Gaussian setting and the bounded setting are considered separately. In Section 2.2, we give concentration results for U to fulfill Assumption 1. In Section 2.3 we give results to fulfill Assumption 2, which are related to the estimation of Σ 1 /2 ∞ and Σ 2 . The proofs of the corresponding results are found in Sections 3.1 to 3.6, respectively.

Main results

We will build a test for the model (2) based on an estimator U of the distance µν 2 , typically a modified U-statistic as defined below. We will first consider a general point of view to deduce bounds on the separation rate when U satisfies certain concentration properties; this will then apply both to the Gaussian and bounded settings.

A general result to upper bound separation rates

As mentioned earlier, from now on we concentrate primarily on the two-sample setting, being understood that upper bounds for the one-sample setting can be deduced readily. In order to define a general framework encompassing as particular cases the more specific settings considered below, in this section we will assume a generic statistical model P for the distribution of the samples X and Y, which we recall we always assume to be independent and i.i.d. with respective squared integrable marginal distributions P X , P Y . We will thus use without comment the fact that a distribution P ∈ P equivalently specifies the marginal distributions P X and P Y of the samples. We will consider the covariance matrices Σ, S of P X , P Y as nuisance parameters influencing the optimal separation distance, and define the sub-models

P Σ,S = {P ∈ P : Cov[P X ] = Σ, Cov[P Y ] = S};
P Σ is defined in an analogous way for the one-sample setting.

The first property we require is a form of 2-sided concentration of U around the target quantity:

Assumption 1. For any (Σ, S) and distribution P ∈ P Σ,S ; for any given α ∈ (0, 1) there exist q 1 = q 1 (Σ, S, α), q 2 = q 2 (Σ, S, α) in R + such that:

P U -µ -ν 2 ≥ µ -ν q 1 + q 2 ≤ α . (7)
Additionally, we will consider the situation where the quantities q 1 , q 2 (which are necessary to find a suitable testing threshold) are not known but must also be estimated from the data; this is the case if the covariance matrices (Σ, S) are unknown. This leads us to our second assumption:

Assumption 2. Suppose Assumption 1 holds, with the notation introduced therein. For any α ∈ (0, 1) there exist two estimators

Q 1 = Q 1 (α) and Q 2 = Q 2 (α)
in R + such that, for any (Σ, S) and distribution P ∈ P Σ,S :

P q 1 (Σ, S, α) -Q 1 (α) ≥ 1 2 q 1 (Σ, S, α) ≤ α , (8)
P q 2 (Σ, S, α) -Q 2 (α) ≥ 1 2 q 2 (Σ, S, α) ≤ α . (9)
(In the "oracle" case where the covariances Σ, S are assumed to be known, of course Assumption 2 is trivially satisfied taking

Q 1 = q 1 , Q 2 = q 2 .)
The following generic result transforms the above assumptions into an estimate of the separation distance for setting (2).

Theorem 3. Let P be a statistical model for setting (2), and U be a statistic.

Let Assumptions 1 and 2 be granted. Given η ≥ 0 and α ∈ (0, 1), let T be the test defined by

T = 1 U -η 2 > 2η Q 1 (α) + 2 Q 2 (α) . ( 10
)
Then for any (Σ, S), provided

δ > 2q 1 + min 2 √ q 2 , 2η -1 q 2 , ( 11
)
it holds, for any distribution P ∈ P Σ,S :

P[T = 1] ≤ 3α , if P satisfies (H 0 (η)); P[T = 0] ≤ 3α , if P satisfies (H 1 (η, δ)).

Concentration properties of the test statistic

The rest of the paper is dedicated to establishing the validity of Assumptions 1 and 2 for the following statistic U (X, Y):

U (X, Y) := 1 n(n -1) n i,j=1 i =j X i , X j + 1 m(m -1) m i,j=1 i =j Y i , Y j - 2 nm n i=1 m j=1 X i , Y j .
(12) Observe that provided expectations µ, ν exist, U (X, Y) is an unbiased estimator of µν 2 . In the KME setting as described in Section 1.3, inner products are replaced by kernel evaluations and the above statistic is the standard unbiased estimate of the squared MMD between P X and P Y . As announced previously, we will concentrate on the following two settings: Definition 4 (Gaussian setting). The samples X and Y are i.i.d. Gaussian in R d of marginal distributions P X = N (µ, Σ) and P Y = N (ν, S), respectively.

In the Gaussian setting, we will assume a finite ambient dimension d for technical reasons: our proofs rely on the Gauss-Lipschitz concentration inequality, which applies in finite dimension. As will appear clearly however, all our results to come are dimension-free in the sense that d never enters the picture, instead only norms of Σ, S come into play. We surmise that our results would apply as well in the same form in the Hilbert-valued setting provided Tr(Σ) and Tr(S) are finite, but did not try to write down a precise approximation argument to this end.

Definition 5 (Bounded setting). The samples X and Y are i.i.d. in a separable Hilbert space H with norm bounded by L > 0. The covariance operators for the marginal sample distributions are denoted Σ and S, respectively; observe that they have finite trace by the boundedness assumption.

Propositions 6 and 9 give concentration bounds for the statistic U , ensuring Assumption 1 in the two above settings. Proposition 6. Assume the Gaussian setting holds and n, m ≥ 2. Then with probability at least 1α,

U -µ -ν 2 ≤ µ -ν q 1 + q 2 , ( 13
)
where U is defined in (12) and

q 1 (Σ, S, α) = 2 Σ op n + S op m u(α) , (14)
q 2 (Σ, S, α) = 32 √ Tr Σ 2 n + √ Tr S 2 m u(α) . ( 15
)
where u(α) :=log α + log 8.

Let us simplify somewhat the above expression when plugged into Theorem 3 in the case of signal detection (1). We also give a matching lower bound (up to constant factor) for the optimal separation distance.

Theorem 7. Consider the signal detection problem (1) and assume the Gaussian setting with covariance matrix Σ. Then the minimum separation distance δ * given by (4) so that the type I and II errors for problem (1) are less that α ∈ (0, 1) for all P ∈ P Σ is upper bounded by

δ * (α, Σ, η) σ n √ u max 1, min d 1 4 * , d * u • σ n η , (16)
where u(α) :=log α + log 60. If d * ≥ 3, then it is lower bounded by Observe that it holds d * ≤ d e , where d e = Tr Σ/σ 2 is the "effective dimensionality" coming into play for signal estimation rates (namely E Xµ 2 1 /2 = σ d e /n, where X is the empirical mean). In the finite d-dimensional case with Σ = I d , it holds d = d e = d * , and the separation ( 16) has been shown to be optimal in the Gaussian setting for η = 0 by [START_REF] Baraud | Non-asymptotic minimax rates of testing in signal detection[END_REF] and for any η ≥ 0 by [START_REF] Blanchard | Minimax Euclidean separation rates for testing convex hypotheses in R d[END_REF]. It exhibits a continuous transition between the signal detection setting (η = 0, δ * ≃ d 1 /4 σ/ √ n) and the hyperplane testing setting (which is equivalent to the 1-dimensional setting by rotational invariance;

δ * (α, Σ, η) ≥ σ n 1 -α 12 max 1, min d 1 4 * , d * (1 -α) • σ n η , (17
η → ∞, δ * ≃ σ/ √ n).
In that particular situation, we observe that the signal separation distance is smaller by a factor d 1 /4 than the signal estimation error, a phenomenon typical of high-dimensional statistics. In the more general setting studied here where Σ can be arbitrary, this difference between rates can be all the more marked, since in addition d * can be much smaller than d e .

We obtain a similar result for the two-sample problem:

Theorem 8. Consider the two-sample mean problem (2) and assume the Gaussian setting with covariance matrices Σ, S. Then the minimum separation distance δ * so that the type I and II errors for problem (2) is less than α ∈ (0, 1) for all P ∈ P Σ,S is upper bounded by

δ * (α, Σ, S, η) σ n,m √ u max 1, min d 1 4 * , d * u • σ n,m η , ( 18
)
where u :=log α + log 60. If d * ≥ 3, then it is lower bounded by Here the effective dimension d * depends on the two covariance matrices Σ and S, weighted by the size of the samples.

δ * (α, Σ, S, η) ≥ σ n,m 1 -α 48 max 1, min d 1 4 * , d * (1 -α) • σ n,m η , (19)
Remark. As mentioned in the introduction, by letting m go to infinity in the two-sample case, we recover the bounds of the one sample case (up to a constant factor). It is worth examining if the converse holds, i.e. if there is an argument to reduce the two-sample problem to the simpler one-sample case (this would simplify some technical aspects of the proofs, somewhat). For the upper bounds on the minimum separation distance, this is the case in some specific situations: for equal sample sizes n = m, the two-sample case can be reduced to the one-sample problem setting by pairing the samples and considering the single sample (X i -Y i ) 1≤i≤n , and one can recover this way in essence the twosample result. If Σ = S, and for general sample sizes, we can also reduce to the single sample with size min(m, n), (X i -Y i ) 1≤i≤min(m,n) , and recover again the two-sample result up to a numerical factor. However a reduction argument in the general case has eluded us. Concerning the lower bound, the argument for the two-sample case indeed hinges on a reduction the one-sample case, by considering the sub-models where one of the two sample means is known, see Section 3.4.

We now turn to the bounded setting.

Proposition 9. Assume the bounded setting holds and n, m ≥ 2. Then with probability at least 1α,

U -µ -ν 2 ≤ µ -ν q 1 + q 2 , ( 20
)
where U is defined in (12) and

q 1 (Σ, S, α) = 2 2 Σ op n + S op m u + 4Lu 3(n ∧ m) , q 2 (Σ, S, α) = 614 √ Tr Σ 2 n + √ Tr S 2 m u + 3708 L 2 u 2 (n ∧ m) 2 , with u(α) = -log α + log 2.
Thus, in the bounded setting we can guarantee that the behavior of the test is qualitatively the same as in the Gaussian setting (see e.g. Theorem 7)and this from a non-asymptotic point view, provided n ∧ m ≥ uL 2 /σ 2 , where σ 2 = Σ op .

A special case of interest is when the data lies on the sphere of radius L, i.e. X i = Y j = L a.s. In this case L 2 = Tr Σ and the above condition can be rewritten n ∧ m ≥ ud e . This situation is met in particular in the KME setting, see Section 1.3, when using a translation-invariant kernel k(z, z

′ ) = k • (z -z ′ ), in which case L 2 = k • (0).

Quantile estimation

Since we are considering the case where Σ, S can be arbitrary in this work, it is natural to assume that these are not known in advance. We study next the estimation of the quantities q 1 and q 2 , in both settings (bounded and Gaussian), in order to check Assumption 2 for our generic theorem. If we can grant that assumption, Theorem 3 guarantees that the separation distance remains qualitatively the same as in the "oracle" situation where they are known. To simplify the exposition, in this section we will present results for the one-sample problem only; similar results, although slightly more technical, can be obtained for the two-sample problem. Thus, we need to have estimators of Σ op and Tr Σ 2more precisely, of their square root.

For q 1 , we will use the empirical covariance operator Σ := Σ(X):

Σ(X) = 1 n n i=1 (X i -µ)(X i -µ) T , (21)
where µ := µ(X) is the empirical mean of the sample X.

Proposition 10 (Gaussian setting). Assume X = (X i ) 1≤i≤n are i.i.d. Gaussian vectors of covariance Σ. For u ≥ 0, with probability at least 1 -3e -u :

Σ 1 2 op -Σ 1 2 op ≤ 3 √ 2 Σ 1 2 op d e n + u n , ( 22
)
where Σ is defined in (21) and d e = Tr Σ/ Σ op .

Proposition 11 (Bounded setting). Assume that X = (X i ) 1≤i≤n are i.i.d. bounded in norm by L and with covariance Σ. For u ≥ 0, with probability at least 1 -2e -u :

Σ 1 2 op -Σ 1 2 op ≤ 4L 2 d e n + 2u n + u 3n ( 23
)
where Σ is defined in (21) and d e = Tr Σ/ Σ op .

These concentration bounds are not sharp in an asymptotic sense, where the main term for the scaling of the deviations is expected to follow that of asymptotic normality for eigenvalues of the empirical covariance operators, as in the classical results of [START_REF] Anderson | An introduction to multivariate statistical analysis[END_REF], but they are largely sufficient for our purposes (see Corollary 14 below). Some refined related nonasymptotic bounds can be found in the recent literature. In particular, [START_REF] Koltchinskii | Concentration inequalities and moment bounds for sample covariance operators[END_REF] derive nonasymptotic results for controlling Σ -Σ in the Gaussian setting, and in the centered case where µ = 0 is known. In fact, in essence the result of our technical Proposition 23 in the proof section (which is like Proposition 11 but in the centered case) can be deduced from the results of [START_REF] Koltchinskii | Concentration inequalities and moment bounds for sample covariance operators[END_REF] by elementary arguments. We decided to include a standalone proof here; while we do rely on the estimates of [START_REF] Koltchinskii | Concentration inequalities and moment bounds for sample covariance operators[END_REF] (or rather on the improved version of van Handel, 2017) for the expectation of the difference, we derive an upper bound on the deviation by a rather direct application of the Gauss-Lipschitz concentration. While [START_REF] Koltchinskii | Concentration inequalities and moment bounds for sample covariance operators[END_REF] also rely on such arguments, their proofs are much more involved, for the reason that they study the norm or the difference while we only are interested in the difference of the (root) norms here. Finally, we also mention very recent results of [START_REF] Jirak | Perturbation bounds for eigenspaces under a relative gap condition[END_REF] for sharp nonasymptotic control of spectral quantities related to Σ, which could also potentially be applied here, though it seems at first glance that a logarithmic dependence in the dimension could enter into play.

For the bounded setting (Proposition 11), the bound (23) could presumably be improved to have Σ op instead of L for the main terms. The results of Theorem 9 of [START_REF] Koltchinskii | Concentration inequalities and moment bounds for sample covariance operators[END_REF] under a sub-Gaussian assumption do not seem to be able to imply Proposition 11, see the more detailed discussion below in Section 2.4.

Turning now to q 2 , we will estimate √ Tr Σ 2 using the following statistic T := T (X), which is an unbiased estimator of Tr Σ 2 :

T (X) := 1 4n(n -1)(n -2)(n -3) i =j =k =l X i -X k , X j -X l 2 . ( 24
)
Proposition 12 (Gaussian setting). Assume X = (X i ) 1≤i≤n are i.i.d. Gaussian vectors of covariance Σ and n ≥ 4. Then for all u ≥ 0:

P T - √ Tr Σ 2 ≥ 30 Tr Σ 2 n u 2 ≤ e 4 e -u , (25)
where T is defined in (24).

Proposition 13 (Bounded setting). Assume that X = (X i ) 1≤i≤n are i.i.d. bounded in norm by L and with covariance Σ and n ≥ 4. Then for all u ≥ 0:

P T - √ Tr Σ 2 ≥ 12L 2 u n ≤ 2e -u . ( 26
)
where T is defined in (24).

Thanks to these concentration results, we can construct estimators of q 1 (Σ, α) and q 2 (Σ, α) satisfying Assumption 2. In the Gaussian setting, we give the following explicit corollary of Propositions 10 and 12; the proof is straightforward and omitted.

Corollary 14 (Gaussian setting). Consider the signal detection problem (1) and assume the Gaussian setting holds. Let α ∈ (0, 1), u = u(α) =log α+log 8, and Q 1 (α) and Q 2 (α) be the statistics defined by

Q 1 (α) = 2 Σ(X) op n u , Q 2 (u) = 32 T (X) n u ,
where Σ is defined in (21) and T in (24). Then for any Σ, provided

n max d e (Σ), u, u 4 , (27)
(we recall d e (Σ) = Tr Σ/ Σ op ), then it holds, for any distribution P ∈ P Σ :

P Q 1 (α) -q 1 (Σ, α) ≤ q 1 (Σ, α)/2 ≤ α, P Q 2 (α) -q 2 (Σ, α) ≤ q 2 (Σ, α)/2 ≤ α,
where q 1 , q 2 are as defined in ( 14),( 15) (with m = ∞).

The condition ( 27) for n is needed to grant Assumption 2: it ensures that the deviations of the estimators Σ 1 /2 op and T coming from Proposition 10 and 12 are smaller than their target quantities Σ 1 /2 op /2 and (Tr Σ 2 ) 1 /2 /2, respectively. The requirement that the size of the sample is larger than the effective dimension d e appears mild.

For the bounded setting and the signal detection problem (1), estimators Q 1 and Q 2 satisfying Assumption 2 can also be constructed in a similar way from Propositions 11 and 13 (details omitted). In the bounded setting, the quantiles q 1 and q 2 of U are composed of two terms, the first (and larger) one gives the dependence in the covariance of the distribution, the second depends on the bound L. This additional term will have to be taken into account, and the condition on n analogous to (27) will involve L. In general this will not be a problem since L or an upper bound on L is supposed to be known, as is the case for instance in the kernel setting (see the concluding discusssion in the previous section). Finally, for the two-sample test problem (2), comparable results can be obtained using the estimators Σ(Y) and T (Y); we omit the details.

Concluding remarks

A technical discussion point: Gaussian, sub-Gaussian, and bounded vectors. The utility of our systematic distinction between the Gaussian and bounded case can be disputed in the light of recent concentration literature (see e.g. [START_REF] Hsu | A tail inequality for quadratic forms of subgaussian random vectors[END_REF][START_REF] Koltchinskii | Concentration inequalities and moment bounds for sample covariance operators[END_REF] and further references therein) deriving results holding for sub-Gaussian random vectors, a seemingly more general setting emcompassing both the Gaussian and bounded settings as particular cases (since bounded variables are sub-Gaussian by Hoeffding's inequality).

This point deserves a specific discussion. The sub-Gaussianity assumption for a vector variable X (assumed centered for simplicity here) often takes the following form: for any unit vector u, denoting X u = X, u , it is assumed that X u ψ2 ≤ C Var[X u ] (where . ψ2 is the Orlicz ψ 2 -norm); or equivalently in terms of Laplace transform,

log(E[exp λ(X u )]) ≤ (C ′ ) 2 λ 2 Var[X u ]/2 for all λ ≥ 0. (28)
A key point is that the factors C or C ′ in those definitions should be independent of u, and they generally appear as global factors in the derived deviation inequalities. If the only information we have is that X is bounded a.s. by L, we see that the factors C or C ′ should be taken of the order of sup u =1 (L/ Var[X u ]) = L Σ -1 1 /2 op , which is not acceptable in a high-dimensional setting, and in particular for the application to KME described in Section 1.3, where one might expect that Σ -1 op can get arbitrarily large or even infinite. Some works (such as Spokoiny and Zhilova, 2013 and the appendix of [START_REF] Spokoiny | Basics of modern mathematical statistics[END_REF] consider settings going beyond sub-Gaussianity, i.e. when ( 28) is only required to hold for λ ≤ M -1 . This allows in principle for more general variables, e.g. chi-squared type statistics or variables admitting Bersteinor Bennett-type control of their Laplace transform, while making the constant C ′ in (28) controlled by a fixed numerical constant. Under this assumption the "first-order" terms are of the correct order, i.e. typically only depend on the variance Σ. Unfortunately, the value of M comes up into additional terms, and since its value has to be independent of u, in the bounded setting the uniformity with respect to u means that M should be again taken of the order of Σ -1 1 /2 op . To summarize, despite our best efforts we were not able to derive from existing general results, working under the (possibly extended) sub-Gaussian assumption, a concentration in the bounded setting that would not involve Σ -1 op , and this is the reason why we treated it separately with tools specific to bounded variables such as the Bousquet-Talagrand inequality. It would be of course of notable interest to obtain results under a general sub-Gaussian assumption sup u =1 X u ψ2 ≤ L, and control deviations only involving various norms of Σ for the main terms, possibly L for smaller-order terms, but not depending on Σ -1 op .

Perspectives. We finally list a few items for future developments.

-It would be interesting to obtain a version of Proposition (11) where the main term does not involve the bound L. -A recent trend of research developed "robust" exponential concentration bounds for estimators of scalars and vectors with minimal moments assumptions (see e.g. [START_REF] Lugosi | Mean estimation and regression under heavy-tailed distributions: A survey[END_REF] for a survey of recent advances). It seems a very interesting question to study if such robust procedures can be pushed to the testing setting and enjoy similar nonasymptotic controls to the Gaussian and bounded settings under much relaxed distributional assumptions. Preliminary calculations seem to indicate that the "median-of-means" (MoM) approach can be applied to U-statistics without particular problems and that Assumption 1 can be granted for MoM versions of U-statistics under the assumption of existing moments of order 4, and presumably Assumption 2 under moments of order 8. -We have analyzed here quantile estimation by direct estimation of unknown quantities coming into the quantile bounds. In practice, quantile estimation by some form of resampling procedure would be often sharper and preferred.

V. Spokoiny also made notable recent contributions to this topic [START_REF] Naumov | Bootstrap confidence sets for spectral projectors of sample covariance[END_REF][START_REF] Spokoiny | Bootstrap confidence sets under model misspecification[END_REF]. In the setting of two-sample testing where the null hypothesis is strict equality, it is possible to obtain tests with exact nonasymptotic level based on permutation tests and variations thereof; see [START_REF] Fromont | Kernels based tests with non-asymptotic bootstrap approaches for two-sample problems[END_REF] for such approaches for testing equality of distributions based on the KME methodology, and [START_REF] Kim | Minimax optimality of permutation tests[END_REF] for recent broad results on minimax optimality for the power of permutation-based tests. Estimating quantiles via bootstrap procedures is also an interesting direction to pursue in setting, in the case where the null hypothesis is based on closeness rather than equality of signals, so that exact permutation tests do not apply; Dette, Kokot, and Aue ( 2020) recently proposed nonstandard bootstrap procedures to tackle this issue. -Lower bounds establishing the optimality of the separation rates appearing have been established in the Gaussian case in Theorem 7. It would be nice find such a lower bound in the bounded case.

Proofs

The proofs of some of the technical results, first stated without justification along the text, can be found in Section 3.7. We first state a standard technical lemma which we will use several times in the following proofs.

Lemma 15. Let a ∈ R + and b ∈ R, then

-min √ b, |b| a ≤ (a 2 + b) + -a ≤ min |b|, |b| 2a . ( 29
)

Proof of Theorem 3

Let us denote D := µν . Under (H 0 ) we have D ≤ η and thus:

E H0 [T ] = P H0 U -η 2 > 2η Q 1 + 2 Q 2 ≤ P H0 U > D 2 + Dq 1 + q 2 + P H0 q 1 -Q 1 > q 1 /2 + P H0 q 2 -Q 2 > q 2 /2 ≤ 3α ,
where we have used Assumptions 1 and 2.

We will prove below that under (H 1 ), we have

P H1 D 2 -Dq 1 (u) -q 2 (u) ≤ η 2 + η2 Q 1 + 2 Q 2 ≤ 2α , (30)
which entails:

P H1 [T = 0] = P H1 U -η 2 ≤ 2η Q 1 + 2 Q 2 ≤ P H1 U ≤ D 2 -Dq 1 -q 2 + P H1 D 2 -Dq 1 -q 2 ≤ η 2 + η2 Q 1 + 2 Q 2 ≤ 3α ,
and the proof is complete. We now prove inequality (30). Let us first solve the following quadratic inequality in Z ≥ 0:

Z 2 -Zq 1 -q 2 ≥ η 2 + 3ηq 1 + 3q 2 . ( 31
)
The equation is satisfied when

Z ≥ q 1 + (2η + 3q 1 ) 2 + 16q 2 2 ;
furthermore, by Lemma 15 and the assumed inequality (11), we have that

q 1 + (2η + 3q 1 ) 2 + 16q 2 2 ≤ η + 2q 1 + min 2 √ q 2 , 2q 2 η ≤ η + δ .
Under (H 1 ), D ≥ η + δ, so D satisfies equation ( 31). We conclude by remarking that, using Assumption 2:

P H1 D 2 -Dq 1 (u) -q 2 (u) ≤ η 2 + η2 Q 1 + 2 Q 2 ≤ P η 2 + η2 Q 1 + 2 Q 2 ≥ η 2 + 3ηq 1 + 3q 2 ≤ 2α . ⊓ ⊔

Proof of Propositions 6 and 9

As much for the Gaussian case as for the bounded case, we will give concentration bounds for the statistic U defined in (12), by decomposing the statistic in four parts. Let us define:

U X := 1 n(n -1) n i,j=1 i =j X i -µ, X j -µ , U Y := 1 m(m -1) m i,j=1 i =j Y i -ν, Y j -ν , U X,Y := 1 nm n i=1 m j=1 X i -µ, Y j -ν , U * := 1 n n i=1 (X i -µ) - 1 m m j=1 (Y j -ν), µ -ν .
We have that

U = µ -ν 2 -2U * + U X + U Y -2U X,Y . ( 32
)
Gaussian setting. We first need some results on Gaussian variables. The first result is a decoupling theorem of [START_REF] Vershynin | High-Dimensional Probability: an introduction with applications to data science[END_REF].

Proposition 16 (Vershynin, 2018, Theorem 6.1.1). Let X 1 , . . . , X n be independent centered and weakly (i.e. Pettis) integrable vectors in a Hilbert space, (a ij ) 1≤i,j≤n a family of real numbers and F : R → R a convex function. Then

E   F i =j a ij X i , X j   ≤ E   F 4 i,j a ij X i , X ′ j   ,
where (X ′ i ) is an independent copy of (X i ). The following lemma is standard; see e.g. [START_REF] Birgé | An alternative point of view on Lepski's method[END_REF], Lemma 8.2.

Lemma 17. Let X a real random variable such that for all 0 < t < b -1 :

log E e tX ≤ (at) 2 1 -bt ,
where a and b are two positive constants. Then, for all t ≥ 0:

P X ≥ 2a √ t + bt ≤ e -t .
Proposition 18. Let X and Y be two independent Gaussian vectors following the distributions N (0, Σ) and N (0, S) respectively. Then for t < (

S op Σ op ) -1 /2 : log E[exp(t X, Y )] ≤ t 2 Tr(SΣ) 2(1 -t S op Σ op ) .
Using Lemma 17, for all u ≥ 0:

P X, Y ≥ 2 Tr(SΣ)u + S op Σ op u ≤ e -u .
We can now prove Proposition 6. The samples X and Y have respective distributions N (µ, Σ) and N (ν, S). We will obtain a concentration inequality for U using its decomposition (32).

Let us first find concentration inequalities for U X and U Y . Using decoupling (see Proposition 16) we have for all t < (4 Σ op ) -1 :

E[exp(tn(n -1)U X )] ≤ E exp 4t n i=1 X i -µ, n i=1 X ′ i -µ ,
where X ′ i are independent copies of the X i s. Then using Proposition 18, it holds with probability at least 1 -2e -u :

n(n -1)|U X | ≤ 4n √ 2 Tr Σ 2 u + Σ op u . ( 33
)
The same method works for U Y . The concentration of U X,Y is directly obtained using Proposition 18. Finally U * is a centered 1-dimensional Gaussian with variance (µν) T Σ n + S m (µν) and we use the classical bound

P |N | ≥ σ √ 2t ≤ 2e -t for N ∼ N (0, σ 2 ).
Thus we obtain that with probability at least 1 -8e -u :

U -µ -ν 2 ≤ 4 n -1 √ 2 Tr Σ 2 u + Σ op u + 4 m -1 √ 2 Tr S 2 u + S op u + 4 √ nm √ 2 Tr ΣSu + ( Σ op S op ) 1 2 u + 2(µ -ν) T Σ n + S m (µ -ν)u .
We conclude by upper bounding the operator norms Σ op and S op by √ Tr Σ 2 and √ Tr S 2 and for the third term we use that

(2 Tr(ΣS))

1 2 ≤ (4 Tr Σ 2 Tr S 2 ) 1 4 ≤ (Tr Σ 2 ) 1 2 + (Tr S 2 ) 1 2 .
We finally use (n -1) -1 ≤ 2n -1 for n ≥ 2 and similarly for m. It is easy to check that the fourth term is upper bounded by q 1 defined in ( 14). It just remains to use that u ≥ 1 to get u ≥ √ u and ( 13).

Bounded setting. The concentration of U is obtained in the bounded setting using a concentration inequality for degenerate U-statistics of [START_REF] Houdré | Exponential inequalities, with constants, for U-statistics of order two. In: Stochastic Inequalities and Applications[END_REF]. We present here a somewhat simplified version suited for our purpose6 .

Theorem 19 (Houdré and Reynaud-Bouret, 2003, Theorem 3.4). Let T 1 , . . . , T N be independent random variables on a probability space (Ω, F , P) with values in a Borel space (T , G). Let

U N = N i=2 i-1 j=1 g i,j (T i , T j ) ,
where g i,j : T × T → R are measurable Borelian functions satisfying

E[g i,j (T i , T j )|T i ] = E[g i,j (T i , T j )|T j ] = 0 .
Let us suppose that the following quantities are finite

A := sup t,t ′ ,i,j |g i,j (t, t ′ )| , B 2 := max    sup t,i i-1 j=1 E[g i,j (t, T j ) 2 ] , sup t,j n i=j+1 E[g i,j (T i , t) 2 ]    , C 2 := N i=2 i-1 j=1 E[g i,j (T i , T j ) 2 ] .
Then for all u > 0:

P U N ≥ 4C( √ 2u + 2 √ 2u) + 202Bu 3/2 + 196Au 2 ≤ 2.77e -u . (34)
Let us prove Proposition 9. We recall that we suppose here that the samples X and Y are both bounded by L. To obtain a deviation inequality for the statistic U , we consider separately the statistics U X + U Y -2U X,Y and then U * .

Using Theorem 19 with N = n+m, T i := X i -µ for 1 ≤ i ≤ n and

T i = Y i -ν for n + 1 ≤ i ≤ n + m, T = {u : u ≤ 4L 2 } and g ij (•, •) =      1 n(n-1) •, • , if 1 ≤ i, j ≤ n, 1 m(m-1) •, • , if n + 1 ≤ i, j ≤ n + m, -1 nm •, • , otherwise,
we get that with probability greater than 1 -5.54e -u :

|U X + U Y -2U X,Y |/2 ≤ 307 √ Tr Σ 2 n + √ Tr S 2 m u + 1854L 2 u 2 . ( 35
)
To obtain the above, we have upper bounded A, B, C by:

A ≤ 8L 2 (n ∧ m) 2 , B 2 ≤ 8L 2 (n ∧ m) 2 Σ op n + S op m , C 2 = 3 2 Tr Σ 2 n + Tr S 2 m ;
then, using that 2 √ ab ≤ a + b and that Σ op ≤ √ Tr Σ 2 , we get (35). For U * , we use Bernstein's inequality (i.e. combining Lemmas 32 and 17) to get that with probability at least 1 -2e -u , it holds:

|U * | ≤ µ -ν 2 Σ op n + S op m u + 2Lu 3n ∧ m . (36)
Combining ( 35) and ( 36), we obtain the claim of Proposition 9. ⊓ ⊔

Proof of Theorem 7

The upper bound is directly obtained using Theorem 3. Assumption 1 is satisfied as a consequence of Proposition 6. We do not consider estimation of nuisance parameters related to the covariance matrix Σ which is assumed to be fixed and known for this result; thus Assumption 2 is trivially satisfied by taking

Q 1 = q 1 (Σ, α), Q 2 = q 2 (Σ, α).
Let us now prove the lower bound (17). The following proof is an adaptation to the non-isotropic Gaussian setting of the proof of Theorem 5.1 in [START_REF] Blanchard | Minimax Euclidean separation rates for testing convex hypotheses in R d[END_REF]. Let α ∈ (0, 1), ans Σ be a positive semidefinite matrix. Without loss of generality, we can assume that Σ is diagonal:

Σ = diag(λ 1 , . . . , λ d ) with λ 1 ≥ . . . ≥ λ d > 0. Let us denote P µ,Σ the distribution of N (µ, Σ) for µ ∈ R d
and introduce the Gaussian mixture distribution:

Q n Σ := 1 2 d-1 m∈M P ⊗n m,Σ , (37)
where

M = (λ 1 v 1 h, . . . , λ d-1 v d-1 h, η)| v ∈ {-1, 1} d-1 .
We take

h 2 := (η+δ) 2 -η 2 Tr Σ 2 -λ 2 d .
Then, for all m ∈ M,

m d = η 2 + (Tr Σ 2 -λ 2 d )h 2 = η + δ .
Let ν = (0, . . . , η), it holds sup P∈H0 P ⊗n (φ = 1) + sup

P∈A δ P ⊗n (φ = 0) ≥ P ⊗n ν,Σ (φ = 1) + Q n Σ (φ = 0) ≥ 1 - 1 2 P ⊗n ν,Σ -Q n Σ TV ≥ 1 - 1 2   R d×n dQ n Σ dP ⊗n ν,Σ 2 dP ⊗n ν,Σ -1   1 2 . ( 38
)
see for instance [START_REF] Baraud | Non-asymptotic minimax rates of testing in signal detection[END_REF]. For a tensor product of Gaussian distributions with fixed, equal covariance, the empirical mean is a sufficient statistic because the Radon-Nikodym derivative of a tensor product of Gaussian measures w.r.t. the Lebesgue measure can be written for x 1 , . . . , x n ∈ R d as

dP ⊗n m,Σ dλ ⊗n (x 1 , . . . , x n ) = φ m,Σ/n (x)F Σ (x 1 , . . . , x n ) ,
where x is the mean of the x i s, φ m,Σ/n is the p.d.f. of a normal N (m, Σ/n) variable, and F Σ is a function of (x 1 , . . . , x n ) which only depends on Σ. Therefore

dQ n Σ dP ⊗n ν,Σ (x 1 , . . . , x n ) = dQ 1 Σ/n dP ν,Σ/n (x),
and thus

R d×n dQ n Σ dP ⊗n ν,Σ 2 dP ⊗n ν,Σ = R d dQ 1 Σ/n dP ν,Σ/n 2 dP ν,Σ/n .
Thus the problem boils down to studying a single Gaussian vector of covariance Σ/n; for the following we will assume n = 1 and replace at the end Σ by Σ/n. Let us compute the densities F ν and Q of these two distributions. For x ∈ R d :

F ν (x) = det Σ(2π) d -1 2 exp - 1 2λ d (x d -η) 2 d-1 i=1 exp - x 2 i 2λ i , and
Q(x) = det Σ(2π) d -1 2 exp - 1 2λ d (x d -η) 2 × 1 2 d-1 vi∈{-1,1} 1≤i≤d-1 d-1 i=1 exp - 1 2λ i (x i -hλ i v i ) 2 = det Σ(2π) d -1 2 exp - 1 2λ d (x d -η) 2 - h 2 2 d-1 i=1 λ i × d-1 i=1 exp - x 2 i 2λ i cosh(hx i ) .
Using that E cosh 2 (aZ) = exp(a 2 σ 2 ) cosh(a 2 σ 2 ) when Z ∼ N (0, σ 2 ), we have that

R d Q(x) 2 F ν (x) dx = det Σ(2π) d -1 2 exp -h 2 d-1 i=1 λ i R exp - 1 2λ d (x d -η) 2 dx d × d-1 i=1 R cosh 2 (hx i ) exp - x 2 i 2λ i dx i = exp -h 2 d-1 i=1 λ i d-1 i=1 exp(h 2 λ i ) cosh(h 2 λ i ) = d-1 i=1 cosh(h 2 λ i ).
By Taylor expansion, we obtain the bound

h 2 λ i ≤ 1 ⇒ cosh(h 2 λ i ) ≤ 1 + e 2 λ 2 i h 4 .
From this and the definition of h we deduce:

log d-1 i=1 cosh(h 2 λ i ) ≤ e 2 (Tr Σ 2 -λ 2 d )h 4 = e 2(Tr Σ 2 -λ 2 d ) (η + δ) 2 -η 2 2 .
The end of the proof follows the same steps as the proof of Theorem 5.1 of [START_REF] Blanchard | Minimax Euclidean separation rates for testing convex hypotheses in R d[END_REF]

P[φ = 1] + sup P∈A δ P[φ = 0] > α .
It follows

δ * ≥ Σ op d * -1s + η 2 -η ≥ 2 -3 2 min s Σ op (d * -1) 1 4 , s Σ op (d * -1) 1 2 η ,
and we obtain the inequality corresponding to the second part of the maximum in the right-hand side of (17) by using that s ≥ (1α) and that d * -1 ≥ 2d * /3 because d * ≥ 3.

Let us prove now that δ * Σ op . Let us consider the eigenvector e 1 associated to the maximum eigenvalue Σ op . Then P (η+δ)e1,Σ ∈ H 0 and P (η+δ)e1,Σ ∈ A δ . Let us denote λ 1 = Σ op /n, we have:

R d dP ⊗n (η+δ)e1,Σ dP ⊗n ηe1,Σ 2 dP ⊗n ηe1,Σ = R d dP (η+δ)e1,Σ/n dP ηe1,Σ/n 2 dP ηe1,Σ/n = e -δ 2 /λ1 √ λ 1 2π R exp - (x -η) 2 2λ 1 exp 2δ(x -η) λ 1 dx = exp 3δ 2 λ 1 . If δ ≤ λ1 3 log(1 + 4(1 -α) 2
), then using ( 38)

sup P∈H0 P[φ = 1] + sup P∈A δ P[φ = 0] > α .
It follows that:

δ * ≥ Σ/n op (1 -α).

Proof of Theorem 8

This proof is similar to the proof of Theorem 7, so some details will be skipped.

As in the one-sample case the upper bound is directly obtained using Theorem 3 and Proposition 6. We just additionally use the following upper bounds:

√ Tr Σ 2 n + √ Tr S 2 m ≤ √ 2 Tr Σ 2 n 2 + Tr S 2 m 2 ≤ √ 2 Tr Σ n + S m 2 ; Σ op n + S op m ≤ 2 max Σ op n , S op m ≤ 2 Σ n + S m op ,
where the last inequality holds because Σ, S are both positive semidefinite. The lower bound in the two-sample case is a direct consequence of the onesample case, by reduction to the case where one of the two sample means is known, say equal to zero. More specifically, let Σ and S be two symmetric positive semidefinite matrices, we consider again the distribution Q n Σ defined in (37). Then Then using the previous results of the proof of Theorem 7 we obtain that

δ * (α) ≥ n -1 Tr Σ 2 -λ 2 d s + η 2 1 2 -η, (39)
with s = 2 e log(1 + 4(1α) 2 ). By the same token we obtain that

δ * (α) ≥ m -1 Tr S 2 -ℓ 2 d s + η 2 1 2 -η, (40)
where ℓ d is the smallest eigenvalue of the matrix S. Because d * ≥ 3, it holds

max n -2 Tr Σ 2 -λ 2 d , m -2 Tr S 2 -ℓ 2 d ≥ 2 3 max n -2 Tr Σ 2 , m -2 Tr S 2 1 2 ≥ 1 6 Tr Σ n + S m 2 ,
and by combining ( 39) and ( 40), we obtain that

δ * (α) ≥ (2 √ 12) -1 σ min √ sd 1 4 * , s σd 1 2 * η ,
where σ = Σ/n + S/m op . We obtain (19) using again that s ≥ 1α.

The last part of the lower bound is obtained as in the one-sample case using first the distributions P ⊗n (η+δ)e1,Σ ⊗ P ⊗m 0,S and P ⊗n ηe1,Σ ⊗ P ⊗m 0,S where e 1 is still the eigenvector associated to the biggest eigenvalue of Σ. We obtain that δ * (α) Σ/n 1 /2 op . By the same token, we obtain that δ * (α) S/m 1 /2 op and conclude the proof using that 2 max( Σ/n op , S/m op ) ≥ Σ/n + S/m op .

Proof of Propositions 10 and 11

We want to obtain a concentration inequality for the estimator Σ op . To this end, we will first study the following:

Σ := 1 n n i=1 (X i -µ)(X i -µ) T , ( 41
)
where µ is the true mean of the sample X. Then we have:

Σ -Σ op = -(µ -µ)(µ -µ) T op = µ -µ 2 . ( 42
)
Gaussian setting. The concentration of Σ

1 /2
op is a consequence of the classical Lipschitz Gaussian concentration property (see e.g. Theorem 3.4 in [START_REF] Massart | Concentration Inequalities and Model Selection[END_REF].

Theorem 20 (Gaussian Lipschitz concentration). Let X = (x 1 , . . . , x d ) be a vector of i.i.d. standard Gaussian variables, and f : R d → R be a L-Lipschitz function with respect to the Euclidean norm. Then for all t ≥ 0:

P[f (X) -E[f (X)] ≥ t] ≤ e -t 2 2L 2 . ( 43
)
The following corollary is a direct consequence of that theorem (we provide a proof in Section 3.7), which will be used to control the term in (42).

Corollary 21. Let X a random Gaussian vector of distribution N (µ, Σ). Then for all u ≥ 0:

P X ≥ µ 2 + Tr Σ + 2 Σ op u ≤ e -u . (44)
We will use the results of [START_REF] Koltchinskii | Concentration inequalities and moment bounds for sample covariance operators[END_REF] giving an upper bound of the expectation of the operator norm of the deviations of Σ from its expectation. The constants come from the improved version given by van Handel (2017).

Theorem 22 [START_REF] Van Handel | Structured random matrices[END_REF]. Let X = (X i ) 1≤i≤n a sample of independent Gaussian vectors of distribution N (0, Σ), then

E Σ -Σ op ≤ Σ op (2 + √ 2) d e n + 2 d e n , (45)
where d e = Tr Σ/ Σ op and Σ is defined in equation (41).

We can now prove a concentration inequality for Σ 1 /2 op .

Proposition 23. Let X = (X i ) 1≤i≤n a sample of independent N (µ, Σ) Gaussian vectors, then for u ≥ 0, with probability at least 1 -2e -u :

Σ 1 2 op -Σ 1 2 op ≤ 2 2 Tr Σ n + 2u Σ op n , ( 46
)
where Σ is defined in (41).

Remark 24. In (46), the lower and upper bounds have been brought together, but the lower bound is in fact slightly better than the upper bound. This is due to the lower bound of the expectation where Tr Σ can be replaced by Σ op , see (49) below.

Proof. We remark that

Σ 1 2 op = sup u d =1 u t Σu = sup u d =1 1 √ n n i=1 u, X i -µ 2 1 2 = sup u d =1 sup v n =1 1 √ n n i=1 u, X i -µ v i dist ∼ sup u d =1 sup v n =1 1 √ n n i=1 u, Σ 1 2 g i v i ,
where (g i ) i=1...n are i.i.d. standard Gaussian vectors and • p for p ∈ N is defined as the Euclidean norm in R p . Let u and v be unit vectors in R d and R n respectively and f u,v : R d×n → R:

f u,v (y) := 1 √ n n i=1 u, Σ 1 2 y i v i , y ∈ R d×n .
These functions are Lipschitz: indeed for all z, y ∈ R d×n we have:

f u,v (y) -f u,v (z) = 1 √ n n i=1 u, Σ 1 2 (y i -z i ) v i ≤ 1 √ n n i=1 Σ 1 2 op y i -z i d |v i | ≤ Σ 1 2 op √ n n i=1 y i -z i 2 d = Σ 1 2 op √ n y -z d×n . (47)
A supremum of Lipschitz functions is Lipschitz, thus we can use the Gaussian Lipschitz concentration (Theorem 20), and get for all x ≥ 0:

P   Σ 1 2 op -E Σ 1 2 op ≥ 2x Σ op n   ≤ e -x , (48)
with the same control for lower deviations.

It remains to upper bound E Σ

1 /2 op -Σ 1 /2
op . For one direction, using Jensen's and triangle inequalities and inequality (29), we get:

E Σ 1 2 op -Σ 1 2 op ≤ Σ op + E Σ -Σ op - Σ op ≤ min   E Σ -Σ op , E Σ -Σ op 2 Σ op   ≤ 2 2 Tr Σ n .
For the last inequality, we have used Theorem 22 for the expectation, then the fact that min (a

√ x + bx) 1 /2 , (a √ x + bx)/2 ≤ max( √ a + b, (a + b)/2) √ x where a = 2 + √ 2, b = 2 and x = d e /n
. This is achieved by treating cases x ≤ 1 and x ≥ 1 separately.

For the other direction, a reformulation of ( 48) is that there exists a random variable g ∼ Exp(1) such that:

Σ 1 2 op ≤ E Σ 1 2 op + 2g Σ op n .
Taking the square then the expectation and then applying Jensen's inequality to the concave function x → (a + b √ x) 2 (a, b ≥ 0), we obtain:

Σ op ≤ E Σ op ≤ E g   E Σ 1 2 op + 2g Σ op n 2   ≤ E Σ 1 2 op + 2 Σ op n 2 ,
and thus

E Σ 1 2 op -Σ 1 2 op ≥ - 2 Σ op n ≥ -2 2 Tr Σ n . ( 49
)
Proof of Proposition 10. It holds

Σ 1 2 op -Σ 1 2 op ≤ Σ 1 2 op -Σ 1 2 op + Σ -Σ 1 2 op .
Then, from (42):

Σ -Σ 1 2
op ≤ µµ . According to Proposition 23 and Corollary 21, we obtain that for u ≥ 0, with probability at least 1 -3e -u :

Σ 1 2 op -Σ 1 2 op ≤ 2 2 Tr Σ n + 2u Σ op n + Tr Σ n + 2 Σ op u n .
So, for u ≥ 0, with probability at least 1 -3e -u :

Σ 1 2 op -Σ 1 2 op ≤ 3 2 Tr Σ n + 2 2u Σ op n .
Bounded setting. We first recall the following concentration result for bounded random vectors in the formulation of [START_REF] Bousquet | A Bennett concentration inequality and its application to suprema of empirical processes[END_REF].

Theorem 25 (Talagrand-Bousquet inequality). Assume (X i ) 1≤i≤n are i.i.d. with marginal distribution P. Let F be a countable set of functions from X to R and assume that all functions f in F are P-measurable, square-integrable, bounded by M and satisfy E[f ] = 0. Then we denote

Z = sup f ∈F n i=1 f (X i ).
Let σ be a positive real number such that σ 2 ≥ sup f ∈F Var[f (X 1 )]. Then for all u ≥ 0, ε > 0 we have:

P Z ≥ E[Z](1 + ε) + √ 2unσ 2 + M u 3 (1 + ε -1 ) ≤ e -u .
where we have used for the variance term:

sup u,v∈B1 E v, Z i Z T i -Σ u 2 ≤ sup u,v∈B1 E v, Z i 2 Z i , u 2 ≤ 4nL 2 Σ op .
We conclude using the upper bound of the expectation from Lemma 27.

Proof of Proposition 11. As in the Gaussian case, we have:

Σ 1 2 op -Σ 1 2 op ≤ Σ 1 2 op -Σ 1 2 op + Σ -Σ 1 2 op .
From Lemma 15 and Proposition 29, we have with probability at least 1e -u :

Σ 1 2 op -Σ 1 2 op ≤ 4L Tr Σ n Σ op + 16 3 L 2 u n ,
where we have used that:

Var[ Z 1 2 ] n ≤ 2L √ Tr Σ √ n . Using Σ -Σ 1 2
op ≤ µµ , and according to Corollary 26, we obtain that for u ≥ 0, with probability at least 1 -2e -u :

Σ 1 2 op -Σ 1 2 op ≤ 4L Tr Σ n Σ op + 16 3 L 2 u n + 2 Tr Σ n + 2 Σ op u n + 4Lu 3n ≤ 8L Tr Σ n Σ op + 4L 2u n + u 3n ,
where we have used for the last inequality that Σ op ≤ 4L 2 . ⊓ ⊔

Proof of Propositions 12 and 13

From a sample X = (X i ) 1≤i≤n of i.i.d. random vectors, we want to estimate Tr Σ 2 where Σ is their common covariance matrix. The statistic T defined in ( 24) is an unbiased estimator of Tr Σ 2 . This statistic is also invariant by translation. constant (∇ µ τ = 0).

If we denote S n the set of permutations of {1, . . . , n}, T can be rewritten as:

T = 1 n! σ∈Sn 1 ⌊n/4⌋ ⌊n/4⌋ i=1 1 4 X σ(4i) -X σ(4i-2) , X σ(4i-1) -X σ(4i-3) 2 ; (52)
namely by symmetry, all the 4-tuples appear the same number of times in the right-hand side, so we just need to divide by the number of terms to obtain the identity (52). We will use this decomposition to obtain a concentration of the statistic T for the Gaussian case and the bounded case, since the inner sum for each fixed permutation is a sum of ⌊n/4⌋ i.i.d. terms.

Gaussian setting. Because the statistic is invariant by translation we can assume without loss of generality that µ = 0. To obtain a deviation inequality for T 1 /2 , we will first find a concentration inequality for T and then use Lemma 15. We obtain concentration via control of moments of T , so we first need some upper bounds on Gaussian moments. The following lemma is proved in Section 3.7.

Lemma 30. Let Z i := X 1 i -X 3 i , X 2 i -X 4 i 2 /4
, where X j i for i = 1, . . . , m and 1 ≤ j ≤ 4 are i.i.d. Gaussian random vectors N (0, Σ). Then for all q ∈ N:

E   1 m m i=1 Z i -Tr Σ 2 2q   ≤ 4 √ 2φq 2 Tr Σ 2 √ m 2q , (53)
where φ = (1 + √ 5)/2 is the golden ratio.

We deduce from this lemma a concentration inequality for T .

Proposition 31. Let (X i ) 1≤i≤n , n ≥ 4 be i.i.d. random vectors with distribution N (µ, Σ). Then for all u ≥ 0:

P T -Tr Σ 2 ≥ 30 u 2 Tr Σ 2 √ n ≤ e 4 e -u , (54)
where T is defined in (24).

Proof. Using Lemma 30, (52) and the convexity of the function x → x 2q , we can upper bound the moments of T :

E T -Tr Σ 2 2q ≤ 4 √ 2φq 2 Tr Σ 2 ⌊n/4⌋ 2q . ( 55
)
Let t ≥ 0 and q ∈ N, then by Markov's inequality

P T -Tr Σ 2 ≥ t ≤ t -2q E T -Tr Σ 2 2q . (56)
Let us choose q as:

q =     e -1 2 √ φ2 1 4 t 1 2 Tr Σ 2 ⌊n/4⌋ -1 2     ,
so that ( 55), ( 56) entail

P T -Tr Σ 2 ≥ t ≤ e -4q .
Let us now take

t = e 2 √ 2φ 4 u 2 Tr Σ 2 ⌊n/4⌋ ≤ 30 u 2 Tr Σ 2 √ n ,
where we have used ⌊n/4⌋ ≥ n/7 for n ≥ 4; we obtain that for all u ≥ 0:

P T -Tr Σ 2 ≥ 30 u 2 Tr Σ 2 √ n ≤ e 4 e -u .

⊓ ⊔

Proposition 12 directly follows from Proposition 31 and Lemma 15.

Bounded setting. As in the Gaussian case, we first obtain a concentration inequality for T and then using Lemma 15, we obtain one for T 1 /2 . We will need the following classical Bernstein's inequality (see for instance Vershynin, 2018, Exercise 2.8.5 for the version below) which gives an upper bound on the Laplace transform of the sum of bounded random variables.

Lemma 32 (Bernstein's inequality). Let (X i ) 1≤i≤m be i.i.d. real centered random variables bounded by B such that

E X 2 1 ≤ σ 2 .
Then for all t < 3/B: log E e t Xi ≤ 1 2

mσ 2 t 2 1 -Bt/3 .
Via Bernstein's inequality we obtain the following result.

Proposition 33. Let (X i ) 1 ≤i≤n , n ≥ 4 be i.i.d. Hilbert-valued random variables with norm bounded by L and covariance Σ, and T defined by (24). Then for all t ≥ 0:

P T -Tr Σ 2 ≥ 8L 2 Tr Σ 2 t n + 10L 4 t n ≤ 2e -t . ( 57
)
where T is defined in (24).

Proof. Let X, X ′ , Y , Y ′ be i.i.d. Hilbert-valued random vectors of expectation µ, covariance Σ and with norm bounded by L, and

Z := X -Y, X ′ -Y ′ 2 /4. Then it holds 0 ≤ Z ≤ 4L 4 , E[Z] = Tr Σ 2 and |Z -E[Z]| ≤ 4L 4 ; Var[Z] ≤ 4L 4 E[Z] = 4L 4 Tr Σ 2 .
Now using the convexity of the exponential function, (52) and then Lemma 32, we can upper bound the Laplace transform of T as follows:

log E e For the other side, we proceed analogously:

T - √ Tr Σ 2 ≥ Tr Σ 2 -8L 2 Tr Σ 2 t n + - √ Tr Σ 2 - 10L 4 t n ≥ -8L 2 t n -L 2 10t n ≥ -12L 2 t n .
⊓ ⊔

Additional proofs

Proof of Lemma 15. This Lemma completes the Lemma 6.1.3 of [START_REF] Blanchard | Minimax Euclidean separation rates for testing convex hypotheses in R d[END_REF]. This is its complete proof. Let a in R + , it is well known that for b ≥ -a 2 : a -|b| ≤ a 2 + b ≤ a + |b| .

On the other hand, suppose that b ≥ 0, the Taylor expansion of the function b → √ a 2 + ba gives that there exists c ∈ (0, b) such that:

a 2 + b -a = b 2 √ a 2 + c ≤ b 2a .
Suppose now that 0 ≥ b ≥ -a 2 , then

a 2 + b ≥ a + b a ⇔ b ≥ 2b + b 2 a 2 ⇔ b ≥ -a 2 .
The equation ( 29) is still true when b < -a 2 because then:

-a ≥ -|b| ≥ -|b| a .

⊓ ⊔

Proof of Proposition 18. Let g be a standard Gaussian random vector in R d , and U T DU be the singular value decomposition of the matrix S 1 /2 ΣS 1 /2 where D = diag(λ , . . . , λ d ). Then we have the following equalities in distribution

Y T ΣY dist ∼ g T S 1 2 ΣS 1 2 g dist ∼ g T U T DU g dist ∼ g T Dg .
The last equality is a consequence of the invariance by rotation of Gaussian vectors. Then for t < 1/ Σ op S op :

E e t X,Y = E e t 2 Σ 1 2 Y 2 2 = E e t 2 g T Dg 2 = E exp t 2 2 d i=1 λ i g 2 i .
Using the independence of the coordinates and thatlog(1x) ≤ x 1-x ≤ x 1-√ x for x < 1 (the first inequality can easily be checked by termwise power series comparison), we obtain:

log E e t X,Y = n i=1 - 1 2 log 1 -t λ i ≤ n i=1 1 2 t 2 λ i 1 -t 2 λ i ≤ 1 2 t 2 Tr(S 1 2 ΣS 1 2 ) 1 -t S 1 2 ΣS 1 2 1 2 op .
We conclude using that Tr(S This function f is Lipschitz with constant Σ 1 2 op . We conclude using Theorem 20 and Jensen's inequality:

E X d ≤ µ 2 d + Tr Σ .

⊓ ⊔

Proof of Corollary 26. We apply Theorem 25, with ε = 1 and the set of functions F = {f u } u H=1 where f u : x ∈ H → x, u H for u ∈ H. We can find a countable subset of the unit sphere because H is separable. Then

Z = sup u H =1 n i=1 X i -µ, u H = n i=1 X i -µ H .
We conclude using that for all u in the unit sphere of H, Var[ X iµ, u H ] ≤ Σ op and | X iµ, u H | ≤ 2L a.s. We use Jensen's inequality to upper bound the expectation:

E[Z] ≤ (n Tr Σ) 1 2 . ⊓ ⊔
Proof of Lemma 27. We upper bound the operator norm with the Frobenius norm. We denote in this proof Z i := X iµ. It holds:

E Σ -Σ op ≤ E Tr Σ -Σ 2 ≤   E Tr 1 n 2 i (Z i Z T i ) 2 + i =j Z i Z T i Z j Z T j -ΣΣ -Σ Σ + Σ 2   1 2 = E Z 4 n - Tr Σ 2 n 1 2 = Var[ Z 2 ] n ≤ 2L √ Tr Σ √ n .

⊓ ⊔

Proof Lemma 30. First let us remark that if X and X ′ are independent N (0, Σ) Gaussian vectors, then

X, X ′ dist ∼ d i=1 λ i g i g ′ i ,
where g i and g ′ i are independent standard Gaussian random variables and the λ i s are the eigenvalues of Σ. Then for q ∈ N, recalling E g 2q i = (2q!)/(2 q q!), E X, X ′ 2q = where we have used (2p)! ≤ 2 2p p! 2 . Using this bound, we upper bound the moments of the Z i s:

E[Z q i ] = 2 -2q E X 1 i -X 3 i , X 2 i -X 4 i 2q ≤ (2q)!(Tr Σ 2 ) q .
We now upper bound the moments of Z i -Tr Σ. Let Z ′ i be an independent copy of Z i , then since E[Z ′ i ] = Tr Σ 2 , by Jensen's inequality

E Z i -Tr Σ 2 2q ≤ E (Z i -Z ′ i ) 2q ≤ 2 2q E Z 2q i ≤ (4q)!(2 Tr Σ 2 ) 2q .
For the odd moments we use that the function (•) 2q+1 is increasing:

-(Tr Σ 2 ) 2q+1 ≤ E Z i -Tr Σ 2 2q+1 ≤ E Z 2q+1 i ≤ (4q + 2)!(Tr Σ 2 ) 2q+1
, so for all q ≥ 0: E Z i -Tr Σ 2 q ≤ (2q)!(2 Tr Σ 2 ) q .

(59)

It remains to upper bound the moments of the sum: 1.

Let us count the number of terms in this last sum. Consider first that we have k non-null terms (p i1 , . . . , p i k ). Their sum is equal to 2q but because these terms are strictly greater than 1, we also have that (p i1 -2) + . . . + (p i k -2) = 2q -2k, where all terms of this sum are nonnegative. The number of k-partitions of 2q-2k is (2q-2k)+(k-1)

k-1

= 2q-k-1 k-1
and then the number of terms in the sum is equal to:

m k=0 m k 2q -k -1 k -1 = m∧q k=0 m k 2q -k -1 k -1 ≤ m q q k=0 2q -k -1 k -1 = m q F (2q -1) ≤ m q φ 2q ,
where F (•) is the Fibonacci sequence and φ = (1 + √ 5)/2 is the golden ratio. So using that (2q)! ≤ (2q) q q q we obtain that

E   1 m m i=1 Z i -Tr Σ 2 2q   ≤ (2φ 2 ) q Tr Σ 2 √ m 2q (2q) 4q . ( 60
)
⊓ ⊔

  ) where σ 2 n := Σ op /n, and d * := Tr Σ 2 / Σ 2 op . (The symbol indicates inequality up to a numerical factor).

  where σ 2 n,m := M n,m op , and d * := Tr M 2 n,m /σ 4 n,m , for M n,m := Σ/n + S/m. (The symbol indicates inequality up to a numerical factor).

  µ d = f (g) , where for y ∈ R d : f (y) = Σ 1 2 y + µ d .

≤

  (2q)!(Tr Σ 2 ) q ,

  Tr Σ 2 t 2 1 -4L 4 t/(3⌊n/4⌋) ,for all t such that the right-hand sise is well defined, i.e. the denominator is strictly positive. Now using Lemma 17, and ⌊n/4⌋ ≥ n/7 for n ≥ 4, for all t ≥ 0 it holdsP T -Tr Σ 2 ≥ 8L 2 Tr Σ 2 tProof of Proposition 13. Assuming the event entering into (58) holds, we will use the inequalities of Lemma 15:

		t T ≤	1 2⌊n/4⌋	4L 4 n +	10L 4 t n	≤ 2e -t .	(58)
	T -	√ Tr Σ 2 ≤ Tr Σ 2 + 8L 2 Tr Σ 2 t n	-	√ Tr Σ 2 +	10L 4 t n
		≤ 4L 2 t n	+ L 2 10t n	≤ 8L 2 t n	.

With the notation Σ p we mean p-Schatten norm. We will freely use in the paper the equivalent notation Σ ∞ = Σ op , Σ 1 = Tr(Σ), Σ 2 2 = Tr(Σ 2 ).

that is, the real random variable Φ(Z) is integrable, which guarantees that the integral of Φ(Z) is well-defined in a strong sense as an element of the Hilbert space; see e.g.[START_REF] Cohn | Measure theory / Donald L. Cohn[END_REF].

In the original result the u deviation term involves an additional constant D and we simply use D ≤ C here.

Acknowledgements. GB acknowledges support from: Deutsche Forschungsgemeinschaft (DFG) -SFB1294/1 -318763901; Agence Nationale de la Recherche (ANR), ANR-19-CHIA-0021-01 "BiSCottE"; the Franco-German University (UFA) through the binational Doktorandenkolleg CDFA 01-18.

The following corollary is a direct consequence of Theorem 25. Some refinement of this result in the same vein (including two-sided deviation control in the uncentered case) can be found in [START_REF] Marienwald | High-dimensional multi-task averaging and application to kernel mean embedding[END_REF] (Proposition 6.2 and Corollary 6.3).

Corollary 26. Let X i for i = 1, . . . , n i.i.d. random vectors bounded by L with expectation µ, covariance Σ in a separable Hilbert space H. Then for u ≥ 0, with probability at least 1e -u :

where Σ is defined in (41).

Remark 28. Using the boundedness of the variables we can upper bound this variance:

Proposition 29. Let (X i ) 1≤i≤n be i.i.d. random vectors in a separable Hilbert space H, with norm bounded by L and covariance Σ, then for any for u ≥ 1, with probability at least 1e -u :

where Σ is defined in (41).

Proof. We denote in this proof Z i := X iµ for 1 ≤ i ≤ n. Let us first remark that if B 1 is the unit ball of H, then:

Since the variables X i have norm bounded by L, it can be assumed equivalently that they take their values in B L = LB 1 , and it holds

is continuous, and the Hilbert space H is separable, the uncountable set B 1 can be replaced by a countable dense subset. Thus we can apply Theorem 25, and obtain that with probability at least 1e -x :


Gaëtan Desrues

Delphine Feuerstein

Thierry Legay

Serge Cazeau

Maxime Sermesant
email: maxime.sermesant@inria.fr

Personal-by-design: a 3D Electromechanical Model of the Heart Tailored for Personalisation

Keywords: Personalisation, Digital twin, Cardiac electromechanical model, Electrophysiology, Electrocardiogram

In this work we present a coupled electromechanical model of the heart for patient-specific simulations, and in particular cardiac resynchronisation therapy. To this end, we propose a fast fully autonomous and flexible pipeline to generate and optimise the data required to run the mechanical simulation. After the meshing of the biventricular segmentation image and the construction of the associated fibres arrangement, we compute the electrical potential propagation in the myocardial tissue from selected onset points on the endocardium. We generate a 12-lead electrocardiogram corresponding to the latter activation map by extrapolating the electrical potential on a virtual torso. This electrical activation is coupled to a mechanical model, featuring a small set of interpretable parameters. We also propose an efficient algorithm to optimise the model parameters, based on patient data. The whole pipeline including a cardiac cycle is computed in 30 minutes, enabling to use this digital twin for diagnosis and therapy planning.

Introduction

Building patient-specific 3D models of the heart can help improving diagnosis and therapy selection for various cardiac diseases. They can for instance help interventional cardiologists in choosing the best pacing method and anticipate the patient response based on the chosen pacing sites [START_REF] Sermesant | Patient-specific electromechanical models of the heart for the prediction of pacing acute effects in CRT: A preliminary clinical validation[END_REF]. However, personalisation, i.e. adapting a generic cardiac model to a specific patient, is still challenging due to theoretical (identifiability) and practical (computation time) issues.

The model presented here combines different modelling approaches and tools that are well suited for personalisation, and results in a computational time compatible with clinical constraints. In this manuscript, we detail the different elements, how we can adjust their parameters, and present simulation results to evaluate their feasibility. The computational cost of each step is also detailed. This pipeline is based on various previous studies that have been combined to propose for the first time in the team a complete and fast method for running a personalised patient-specific simulation of a beating heart, including:

-a local anatomical correction tool, -a realistic approach to simulate the Purkinje activation and electrical propagation, -algorithms to quickly generate body surface electrocardiograms, -a ready-to-use mechanical model, -a sensitivity analysis to assess the model's observability and -algorithms to personalise the models' parameters based on patient data. Leads location on a virtual torso. f. 12-lead ECG. This pipeline generates all the data required to run the mechanical simulation

Anatomical Model

We suppose here that the patient biventricular myocardium was already segmented from 3D images. This is more and more available thanks to deep learning approaches. Note that the atria are not included in our 3D model, but their interactions with the ventricles are modelled in different ways: a delay is introduced in the electrical model to simulate the propagation of the electrical signal from the sino-atrial to the auriculo-ventricular node; the atrial pressure is modelled as a sigmoid curve; and the atrial contraction is taken into account in the boundary conditions, as a spring force on the basal tetrahedra.

Meshing and Labelling

The binary mask from the segmentation (Fig. 2a) is resampled into 1 × 1 × 1mm 3 voxels and contains the labels for left and right ventricles. The labels for endocardium, myocardium and epicardium are built with a ray-tracing method, emerging from each ventricle barycentre. A tetrahedral mesh is built using the remeshing software MMG [START_REF] Dapogny | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF] (Fig. 2b), resulting in a mesh of approximately 90k tetrahedra and 20k vertices. In addition to the generated mesh, tools were developed to virtually modify the anatomy, in order to increase the available number of healthy and pathological cases, and to correct any inaccuracies from the segmentation method. For instance, on Fig. 1a, the right ventricle volume has been reduced by 20%. Several pathological cases could be generated similarly from a healthy heart geometry, for example dilated or hypertrophic cardiomyopathies.

Cardiac Fibres

The spatial orientation of the muscle fibres plays a major role in the excitation and contraction of the heart [START_REF] Streeter | Fiber orientation in the canine left ventricle during diastole and systole[END_REF]. We suppose that the fibre helix angle varies from -α to +α across the myocardial wall and remains constant throughout the cardiac cycle. We neglect the angle variations in any other direction.

To compute the fibres direction, we need to build the local coordinate system in each voxel. To this end, we first compute the tissue thickness using a simple diffusion equation on the endocardium and epicardium as in [START_REF] Yezzi | An eulerian pde approach for computing tissue thickness[END_REF]. The resultant vector field gives the radial direction. The other directions are then easy to obtain. The fibre directions are given by the circumferential basis vector, rotated by an angle α f ∈ [-α, +α] along the radial direction (Fig. 1b). In this framework, the angle variation is chosen to linearly increase from the endocardium (α f = -80 o ) to the epicardium (α f = 80 o ). Note that the septum is considered as left ventricle for the fibres generation.

AHA Regions

The 17 American Heart Association (AHA) segments are defined as a standardised segmentation of the left ventricle. We can extend this to 29 regions (Fig. 2b) for the whole myocardial tissue. These segments are mostly used in the mechanical simulation. They can also give a good insight on regional properties such as stress and strain, thus help in assessing intraventricular dyssynchrony.

To build these regions for each ventricle, we let an ellipsoid fit the ventricle mesh vertices and use the spherical coordinates to generate the regions on the ellipsoid. Finally, the segments are projected onto the tetrahedral mesh [START_REF] Cedilnik | Fast personalized electrophysiological models from computed tomography images for ventricular tachycardia ablation planning[END_REF]. At the base of the mesh, a stiffer and less conductive region is introduced to model the valves fibrous tissue (beige region at the top of the mesh in Fig. 2b).

Electrophysiology simulation

Electrical activation

At each cardiac cycle, the mechanical contraction of the heart is driven by the electrical activation. The electrical signal emerging from the auriculo-ventricular node is conducted along the bundle of His and the Purkinje network to the endocardium (Fig. 1c). The numerous termination points of the Purkinje fibres on the endocardium are reduced to a dozen of points (Fig. 1d) and selected matching a real endocardial map (Fig. 3). To account for the fast potential propagation due to the Purkinje fibres, a thin endocardial layer is assigned a higher isotropic conductivity. These onset points are activated between 0 and 21ms, t = 0 marking the beginning of the QRS complex. In the remaining tissue, the potential propagation is computed using a fast marching method and solving the anisotropic Eikonal equation [START_REF] Sermesant | An anisotropic multi-front fast marching method for realtime simulation of cardiac electrophysiology[END_REF] for each mesh vertex:

σ √ ∇T T D∇T = 1
After a short period of time (duration of the QRS), all the myocardial cells are activated and contract until the repolarisation wave arrives. The action potential duration (APD) is linearly interpolated from the endocardium to the epicardium (shorter on the epicardium), using a wall depth map.

ECG generation

The ECG records the bio-electric activity of the cardiac cells by measuring the electrical potential at the 6 precordial leads placed on the patient's torso and 3 (or 4) limb leads (Fig. 1e). In order to simulate this electrical activity, we consider each voxel of the image as a dipole of current density j eq = -σ∇v where ∇v is the spacial gradient of the potential v [START_REF] Cedilnik | Fully Automated Electrophysiological Model Personalisation Framework from CT Imaging[END_REF]. By using the chain rule, we obtain:

j eq = -σ ∂v ∂T ∇T
where σ is the local conductivity, ∇T the gradient of the activation map and ∂v ∂T is given by solving the Mitchell-Schaeffer model [START_REF] Mitchell | A two-current model for the dynamics of cardiac membrane[END_REF] using a forward Euler scheme. We suppose here that the body is a homogeneous material with constant conductivity σ T . The electrical potential at a distance r from the source is developed in [START_REF] Malmivuo | Bioelectromagnetism -Principles and Applications of Bioelectric and Biomagnetic Fields[END_REF] and is given by:

Φ(r) = 1 4πσ T V j eq • ∇ 1 r dV
We finally derive the potential contribution of each voxel with the Einthoven triangle to obtain the augmented leads and plot the electrocardiogram Fig. 1f.

Cardiac mechanics

After the electrical depolarisation, the muscle cells in the heart release ions that lead to sarcomere shortening and active muscle contraction. These cardiomyocites are surrounded by an extracellular matrix that accounts for the passive material. We use here the Bestel-Clement-Sorine [START_REF] Bestel | A biomechanical model of muscle contraction[END_REF] model, further improved by [START_REF] Chapelle | An energy-preserving muscle tissue model: formulation and compatible discretizations[END_REF], based on a multi-scale physiological description of the myocardial muscle function.

Fig. 4: Mesh deformation and active stress over a cardiac cycle

Bestel-Clement-Sorine model

In this model, the heart is described as a passive isotropic Mooney-Rivlin material with 3 parameters, which accounts for the elasticity and friction in the cardiac extracellular matrix (mainly collagen) surrounding the fibres.

The electrical stimulation is derived from the Eikonal activation map computed beforehand and is coupled to the active orthotropic contraction part, which has 7 main parameters and accounts for the active stress (Fig. 4) along the cardiac fibres and the elasticity between sarcomeres and Z-discs. Note that combining the Mooney-Rivlin isotropic material with the elasticity of the Z-discs, the material is globally considered transversely isotropic.

The cardiac cycle is decomposed into four phases: filling, isovolumetric contraction, ejection and isovolumetric relaxation (Fig. 5 and Fig. 6). We couple a haemodynamic model implementing those phases and use the four-element Windkessel model to compute the arterial pressures.

Moreover, we performed a sensivity analysis (Table 1) on the mechanical model using a Latin Hypercube Sampling method, coupled with the computation of the Partial Ranked Correlation Coefficient [START_REF] Marino | A methodology for performing global uncertainty and sensitivity analysis in systems biology[END_REF] for a collection of relevant model outputs, selected from [START_REF] Cazeau | Statistical ranking of electromechanical dyssynchrony parameters for CRT[END_REF].

The four first parameters evaluated are the four elements of the Windkessel model and thus are not, or poorly, correlated with non-ejection indexes such as LPEI or the isovolumetric contraction time. Since they control the shape of the arterial and ventricular pressures during the ejection, a strong correlation on the The closing time of the aortic and pulmonary valves is also depending on those pressures, thus the correlation observed on the isovolumetric relaxation time. We can notice that the relaxation rate, the myocardial passive stiffness and the fibres extremum angle do not seem to impact these model outputs. Finally, the contractility is strongly correlated with all the considered indexes.

Personalisation method

This model can be customised with only a few set of parameters and is therefore adapted to personalisation. To this end, we use the Covariance Matrix Adaptation Evolution Strategy (CMA-ES), an evolutionary algorithm for difficult non-linear non-convex black-box optimisation problems [START_REF] Hansen | The CMA evolution strategy: A comparing review[END_REF]. The parameters are sampled according to a multi-variate normal distribution. The covariance matrix of this distribution is iteratively updated such as the new set of parameters minimises the cost function. The latter is built such that the output of the model with the current parameters fits the patient data. Moreover, this personalisation method is embarrassingly parallel, as all simulations are independent within an iteration.

Preliminary results on the mechanical model are presented in Fig. 7. For these results, only the Windkessel and the contractility parameters have been included in the CMA-ES framework, thus spanning a five-dimensional parameter space.

Anatomy

Depending on the available patient data, it is possible to generate a mesh that meets several constraints: axis length, volumes, wall thickness, etc.. If the anatomy of the heart is poorly defined, it is possible to adjust the geometry of mesh based on the patient's cardiomyopathy.

Electrophysiology

In order to match the QRS axis and the QRS duration of each patient, the following parameters are being optimised:

-depolarisation delay of the Purkinje breakthrough points (one parameter for all termination points or one per point), -position on the endocardium of these points (two parameters per point) and -endocardial and myocardial conductivities (two parameters).

Mechanics

All the indexes presented in Table 1 (columns of the table) will be used for the mechanical personalisation. Again, the following parameters (rows of the table) are being optimised:

-four elements of the Windkessel model. Note that the parameters for the right ventricle are adjusted according to those of the left ventricle and are not counted here (four parameters), -contraction and relaxation parameters (two parameters), -Mooney-Rivlin material parameters (three parameters) and -active parameters, sarcomere contractility and viscosity, Z-discs elasticity (three parameters). This personalisation method has shown great promises in our preliminary tests, in terms of computational cost and global error. We are currently implementing the CMA-ES method for both the electrophysiological and mechanical frameworks on a computer cluster. This algorithm converges in only a few iterations, for example 30 iterations if the population size is ten times the dimension of the parameter space. This space is spanned by a user-defined range of realistic values for each parameter.

With this configuration, the electrophysiological model can be personalised in less than two hours and the mechanical model in less than twelve hours on a computer cluster.

Conclusion

In this manuscript, we have presented a fast (see Table 2) pipeline for patientspecific electromechanical simulations, using a model relying on a small set of interpretable parameters, allowing for extensive personalisation. This pipeline is built upon three submodels: an anatomical one, allowing local corrections, as well as mesh and fibres generation; an electrophysiological one, generating ECGs and simulating the complex activation and propagation of the electrical potential; and a mechanical submodel. We also lead a sensitivity analysis to assess the model's observability, and developed methods to efficiently fit the model output to patient data. Such approach opens up possibilities in order to use cardiac modelling within clinical applications.

Fig. 1 :

 1 Fig. 1: Pipeline for patient-specific simulation. a. Initial and modified geometries. b. Fibre directions. c. Activation points and activation times on the endocardial surfaces. d. Activation map. e.Leads location on a virtual torso. f. 12-lead ECG. This pipeline generates all the data required to run the mechanical simulation

Fig. 2 :

 2 Fig. 2: Generated topologies

Fig. 3 :

 3 Fig. 3: Real endocardial activation map (courtesy of Dr.Mouhoub) and their computed homologous. The computed map (white background) shows both ventricles endocardium (no atria) and the fast isotropic propagation wave.

Fig. 5 :Fig. 6 :

 56 Fig. 5: Left ventricle pressure and volume curves extracted from the simulation

Fig. 7 :

 7 Fig. 7: Convergence of the CMA-ES features error (%) on a mechanical personalisation for two (left) and four (right) patient indexes

Table 1 :

 1 Results of the Sensivity Analysis. EF: Ejection Fraction, LVET: Left Ventricle Ejection Time, LPEI: Left ventricle Pre-Ejection Interval, SD: Systole Duration, Iso.CT: Isovolumetric Contraction Time, Iso.RT: Isovolumetric Relaxation Time, IVD: InterVentricular Delay. R p : arterial peripheral resistance, Z c : characteristic impedance of the artery, LL: total arterial inertance, C: total arterial compliance, K rs : relaxation rate, K atp : contraction rate, K: stiffness, σ: contractility, α: maximum fibres angle. Green: p-value<0.001, Yellow: p-value<0.01, Orange: p-value<0.1, Red: p-value>0.1.

		EF	LVET	LPEI	LPEI/LVET	SD	Iso.CT	Iso.RT	IVD
	R p	-0.42222	0.23128	0.02273	-0.11963	0.18141	0.02645	0.64766	-0.00766
	Z c	-0.03788	0.21657	0.00576	-0.10319	0.14087	0.00125	0.05784	-0.00157
	LL	0.24441	0.58033	-0.00847	-0.36903	0.34962	-0.01144	-0.349	-0.00948
	τ	0.32062	-0.29936	0.10974	0.304	-0.19011	0.10257	-0.48444	0.18499
	K rs	-0.02571	-0.03178	0.01345	0.01319	-0.02072	0.00205	-0.01778	0.02944
	K atp	0.07897	-0.02247	-0.4114	-0.39651	-0.25082	-0.41028	-0.00922	-0.30771
	K	-0.02688	0.00409	0.03267	0.02106	0.00687	0.03839	0.02471	0.06666
	σ	0.72809	-0.25892	-0.81655	-0.64763	-0.62813	-0.81896	-0.32104	-0.79371
	α	-0.02548	-0.03125	0.02816	0.00826	-0.00697	0.01788	-0.00138	0.096
					1

Table 2 :

 2 Computation time of the different steps of the pipeline. Computed on an Intel Core i5-4430 CPU @ 3.00GHz machine.

	Anatomical transformations	Meshing Labelling	Fibres direction	Activation map	ECG generation	One cardiac cycle
	10s-1min	30s	1min	15s	30s	1min	20min

Acknowledgments

This work has been supported by the French government through the National Research Agency (ANR) Investments in the Future 3IA Côte d'Azur (ANR-19-P3IA-000) and by Microport CRM funding.