<figure><img src="image_1.png" /></figure>

HAL
open science

# Optimal multiple change-point detection for high-dimensional data

Emmanuel Pilliat, Alexandra Carpentier, Nicolas Verzelen

► To cite this version:

Emmanuel Pilliat, Alexandra Carpentier, Nicolas Verzelen. Optimal multiple change-point detection for high-dimensional data. 2022. hal-03004860v3

HAL Id: hal-03004860

https://hal.science/hal-03004860v3

Preprint submitted on 7 Dec 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Optimal multiple change-point detection for high-dimensional data

Emmanuel Pilliat, Alexandra Carpentier, Nicolas Verzelen

► To cite this version:

Emmanuel Pilliat, Alexandra Carpentier, Nicolas Verzelen. Optimal multiple change-point detection for high-dimensional data. Electronic Journal of Statistics, 2023, 17 (1), pp.1240-1315. 10.1214/23-EJS2126. hal-03004860v4

HAL Id: hal-03004860

https://hal.science/hal-03004860v4

Submitted on 6 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Optimal multiple change-point detection for high-dimensional data

Emmanuel Pilliat¹, Alexandra Carpentier², and Nicolas Verzelen³

¹Université de Montpellier, Montpellier, France

²Institut für Mathematik, Universität Potsdam, Germany

³INRAE Montpellier, Montpellier, France

December 7, 2022

## Abstract

This manuscript makes two contributions to the field of change-point detection. In a general change-point setting, we provide a generic algorithm for aggregating local homogeneity tests into an estimator of change-points in a time series. Interestingly, we establish that the error rates of the collection of tests directly translate into detection properties of the change-point estimator. This generic scheme is then applied to various problems including covariance change-point detection, nonparametric change-point detection and sparse multivariate mean change-point detection. For the latter, we derive minimax optimal rates that are adaptive to the unknown sparsity and to the distance between change-points when the noise is Gaussian. For sub-Gaussian noise, we introduce a variant that is optimal in almost all sparsity regimes.

## 1 Introduction

Change-point detection has a long history since the seminal work of Wald [39] that lead to flourishing lines (see [31, 36] for recent surveys). Earlier contributions focused on the problems of detecting and localizing change-points in a univariate time series. Spurred by applications in genomics [32] and finance, there has been a recent trend in the literature towards the analysis of more complex time series for instance in a high-dimensional linear space [21] or even belonging to a non-Euclidean space [8].

In this work, we study high-dimensional time series whose mean may change possibly on a few number of coordinates. See the introduction of [46] for an account of possible applications and practical motivations. In particular, we build a procedure which is able to detect and localize change-points under minimal assumptions on the height of these change-points. Along the way towards this optimal procedure, we define and analyze a scheme for general change-point problems that aggregates a collection of local tests into an estimator change-points. This generic scheme is of independent interest and easily allows to derive optimal change-point procedure in other complex settings such as covariance change-points problems or nonparametric change-point problems. In this introduction, we first describe this generic scheme before turning to our results in high-dimensional sparse change-point detection and finally discussing other applications.

### 1.1 General change-point setting

In the most general form of a change-point problem, we consider a random sequence $Y = (y_1, y_2, \dots, y_n)$ in some measured space $\mathcal{Y}^n$ and, for $t = 1, \dots, n$, we write $\mathbb{P}_t$ for the marginal
# Optimal multiple change-point detection for high-dimensional data*

Emmanuel Pilliat¹, Alexandra Carpentier² and Nicolas Verzelen³

¹IMAG, Univ Montpellier, CNRS, Montpellier, France.
e-mail: emmanuel.pilliat@umontpellier.fr

²Institut für Mathematik, Universität Potsdam, Germany.
e-mail: carpentier@uni-potsdam.de

³INRAE, Mistea, Institut Agro, Univ Montpellier, Montpellier, France.
e-mail: nicolas.verzelen@inrae.fr

**Abstract:** This manuscript makes two contributions to the field of change-point detection. In a general change-point setting, we provide a generic algorithm for aggregating local homogeneity tests into an estimator of change-points in a time series. Interestingly, we establish that the error rates of the collection of tests directly translate into detection properties of the change-point estimator. This generic scheme is then applied to various problems including covariance change-point detection, nonparametric change-point detection and sparse multivariate mean change-point detection. For the latter, we derive minimax optimal rates that are adaptive to the unknown sparsity and to the distance between change-points when the noise is Gaussian. For sub-Gaussian noise, we introduce a variant that is optimal in almost all sparsity regimes.

**Keywords and phrases:** Multivariate time series, minimax rate, Gaussian noise.

Received February 2022.

## 1. Introduction

Change-point detection has a long history since the seminal work of Wald [39] that lead to flourishing lines (see [31, 36] for recent surveys). Earlier contributions focused on the problems of detecting and localizing change-points in a univariate time series. Spurred by applications in genomics [32] and finance, there has been a recent trend in the literature towards the analysis of more

<footnote>*The work of A. Carpentier is partially supported by the Deutsche Forschungsgemeinschaft (DFG) Emmy Noether grant MuSyAD (CA 1488/1-1), by the DFG - 314838170, GRK 2297 MathCoRe, by the FG DFG, by the DFG CRC 1294 ‘Data Assimilation’, Project A03, by the Forschungsgruppe FOR 5381 “Mathematical Statistics in the Information Age – Statistical Efficiency and Computational Tractability”, Project TP 02, by the Agence Nationale de la Recherche (ANR) and the DFG on the French-German PRCI ANR ASCAI CA 1488/4-1 “Aktive und Batch-Segmentierung, Clustering und Seriation: Grundlagen der KI” and by the UFA-DFH through the French-German Doktorandenkolleg CDFA 01-18 and by the SFI Sachsen-Anhalt for the project RE-BCI. The work of E. Pilliat and N. Verzelen has been partially supported by ANR-21-CE23-0035 (ASCAI).</footnote>
distribution of $y_t$. We are also given a functional $\Gamma$ mapping the probability distribution $P_t$ to some space $\mathcal{V}$. Then, the purpose of change-point detection is to detect changes in the sequence $(\Gamma(\mathbb{P}_1), \Gamma(\mathbb{P}_2), \dots, \Gamma(\mathbb{P}_n))$ in $\mathcal{V}^n$ and to estimate the positions of these changes. This setting is really general and does not require that the random variables $(y_t)$ are independent.

Let us shortly explain how this general framework encompasses most offline change-point detection problems. In the Gaussian mean univariate change-point setting, we have $\mathcal{Y} = \mathbb{R}$, the distribution $P_t$ corresponds to the normal distribution with mean $\theta_t \in \mathbb{R}$ and variance $\sigma^2$ and $\Gamma(\mathbb{P}_t) = \theta_t$. In the (heteroscedastic) mean univariate change-point problem, the distribution $P_t$ is not necessarily Gaussian and, in particular, the variance of $y_t$ is allowed to vary with $t$. Still, one is only interested in detecting variations of $\Gamma(\mathbb{P}_t) = \int x d\mathbb{P}_t = \mathbb{E}[y_t]$. By contrast, in the *variance* univariate change-point problems, one focuses on changes in the variance of $y_t$. This can be done by taking $\Gamma(\mathbb{P}_t) = \int x^2 d\mathbb{P}_t - [\int x d\mathbb{P}_t]^2 = \text{Var}(y_t)$. If one is interested in possibly nonparametric changes in the distributions, then the functional $\Gamma$ is simply taken to be the identity map. In semi-parametric quantile change-point detection [22], the univariate distributions $P_t$ can be arbitrary whereas $\Gamma(\mathbb{P}_t)$ is a quantile of $P_t$.

To further formalize the change-point detection problem in the sequence $(\Gamma(\mathbb{P}_1), \Gamma(\mathbb{P}_2), \dots, \Gamma(\mathbb{P}_n))$, we define an integer $0 \le K \le n-1$ and a vector of integers $\tau = (\tau_1, \dots, \tau_K)$ satisfying $1 = \tau_0 < \tau_1 < \dots < \tau_K < \tau_{K+1} = n+1$ such that $\Gamma(\mathbb{P}_t)$ is constant over each interval $[\tau_k, \tau_{k+1} - 1]$ and $\Gamma(\mathbb{P}_{\tau_k - 1}) \ne \Gamma(\mathbb{P}_{\tau_k})$. Hence, $\tau_k$ corresponds to the *position* of the $k^{th}$ change-point. We shall often refer to $\tau_k$ as a *change-point*. Equipped with this notation, we are interested in building an estimator $\hat{\tau} = (\hat{\tau}_1, \dots, \hat{\tau}_{\hat{K}})$ of $\tau$ from the time series $Y$. Here, $\hat{\tau}_1, \dots, \hat{\tau}_{\hat{K}}$ correspond to the *estimated change-points* of $\tau$ and $\hat{K}$ to the number of the estimated change-points.

### 1.1.1 Desirable Guarantees of an estimator.

Before describing the generic scheme for estimating $\tau$, let us first formalize the desired properties of a good change-point procedure. Informally, the primary objectives are to detect most if not all change-points while estimating no (or at least very few) spurious change-points.

Regarding the latter objective, it is usually required that the number of change-points $K$ is not overestimated by $\hat{\tau}$. Here, we require a slightly stronger local property introduced in [38]. An estimator $\hat{\tau}$ of size $\hat{K}$ is said to detect no spurious change-points (NoSp) if

$$ \left\{ \left| \begin{array}{l} \{\hat{\tau}_{k'}, 1 \le k' \le \hat{K}\} \cap \left[ \tau_k - \frac{\tau_k - \tau_{k-1}}{2}, \tau_k + \frac{\tau_{k+1} - \tau_k}{2} \right] \\ \{\hat{\tau}_{k'}, 1 \le k' \le \hat{K}\} \subset \left[ \tau_1 - \frac{\tau_1 - 1}{2}, \tau_K + \frac{n+1-\tau_K}{2} \right] \end{array} \right| \le 1, \quad \text{for all } 1 \le k \le K; \right. \tag{1} $$

The second condition simply ensures that no change-point is estimated near the boundaries of the time series. The first condition entails that, for each change-point $\tau_k$ there is at most one estimated change-point $\hat{\tau}_k$ in the interval $[\tau_k - (\tau_k - \tau_{k-1})/2, \tau_k + (\tau_{k+1} - \tau_k)/2]$. In other words, (NoSp) requires that, on each sub-interval, the number of change-points is not overestimated.

Let us now formalize the objective of detecting the change-points. In this work, we consider as in [38] realistic settings where some change-points are so close or their heights are so small that they are impossible to detect. As a consequence, we can only hope to detect the subset of significant change-points. In what follows, we define a subset $\mathcal{K}^* \subset [K]$ of change-point indices that correspond to significant change-points. Obviously, the significance of a particular change-point is relative to the problem under consideration - data distribution, nature of change-points - and the definition is problem dependent. As an example, we define in the next subsection the suitable notion of energy
complex time series for instance in a high-dimensional linear space [21] or even belonging to a non-Euclidean space [8].

In this work, we study high-dimensional time series whose mean may change possibly on a few number of coordinates. See the introduction of [46] for an account of possible applications and practical motivations. In particular, we build a procedure which is able to detect and localize change-points under minimal assumptions on the height of these change-points. Along the way towards this optimal procedure, we define and analyze a scheme for general change-point problems that aggregates a collection of local tests into an estimator change-points. This generic scheme is of independent interest and easily allows to derive optimal change-point procedure in other complex settings such as covariance change-points problems or nonparametric change-point problems. In this introduction, we first describe this generic scheme before turning to our results in high-dimensional sparse change-point detection and finally discussing other applications.

## 1.1. General change-point setting

In the most general form of a change-point problem, we consider a random sequence $Y = (y_1, y_2, \dots, y_n)$ in some measured space $\mathcal{Y}^n$ and, for $t = 1, \dots, n$, we write $\mathbb{P}_t$ for the marginal distribution of $y_t$. We are also given a functional $\Gamma$ mapping the probability distribution $\mathbb{P}_t$ to some space $\mathcal{V}$.

Then, the purpose of change-point detection is to detect changes in the sequence $(\Gamma(\mathbb{P}_1), \Gamma(\mathbb{P}_2), \dots, \Gamma(\mathbb{P}_n))$ in $\mathcal{V}^n$ and to estimate the positions of these changes. This setting is really general and does not require that the random variables $(y_t)$ are independent.

Let us shortly explain how this general framework encompasses most offline change-point detection problems. In the Gaussian mean univariate change-point setting, we have $\mathcal{Y} = \mathbb{R}$, the distribution $\mathbb{P}_t$ corresponds to the normal distribution with mean $\theta_t \in \mathbb{R}$ and variance $\sigma^2$ and $\Gamma(\mathbb{P}_t) = \theta_t$. In the (heteroscedastic) mean univariate change-point problem, the distribution $\mathbb{P}_t$ is not necessarily Gaussian and, in particular, the variance of $y_t$ is allowed to vary with $t$. Still, one is only interested in detecting variations of $\Gamma(\mathbb{P}_t) = \int x d\mathbb{P}_t = \mathbb{E}[y_t]$. By contrast, in the *variance* univariate change-point problems, one focuses on changes in the variance of $y_t$. This can be done by taking $\Gamma(\mathbb{P}_t) = \int x^2 d\mathbb{P}_t - [\int x d\mathbb{P}_t]^2 = \text{Var}(y_t)$. If one is interested in possibly nonparametric changes in the distributions, then the functional $\Gamma$ is simply taken to be the identity map. In semi-parametric quantile change-point detection [22], the univariate distributions $\mathbb{P}_t$ can be arbitrary whereas $\Gamma(\mathbb{P}_t)$ is a quantile of $\mathbb{P}_t$.

To further formalize the change-point detection problem in the sequence $(\Gamma(\mathbb{P}_1), \Gamma(\mathbb{P}_2), \dots, \Gamma(\mathbb{P}_n))$, we define an integer $0 \le K \le n-1$ and a vector of integers $\tau = (\tau_1, \dots, \tau_K)$ satisfying $1 = \tau_0 < \tau_1 < \dots < \tau_K < \tau_{K+1} = n+1$ such that $\Gamma(\mathbb{P}_t)$ is constant over each interval $[\tau_k, \tau_{k+1} - 1]$ and $\Gamma(\mathbb{P}_{\tau_k - 1}) \ne \Gamma(\mathbb{P}_{\tau_k})$. Hence, $\tau_k$ corresponds to the *position* of the $k^{th}$ change-point. We shall often refer to $\tau_k$ as a *change-point*. Equipped with this notation, we are interested
and significance of a change-point in the mean multivariate change-point setting. In Section 6, we formalize this notion for covariance and univariate nonparametric change-point problems. In light of this discussion, the second guarantee we aim for is to detect all significant change-points. A change-point $\tau_k$ is said to be detected if there is at least one estimated change-point $\hat{\tau}_l$ in the interval $[\tau_k - (\tau_k - \tau_{k-1})/2, \tau_k + (\tau_{k+1} - \tau_k)/2]$. Equivalently, this means that at least one of the estimated change-points is closer to $\tau_k$ than to any other true change-point.

Aside from (NoSp) and (detect) properties, one may additionally aim at localizing the change-points as well as possible – see the discussions in [41]. Given a specific change-point $\tau_k$ detected by an estimator $\hat{\tau}$, its localization error $d_{H,1}(\hat{\tau}, \tau_k)$ is defined by

$$d_{H,1}(\hat{\tau}, \tau_k) = \min_{l=1,...,|\hat{\tau}|} |\hat{\tau}_l - \tau_k|,$$

which is the smallest distance between $\tau_k$ and one of the estimated change-points. While this work mainly focused on the detection problem, we shall also provide localization bounds along the way.

### 1.1.2 A generic roadmap for change-point detection.

In this manuscript, our first contribution is a generic procedure for aggregating a collection of tests into an estimator $\hat{\tau}$ of $\tau$. For two positive integers $(l, r)$, we consider the time interval $[l-r, l+r)$. Suppose we are given a collection $\mathcal{G}$ of such $(l, r)$. For each $(l, r) \in \mathcal{G}$, we are also given a homogeneity test $T_{l,r}$ of the null hypothesis $\mathcal{H}_0: \{(\Gamma(\mathbb{P}_t))\}$ is constant over the segment $[l-r, l+r)$. This hypothesis is equivalent to the absence of any change-point on the interval $(l-r, l+r)$. Given such a collections of homogeneity tests $(T_{l,r})$, $(l, r) \in \mathcal{G}$, we build in this manuscript an estimator $\hat{\tau}$ that satisfies the following properties. If the multiple testing procedure does not reject any true null hypothesis (no false positives), then $\hat{\tau}$ does not estimate any spurious change-point, that is, it satisfies (No Sp). Furthermore, any change-point $\tau_k$ that is detected by some test $T_{\bar{\tau}_k, \bar{r}_k}$, where $\bar{\tau}_k$ is close enough to $\tau_k$ and $\bar{r}_k$ is small enough is detected by the estimator $\hat{\tau}$. In other words, we establish a completely generic result that translates properties of the multiple testing procedure into detection properties. Thus, the construction of a change-point procedure boils down to building a suitable multiple testing procedure $(T_{l,r})$, $(l, r) \in \mathcal{G}$ whose family-wise error rate (FWER) is controlled, while being able to detect all the significant change-points. In turn, this allows us to reduce the problem of change-point detection under minimal distance between the change-points to the well-established field of minimax testing.

### 1.1.3 Related Work and possible applications.

In the last years, there has been a growing interest into the extension of univariate mean change-point procedures such as wild binary segmentation (WBS) [14] to other problems such as covariance change-point [40], network change-point [41], or nonparametric change-point [33]. For each of these problems (and for others), it turns out that the general ideas of WBS can be instantiated. However, for each setting, the proofs need to be fully adapted in a case by case manner. Besides, the resulting procedures are only optimal up to logarithmic terms.

Recently, Chan and Chen [5] and Kovács et al. [24] have introduced bottom-up aggregation procedures for mean change-point segmentation (see also [25] for localization improvements). Moreover, Kovács et al. [24, 25] illustrate the numerical performances to other change-point models, such as graphical models or multivariate mean-change point models. In fact, one may extend their procedures to generic problems, but the theoretical guarantees are only provided for univariate models and it remains unclear whether one can extend them beyond very specific cases.
in building an estimator $\hat{\tau} = (\hat{\tau}_1, \dots, \hat{\tau}_K)$ of $\tau$ from the time series $Y$. Here, $\hat{\tau}_1, \dots, \hat{\tau}_K$ correspond to the estimated change-points of $\tau$ and $K$ to the number of the estimated change-points.

### 1.1.1. Desirable guarantees of an estimator

Before describing the generic scheme for estimating $\tau$, let us first formalize the desired properties of a good change-point procedure. Informally, the primary objectives are to detect most if not all change-points while estimating no (or at least very few) spurious change-points.

Regarding the latter objective, it is usually required that the number of change-points $K$ is not overestimated by $\hat{\tau}$. Here, we require a slightly stronger local property introduced in [38]. An estimator $\hat{\tau}$ of size $K$ is said to detect no spurious change-points (NoSp) if

$$ \left\{ \begin{array}{l} \left| \{ \hat{\tau}_{k'}, 1 \le k' \le K \} \cap \left[ \tau_k - \frac{\tau_{k-1} - \tau_{k'}}{2}, \tau_k + \frac{\tau_{k+1} - \tau_k}{2} \right] \right| \le 1 \quad \text{for all } 1 \le k \le K; \\ \{ \hat{\tau}_{k'}, 1 \le k' \le K \} \subset \left[ \tau_1 - \frac{\tau_1 - 1}{2}, \tau_K + \frac{n+1-\tau_K}{2} \right]. \end{array} \right. \qquad (1) $$

The second condition simply ensures that no change-point is estimated near the boundaries of the time series. The first condition entails that, for each change-point $\tau_k$ there is at most one estimated change-point $\hat{\tau}_k$ in the interval $[\tau_k - (\tau_k - \tau_{k-1})/2, \tau_k + (\tau_{k+1} - \tau_k)/2]$. In other words, (NoSp) requires that, on each sub-interval, the number of change-points is not overestimated.

Let us now formalize the objective of detecting the change-points. In this work, we consider as in [38] realistic settings where some change-points are so close or their heights are so small that they are impossible to detect. As a consequence, we can only hope to detect the subset of significant change-points. In what follows, we define a subset $K^* \subset [K]$ of change-point indices that correspond to significant change-points. Obviously, the significance of a particular change-point is relative to the problem under consideration – data distribution, nature of change-points – and the definition is problem dependent. As an example, we define in the next subsection the suitable notion of energy and significance of a change-point in the mean multivariate change-point setting. In Section 6, we formalize this notion for covariance and univariate nonparametric change-point problems. In light of this discussion, the second guarantee we aim for is to detect all significant change-points. A change-point $\tau_k$ is said to be detected if there is at least one estimated change-point $\hat{\tau}_l$ in the interval $[\tau_k - (\tau_k - \tau_{k-1})/2, \tau_k + (\tau_{k+1} - \tau_k)/2]$. Equivalently, this means that at least one of the estimated change-points is closer to $\tau_k$ than to any other true change-point.

Aside from (NoSp) and (detect) properties, one may additionally aim at localizing the change-points as well as possible – see the discussions in [41]. Given a specific change-point $\tau_k$ detected by an estimator $\hat{\tau}$, its localization error $d_{H,1}(\hat{\tau}, \tau_k)$ is defined by

$$ d_{H,1}(\hat{\tau}, \tau_k) = \min_{l=1,\dots,|\hat{\tau}|} |\hat{\tau}_l - \tau_k|, $$
In contrast, it is quite straightforward to adapt our generic procedure to any new setting once suitable homogeneity multiple tests have been crafted. As the most prominent example, we consider the sparse high-dimensional mean change-point detection and establish the optimality of our procedure – see the next subsection for details. In Section 6, we also handle the covariance change-point detection and the univariate nonparametric change-point detection problems. In each case, we pinpoint the first tight minimal conditions for detection.

Besides, we could apply our strategy to other problems such changes in auto-regressive models [43], changes in the inverse covariance matrix of $y_i$ [17, 24] or changes in a high-dimensional regression model [34]. All such change-point problems can be addressed through the construction and careful analysis of two-sample tests for auto-regressive models, inverse covariance matrices, and linear regression models respectively. Similarly, we can build Kernel change-point procedures [1, 16] from kernel two-sample tests [18].

## 1.2 Sparse Multivariate Change-point Setting

As explained above, our primary application of our generic scheme is the multivariate mean change-point detection problem with sparse variations where one observes a time series $Y = (y_1, \dots, y_n) \in \mathbb{R}^{p \times n}$ with unknown means $\Theta = (\theta_1, \dots, \theta_n) \in \mathbb{R}^{p \times n}$ so that we have the decomposition

$$y_t = \theta_t + \varepsilon_t \quad t = 1, \dots, n, \qquad (2)$$

where the noise matrix $\varepsilon = (\varepsilon_1, \dots, \varepsilon_n)$ is made of independent and mean zero random vectors of size $p$. In this manuscript, we make two distributional assumptions on the noise. Either we suppose that all random vectors $\varepsilon_i$ follow independent normal distribution with variance $\sigma^2 \mathbf{I}_p$ (see Section 3) or that the components of $\varepsilon_i$ follow independent sub-Gaussian distributions with variance $\sigma^2$ (see Section 4). In either case, we assume that $\sigma^2$ is known.

Here, we are interested in the variations of the mean vector $\theta_t$ so that, relying on the formalism of the previous subsection, we have $\Gamma(\mathbb{P}_t) = \theta_t$. Considering the vector of change-points $\tau = (\tau_1, \dots, \tau_K)$, we can define $K + 1$ vectors $\mu_0, \dots, \mu_K$ in $\mathbb{R}^p$ satisfying $\mu_k \neq \mu_{k+1}$ for all $k = 0, \dots, K - 1$ such that

$$\theta_t = \sum_{k=0}^{K} \mu_k \mathbf{1}_{\tau_k \le t < \tau_{k+1}} .$$

Equivalently, $\mu_k$ is the constant mean of $y$ over the interval $[\tau_k, \tau_{k+1} - 1]$. The difference $\mu_k - \mu_{k-1}$ in $\mathbb{R}^p$ measures the variation of $\Theta$ at the change-point $\tau_k$ and can possibly have many null coordinates. In this possibly sparse multi-dimensional setting, the significance of a change-point is measured through three quantities $\Delta_k$, $r_k$, and $s_k$. First, the height $\Delta_k$ of the change-point $\tau_k$ is defined as the Euclidean norm of the signal difference. The length $r_k$ of the change-point $\tau_k$ is the minimal distance from $\tau_k$ to another change-point, $\tau_{k-1}$ or $\tau_{k+1}$. More precisely,

$$\Delta_k = \|\mu_k - \mu_{k-1}\|; \quad r_k = \min(\tau_{k+1} - \tau_k, \tau_k - \tau_{k-1}). \qquad (3)$$

As a simple example, Figure 1 depicts a one dimensional piece-wise constant sequence $\Theta$ with 3 change-points illustrating the setting presented above. In the univariate change-point literature (e.g. [7, 14, 15]) the height and the length of a change-point characterize the significance of a change-point. In the multivariate setting, where the change-points can be sparse, meaning the number of non null coordinates of the vector $\mu_k - \mu_{k-1}$ is possibly small, one also considers the sparsity $s_k$ of change-point $\tau_k$, defined by

$$s_k = \|\mu_k - \mu_{k-1}\|_0, \qquad (4)$$

where, for any $v \in \mathbb{R}^p$, $\|v\|_0 = \sum_{1 \le i \le p} \mathbf{1}\{v_i \ne 0\}$.
which is the smallest distance between $\tau_k$ and one of the estimated change-points. While this work mainly focused on the detection problem, we shall also provide localization bounds along the way.

### 1.1.2. A generic roadmap for change-point detection

In this manuscript, our first contribution is a generic procedure for aggregating a collection of tests into an estimator $\hat{\tau}$ of $\tau$. For two positive integers $(l, r)$, we consider the time interval $[l - r, l + r)$. Suppose we are given a collection $\mathcal{G}$ of such $(l, r)$. For each $(l, r) \in \mathcal{G}$, we are also given a homogeneity test $T_{l,r}$ of the null hypothesis $H_0: \{(\Gamma(\mathbb{P}_t))\}$ is constant over the segment $[l - r, l + r)$. This hypothesis is equivalent to the absence of any change-point on the interval $(l - r, l + r)$. Given such a collection of homogeneity tests $(T_{l,r}), (l,r) \in \mathcal{G}$, we build in this manuscript an estimator $\hat{\tau}$ that satisfies the following properties. If the multiple testing procedure does not reject any true null hypothesis (no false positives), then $\hat{\tau}$ does not estimate any spurious change-point, that is, it satisfies (No Sp). Furthermore, any change-point $\tau_k$ that is detected by some test $T_{\tau_k, \bar{\tau}_k}$, where $\bar{\tau}_k$ is close enough to $\tau_k$ and $\bar{\tau}_k$ is small enough is detected by the estimator $\hat{\tau}$. In other words, we establish a completely generic result that translates properties of the multiple testing procedure into detection properties. Thus, the construction of a change-point procedure boils down to building a suitable multiple testing procedure $(T_{l,r}), (l,r) \in \mathcal{G}$ whose family-wise error rate (FWER) is controlled, while being able to detect all the significant change-points. In turn, this allows us to reduce the problem of change-point detection under minimal distance between the change-points to the well-established field of minimax testing.

### 1.1.3. Related Work and possible applications

In the last years, there has been a growing interest into the extension of univariate mean change-point procedures such as wild binary segmentation (WBS) [14] to other problems such as covariance change-point [40], network change-point [41], or nonparametric change-point [33]. For each of these problems (and for others), it turns out that the general ideas of WBS can be instantiated. However, for each setting, the proofs need to be fully adapted in a case by case manner. Besides, the resulting procedures are only optimal up to logarithmic terms.

Recently, Chan and Chen [5] and Kovács et al. [24] have introduced bottom-up aggregation procedures for mean change-point segmentation (see also [25] for localization improvements). Moreover, Kovács et al. [24, 25] illustrate the numerical performances to other change-point models, such as graphical models or multivariate mean-change point models. In fact, one may extend their procedures to generic problems, but the theoretical guarantees are only provided for univariate models and it remains unclear whether one can extend them beyond very specific cases.
<figure><img src="image_2.png" /><figcaption>Figure 1: An example of a piece-wise constant sequence $\Theta$ with 3 change-points and $p = 1$.</figcaption></figure>

### 1.2.1 Two-sample tests and CUSUM statistics

Our objective is to detect and recover positions $(\tau_k)_{k \le K}$ under minimal conditions on the change-point height $\Delta_k$, change-point length $r_k$ and sparsity $s_k$. In view of the generic change-point procedure discussed in the previous subsection, this mainly boils down to building suitable tests of the assumptions $\{\Theta$ is constant over $[l - r, l + r)\}$ versus $\{\Theta$ is not constant on this segment$\}$. Following the literature on binary and wild binary segmentation, we consider the CUSUM statistic

$$C_{l,r}(Y) = \sqrt{\frac{r}{2\sigma^2}} \left( \frac{1}{r} \sum_{i=l}^{l+r-1} y_i - \frac{1}{r} \sum_{i=l-r}^{l-1} y_i \right).$$

This statistic computes the normalized difference of empirical mean of $y_i$ on $[l - r, l)$ and $[l, l + r)$. If the noise is Gaussian and if $\Theta$ is constant on $[l - r, l + r)$, then $C_{l,r}(Y)$ simply follows a standard $p$-dimensional normal distribution. To simplify, consider a specific instance of our testing problem where we want to test whether $\{\Theta$ is constant over $[l - r, l + r)\}$ versus $\{\Theta$ contains exactly one change-point at $l$ on the segment $[l - r, l + r)\}$. This corresponds to a two-sample mean testing problem, for which the CUSUM statistic $C_{l,r}(Y)$ is a sufficient statistic if the noise is Gaussian. Then, given $C_{l,r}(Y)$, one wants to test whether its expectation is 0 (no change-point on $[l - r, l + r)$) versus its expectation is non-zero but is $s$-sparse for some unknown $s$. This classical detection problem is well understood [11] and it is well known that a combination of a $\chi^2$-type test with a higher-criticism-type test is optimal. Here, the challenge stems from the fact that we do not want to perform a single such test, but a large collection of tests over a collection of $(l, r) \in \mathcal{G}$.

### 1.2.2 Our contribution

As usual in the mean change-point literature, we consider the energy $r_k \Delta_k^2$ of the change-point $\tau_k$. Up to a possible factor in $[1/2, 1]$, $r_k \Delta_k^2$ is the square distance between $\Theta$ and its projection on the space of vectors $\Theta'$ with change-point at $(\tau_1, \dots, \tau_{k-1}, \tau_{k+1}, \dots, \tau_K)$ –see e.g. [38] for a discussion in the univariate setting. In other words, the energy $r_k \Delta_k^2$ characterizes the significance
In contrast, it is quite straightforward to adapt our generic procedure to any new setting once suitable homogeneity multiple tests have been crafted. As the most prominent example, we consider the sparse high-dimensional mean change-point detection and establish the optimality of our procedure – see the next subsection for details. In Section 6, we also handle the covariance change-point detection and the univariate nonparametric change-point detection problems. In each case, we pinpoint the first tight minimal conditions for detection.

Besides, we could apply our strategy to other problems such changes in auto-regressive models [43], changes in the inverse covariance matrix of $y_i$ [17, 24] or changes in a high-dimensional regression model [34]. All such change-point problems can be addressed through the construction and careful analysis of two-sample tests for auto-regressive models, inverse covariance matrices, and linear regression models respectively. Similarly, we can build Kernel change-point procedures [1, 16] from kernel two-sample tests [18].

## 1.2. Sparse multivariate change-point setting

As explained above, our primary application of our generic scheme is the multivariate mean change-point detection problem with sparse variations where one observes a time series $Y = (y_1, \dots, y_n) \in \mathbb{R}^{p \times n}$ with unknown means $\Theta = (\theta_1, \dots, \theta_n) \in \mathbb{R}^{p \times n}$ so that we have the decomposition

$$y_t = \theta_t + \varepsilon_t \quad t = 1, \dots, n, \qquad (2)$$

where the noise matrix $\varepsilon = (\varepsilon_1, \dots, \varepsilon_n)$ is made of independent and mean zero random vectors of size $p$. In this manuscript, we make two distributional assumptions on the noise. Either we suppose that all random vectors $\varepsilon_i$ follow independent normal distribution with variance $\sigma^2 \mathbf{I}_p$ (see Section 3) or that the components of $\varepsilon_i$ follow independent sub-Gaussian distributions with variance $\sigma^2$ (see Section 4). In either case, we assume that $\sigma^2$ is known.

Here, we are interested in the variations of the mean vector $\theta_t$ so that, relying on the formalism of the previous subsection, we have $\Gamma(\mathbb{P}_t) = \theta_t$. Considering the vector of change-points $\tau = (\tau_1, \dots, \tau_K)$, we can define $K+1$ vectors $\mu_0, \dots, \mu_K$ in $\mathbb{R}^p$ satisfying $\mu_k \neq \mu_{k+1}$ for all $k = 0, \dots, K-1$ such that

$$\theta_t = \sum_{k=0}^{K} \mu_k \mathbf{1}_{\tau_k \le t < \tau_{k+1}} .$$

Equivalently, $\mu_k$ is the constant mean of $y$ over the interval $[\tau_k, \tau_{k+1} - 1]$. The difference $\mu_k - \mu_{k-1}$ in $\mathbb{R}^p$ measures the variation of $\Theta$ at the change-point $\tau_k$ and can possibly have many null coordinates. In this possibly sparse multi-dimensional setting, the significance of a change-point is measured through three quantities $\Delta_k, r_k$, and $s_k$. First, the height $\Delta_k$ of the change-point $\tau_k$ is defined as the Euclidean norm of the signal difference. The length $r_k$ of the change-point $\tau_k$ is the minimal distance from $\tau_k$ to another change-point, $\tau_{k-1}$ or $\tau_{k+1}$. More
of the change-point $\tau_k$. In Section 3, we introduce a multi-scale change-point detection procedure detecting any change-point $\tau_k$ whose energy is higher, up to a numerical constant, than $\sigma^2 s_k \log(1+\frac{\sqrt{p}}{s_k}\sqrt{\log(n/r_k)}) + \sigma^2 \log(n/r_k)$. This result is valid for arbitrary length $r_k$ and sparsity $s_k$, and does not require the knowledge of these two quantities. In summary, our procedure does not estimate any spurious change-point (**NoSp**) and **detects** all the change-points whose energy are higher than the latter threshold. In Section 5, we establish that, as soon as the unknown number $K$ of the change-points is larger than 1, the condition $\sigma^2 s_k \log(1 + \frac{\sqrt{p}}{s_k}\sqrt{\log(n/r_k)}) + \sigma^2 \log(n/r_k)$ on the energy is tight with respect to $n, p, r_k$ and $s_k$, in the sense that no procedure achieving (**NoSp**) is able to detect with high probability a change-point whose energy is smaller (up to some constant) than the latter threshold. In Section 4, we consider the more general setting where the noise is $L$-sub-gaussian with known variance, and we establish a similar result to the Gaussian case up to a logarithmic loss in some regimes. Finally, we illustrate in Section 8 the behavior of our procedure on numerical experiments.

### 1.2.3 Related work

For dense change-points ($s_k = p$) but with unknown covariance for the noise, Wang et al. [45] (see also [44]) study the behavior of a procedure based on $U$-statistics of the CUSUM. Jirak [21] and Yu and Chen [48] introduce binary segmentation procedures based on the $l_\infty$ norm of the CUSUMs. Although those work explicitly characterize the asymptotic distribution of the test statistics and, for some of them, allow temporal dependencies in the data, the corresponding energy requirements for change-point detection are either not studied or turn out to be suboptimal.

Closest to our work, Chan and Chen [5] study a bottom-up approach to detect change-points of a Gaussian multivariate time series in an asymptotic setting. More specifically, the authors consider an asymptotic regime where the size of the time series is exponential in the dimension: $n = e^{p\zeta}$ with $\zeta \in (0, 1)$. The authors also assume that the number $K$ of change-points remains finite when $n, p \to \infty$ and that the minimal sparsity $s$ of these change-points is polynomial is $p$. In this specific regime, their procedures provably recover change-points under a near-minimal (up to logarithmic factors with respect to $n$) condition on the energy. In contrast, our results provide non-asymptotic and tight results for all scaling with respect to $n$ and $p$, allow for arbitrarily large number $K$ of change-points and allow for the presence of non-significant change-points. In the same specific asymptotic setting, [20] introduce a so called score test statistic used in a change-point detection procedure which is shown to achieve the same performance as [5] in the gaussian model but also handle Poisson observations.

Recently, Liu et al. [28] have characterized the optimal detection rate of a possibly sparse change-point in the specific case where there is at most one change-point, but the optimal rates are significantly slower in the multiple change-point setting. See also [12] and [9] for earlier results. Wang and Samworth [46] have proposed the INSPECT method based on sparse projection to handle sparse change-points, but INSPECT provably detects the change-points under strong assumption on the energy; see Section 3 for a precise comparison.

In the univariate setting ($p = 1$), minimal energy requirements for change-point detection are well understood [13, 15, 38, 42] and are nearly achieved by a wide range of procedures including penalized least-square and multi-scale tests methods.
<figure><img src="image_2.png" /><figcaption>FIG 1. An example of a piece-wise constant sequence $\Theta$ with 3 change-points and $p = 1$.</figcaption></figure>

precisely,

$$ \Delta_k = \|\mu_k - \mu_{k-1}\|; \quad r_k = \min(\tau_{k+1} - \tau_k, \tau_k - \tau_{k-1}). \quad (3) $$

As a simple example, Figure 1 depicts a one dimensional piece-wise constant sequence $\Theta$ with 3 change-points illustrating the setting presented above. In the univariate change-point literature (e.g. [7, 14, 15]) the height and the length of a change-point characterize the significance of a change-point. In the multivariate setting, where the change-points can be sparse, meaning the number of non null coordinates of the vector $\mu_k - \mu_{k-1}$ is possibly small, one also considers the sparsity $s_k$ of change-point $\tau_k$, defined by

$$ s_k = \|\mu_k - \mu_{k-1}\|_0, \quad (4) $$

where, for any $v \in \mathbb{R}^p$, $\|v\|_0 = \sum_{1 \le i \le p} 1\{v_i \ne 0\}$.

### 1.2.1. Two-sample tests and CUSUM statistics

Our objective is to detect and recover positions $(\tau_k)_{k \le K}$ under minimal conditions on the change-point height $\Delta_k$, change-point length $r_k$ and sparsity $s_k$. In view of the generic change-point procedure discussed in the previous subsection, this mainly boils down to building suitable tests of the assumptions $\{\Theta$ is constant over $[l - r, l + r)\}$ versus $\{\Theta$ is not constant on this segment\}$. Following the literature on binary and wild binary segmentation, we consider the CUSUM statistic

$$ C_{l,r}(Y) = \sqrt{\frac{r}{2\sigma^2}} \left( \frac{1}{r} \sum_{i=l}^{l+r-1} y_i - \frac{1}{r} \sum_{i=l-r}^{l-1} y_i \right). $$
## 2 A Generic algorithm for multiscale change-point detection on a grid

In this section, we study the problem of change-point detection in the general setting defined in Section 1.1. We introduce a bottom-up algorithm that aggregates a collection of homogeneity tests, performed at many positions, and for many scales, of our data. Then, we establish that, under some conditions on these tests, the procedure detects significant change-points.

### 2.1 Grid and multiscale statistics

Since our purpose is to translate a collection of local tests $T = (T_{l,r})_{(l,r) \in \mathcal{G}}$ indexed by a grid $\mathcal{G}$ into a change-point detection procedure, we first need to formalize what we mean by a grid. Henceforth, we call a grid $\mathcal{G}$ of $[n]$ a collection of locations and scales where a scale $r$ is a positive integer smaller or equal to $\lfloor n/2 \rfloor$ and a location $l$ is an integer between $r+1$ and $n-r$. This couple $(l, r)$ refers to the segment $[l-r, l+r)$ centered at $l$ and with radius $r$. Formally, $\mathcal{G}$ is therefore a subset of $J_n = \{(l,r): r = 1, \dots, \lfloor \frac{n}{2} \rfloor \text{ and } l = r+1, \dots, n-r+1\}$. Given a grid $\mathcal{G}$, we call $\mathcal{R}$ its collection of scales, that is $\mathcal{R} = \{r: \exists l \text{ s.t. } (l,r) \in \mathcal{G}\}$. Finally, for a scale $r \in \mathcal{R}$, $D_r$ stands for the corresponding collection of locations, that is $D_r = \{l: (l,r) \in \mathcal{G}\}$. Although we do not make any assumption on the grid $\mathcal{G}$ for the time being, we will mainly consider two specific grids in this section: the complete grid $G_F = J_n$ and the dyadic grid $G_D$ defined by $\mathcal{R} = \{1, 2, 4, \dots, 2^{\lfloor \log_2(n) \rfloor - 1}\}$, $D_1 = [2, n]$, and

$$ D_r = \left\{ r + 1, 3\lfloor r/2 \rfloor + 1, 4\lfloor r/2 \rfloor + 1, \dots, \left( \frac{n}{\lfloor r/2 \rfloor} - 2 \right) \lfloor \frac{r}{2} \rfloor + 1, n - r + 1 \right\} \quad \text{for } r \in \mathcal{R} \setminus \{1\}. \quad (5) $$

See Figure 2 for a visual representation of the dyadic grid. At some points, we shall also mention a-adic grids $G_a$. For any $a \in (0, 1)$, $G_a$ is defined by $\mathcal{R} = \{1, \lfloor a^{-1} \rfloor, \lfloor a^{-2} \rfloor, \dots, \lfloor a^{1-\lfloor \log(n)/ \log(a) \rfloor} \rfloor\}$ and $D_r$ as in (5). Interestingly, the cardinality of the dyadic grid or more generally of the a-adic grid is order $O(n)$, whereas the complete grid $G_D$ is quadratic.

Grids are reminiscent of the c-normal systems of intervals introduced by Nemirovsky [30] (see also [27] for a definition) although our definition allows for non-necessarily normal intervals.

Given a fixed grid $\mathcal{G}$, a multiscale test is simply a collection of test $T = (T_{l,r})_{(l,r) \in \mathcal{G}}$ indexed by the elements of $\mathcal{G}$, which amounts to testing at all scales $r \in \mathcal{R}$ and all locations $l \in D_r$ whether the functional $\Gamma(\mathbb{P}_t)$ is constant over the segment $[l-r, l+r)$. Equivalently, $T_{l,r}$ tests whether there exists a change-point in $[l-r+1, l+r-1]$.

### 2.2 From a multiscale test to a change-point detection procedure

Our purpose is to introduce a generic procedure to translate a multiscale procedure into a vector of change-points. Intuitively, if, for some $(l, r) \in \mathcal{G}$, we have $T_{l,r} = 1$, then the functional $\Gamma(\mathbb{P}_t)$ is certainly not constant over $[l-r, l+r)$ which entails that there is possibly at least one change-point in $[l-r+1, l+r-1]$. As a consequence, the multiscale test gives a collection $\mathcal{I}(T) = \{[l-r+1, l+r-1] \text{ s.t. } T_{l,r} = 1\}$ of intervals that tentatively contain at least one change-point.

If all these intervals were disjoint, then one simply would take $\hat{\tau}$ as the sequence of centers of these intervals. Unfortunately, when two intervals $[l_1 - r_1 + 1, l_1 + r_1 - 1]$ and $[l_2 - r_2 + 1, l_2 + r_2 - 1]$ in $\mathcal{I}(T)$ have a non-empty intersection, one cannot necessarily decipher whether there is only one change-point in the intersection of both intervals or if each interval contains a specific change-point. Hence, our general objective is to transform the collection $\mathcal{I}(T)$ into a collection of non-intersecting intervals by either discarding or merging some of them.
This statistic computes the normalized difference of empirical mean of $y_i$ on $[l-r, l)$ and $[l, l+r)$. If the noise is Gaussian and if $\Theta$ is constant on $[l-r, l+r)$, then $C_{l,r}(Y)$ simply follows a standard $p$-dimensional normal distribution. To simplify, consider a specific instance of our testing problem where we want to test whether $\{\Theta$ is constant over $[l-r, l+r)\}$ versus $\{\Theta$ contains exactly one change-point at $l$ on the segment $[l-r, l+r)\}$. This corresponds to a two-sample mean testing problem, for which the CUSUM statistic $C_{l,r}(Y)$ is a sufficient statistic if the noise is Gaussian. Then, given $C_{l,r}(Y)$, one wants to test whether its expectation is 0 (no change-point on $[l-r, l+r)$) versus its expectation is non-zero but is $s$-sparse for some unknown $s$. This classical detection problem is well understood [11] and it is well known that a combination of a $\chi^2$-type test with a higher-criticism-type test is optimal. Here, the challenge stems from the fact that we do not want to perform a single such test, but a large collection of tests over a collection of $(l,r) \in \mathcal{G}$.

### 1.2.2. Our contribution

As usual in the mean change-point literature, we consider the energy $r_k\Delta_k^2$ of the change-point $\tau_k$. Up to a possible factor in $[1/2, 1]$, $r_k\Delta_k^2$ is the square distance between $\Theta$ and its projection on the space of vectors $\Theta'$ with change-point at $(\tau_1, \dots, \tau_{k-1}, \tau_{k+1}, \dots, \tau_K)$ – see e.g. [38] for a discussion in the univariate setting. In other words, the energy $r_k\Delta_k^2$ characterizes the significance of the change-point $\tau_k$. In Section 3, we introduce a multi-scale change-point detection procedure detecting any change-point $\tau_k$ whose energy is higher, up to a numerical constant, than $\sigma^2 s_k \log(1 + \frac{\sqrt{p}}{s_k} \sqrt{\log(n/r_k)}) + \sigma^2 \log(n/r_k)$. This result is valid for arbitrary length $r_k$ and sparsity $s_k$, and does not require the knowledge of these two quantities. In summary, our procedure does not estimate any spurious change-point (**NoSp**) and **detects** all the change-points whose energy are higher than the latter threshold. In Section 5, we establish that, as soon as the unknown number $K$ of the change-points is larger than 1, the condition $\sigma^2 s_k \log(1 + \frac{\sqrt{p}}{s_k} \sqrt{\log(n/r_k)}) + \sigma^2 \log(n/r_k)$ on the energy is tight with respect to $n$, $p$, $r_k$ and $s_k$, in the sense that no procedure achieving (**NoSp**) is able to detect with high probability a change-point whose energy is smaller (up to some constant) than the latter threshold. In Section 4, we consider the more general setting where the noise is $L$-sub-gaussian with known variance, and we establish a similar result to the Gaussian case up to a logarithmic loss in some regimes. Finally, we illustrate in Section 8 the behavior of our procedure on numerical experiments.

### 1.2.3. Related work

For dense change-points ($s_k = p$) but with unknown covariance for the noise, Wang et al. [45] (see also [44]) study the behavior of a procedure based on $U$-statistics of the CUSUM. Jirak [21] and Yu and Chen [48] introduce binary segmentation procedures based on the $l_\infty$ norm of the CUSUMs. Although those
<figure><img src="image_3.png" /><figcaption>Figure 2: The dyadic grid is represented as follows : for each $r = 2^i$ and $l \in D_r$, we draw the interval $[l - r + 1, l + r - 1]$ at position $(l, \log_2(r))$.</figcaption></figure>

We propose the following bottom-up iterative procedure for building a collection of non-intersecting intervals. Start with $\mathcal{T}_0 = \mathcal{S}_0 = \emptyset$. For any scale $r \in \mathcal{R}$, we compute the collections $\mathcal{S}_r$ of intervals of scale $r$ and the collection $\mathcal{T}_r$ of locations based on the following

$$ \mathcal{T}_r = \left\{ l \in D_r, \ T_{l,r} = 1 \quad \text{and} \quad [l-r+1, l+r-1] \cap \left( \bigcup_{r' < r, r' \in \mathcal{R}} \mathcal{S}_{r'} \right) = \emptyset \ ; \right\} \\ \mathcal{S}_r = \bigcup_{l \in \mathcal{T}_r} [l-r+1, l+r-1] . $$

The sets $\mathcal{T}_1$ and $\mathcal{S}_1$ are made of all positions $l$ such that $T_{l,1} = 1$. More generally, $\mathcal{T}_r$ contains all locations $l$ such that $T_{l,r} = 1$ and the corresponding interval $[l - r + 1, l + r - 1]$ does not intersect with any of the detected intervals at a smaller scale $r' < r$. The set $\mathcal{S}_r$ contains all intervals associated to $\mathcal{T}_r$.

One can easily check that $\mathcal{S} = \bigcup_r \mathcal{S}_r$ is a union of closed non-intersecting intervals. Denote $\mathcal{C} = \{C_1, \dots, C_{\hat{K}}\}$ the partition of $\mathcal{S}$ into connected components such that, for all $1 \le i < j \le \hat{K}$, $\max C_i < \min C_j$. Finally, we estimate the vector of change-points $\hat{\tau}$ by taking the center of each segment $C_k$. In other words, we take $\hat{\tau}_k := \frac{1}{2}(\min C_k + \max C_k)$ for any $1 \le k \le \hat{K}$. This bottom-up aggregation procedure is summarized in Algorithm 1 and illustrated in Figure 3 below.

**Remark:** If, for some $r \in \mathcal{R}$ and some $l_1 < l_2 \in D_r$, we have $T_{l_1,r} = 1$, $T_{l_2,r} = 1$, and $l_1 + r - 1 \ge l_2 - r + 1$, then $\mathcal{S}_r$ contains the segment $[l_1 - r + 1, l_2 + r - 1]$. In other words, our aggregation procedure merges two intervals if and only if they correspond to the same scales. In Section A, we also introduce a variant of the algorithm where, instead of merging these two intersecting with identical scale, we discard one of them.

**Computational Cost.** A naive implementation of Algorithm 1 - and also of Algorithm 2 defined in Appendix - requires to compute all tests $T_{l,r}$ on the grid, whereas the aggregation procedure only needs to compute a number of tests $T_{l,r}$ proportional to the size of the grid. More precisely, if the
work explicitly characterize the asymptotic distribution of the test statistics and, for some of them, allow temporal dependencies in the data, the corresponding energy requirements for change-point detection are either not studied or turn out to be suboptimal.

Closest to our work, Chan and Chen [5] study a bottom-up approach to detect change-points of a Gaussian multivariate time series in an asymptotic setting. More specifically, the authors consider an asymptotic regime where the size of the time series is exponential in the dimension: $n = e^{p\zeta}$ with $\zeta \in (0, 1)$. The authors also assume that the number $K$ of change-points remains finite when $n, p \to \infty$ and that the minimal sparsity $s$ of these change-points is polynomial is $p$. In this specific regime, their procedures provably recover change-points under a near-minimal (up to logarithmic factors with respect to $n$) condition on the energy. In contrast, our results provide non-asymptotic and tight results for all scaling with respect to $n$ and $p$, allow for arbitrarily large number $K$ of change-points and allow for the presence of non-significant change-points. In the same specific asymptotic setting, [20] introduce a so called score test statistic used in a change-point detection procedure which is shown to achieve the same performance as [5] in the gaussian model but also handle Poisson observations.

Recently, Liu et al. [28] have characterized the optimal detection rate of a possibly sparse change-point in the specific case where there is at most one change-point, but the optimal rates are significantly slower in the multiple change-point setting. See also [12] and [9] for earlier results. Wang and Samworth [46] have proposed the INSPECT method based on sparse projection to handle sparse change-points, but INSPECT provably detects the change-points under strong assumption on the energy; see Section 3 for a precise comparison.

In the univariate setting ($p = 1$), minimal energy requirements for change-point detection are well understood [13, 15, 38, 42] and are nearly achieved by a wide range of procedures including penalized least-square and multi-scale tests methods.

## 2. A Generic algorithm for multiscale change-point detection on a grid

In this section, we study the problem of change-point detection in the general setting defined in Section 1.1. We introduce a bottom-up algorithm that aggregates a collection of homogeneity tests, performed at many positions, and for many scales, of our data. Then, we establish that, under some conditions on these tests, the procedure detects significant change-points.

### 2.1. Grid and multiscale statistics

Since our purpose is to translate a collection of local tests $T = (T_{l,r})_{(l,r) \in \mathcal{G}}$ indexed by a grid $\mathcal{G}$ into a change-point detection procedure, we first need to formalize what we mean by a grid. Henceforth, we call a grid $\mathcal{G}$ of $[n]$ a collection of locations and scales where a scale $r$ is a positive integer smaller or equal to $\lfloor n/2 \rfloor$.
**Data:** $y_t, t = 1...n$ and local test statistics $(T_{l,r})_{(l,r) \in \mathcal{G}}$

**Result:** $(\hat{\tau}_k)_{k \le \hat{K}}$

$\mathcal{T}_r, \mathcal{S}_r = \emptyset$ for all $r \in \mathcal{R}$ and $\mathcal{S} = \emptyset$;

for increasing $r \in \mathcal{R}$ do

for $l \in \mathcal{D}_r$ s.t. $T_{l,r} = 1$ do
    if $[l-r+1, l+r-1] \cap \mathcal{S} = \emptyset$ then
        $\mathcal{T}_r \leftarrow \mathcal{T}_r \cup \{l\}$;
        $\mathcal{S}_r \leftarrow \mathcal{S}_r \cup [l-r+1, l+r-1]$;
    end
end
$\mathcal{S} = \mathcal{S} \cup \mathcal{S}_r$;

end

Let $(C_k)_{k=1,...,\hat{K}}$ be the connected components of $\mathcal{S}$ sorted in increasing order;

return $(\hat{\tau}_k = \frac{1}{2}(\min C_k + \max C_k))_{k=1,...,\hat{K}}$

**Algorithm 1:** Bottom-up aggregation procedure of multiscale tests

computational cost of $T_{l,r}$ is $\Lambda_{l,r}$ for each $(l, r)$ in the grid $\mathcal{G}$, then the aggregation procedure requires $O(\sum_{(l,r) \in \mathcal{G}} \Lambda_{l,r})$ computations. If for all $(l, r)$, the cost $\Lambda_{l,r}$ is proportional to $r$, that is $\Lambda_{l,r} = O(r\Lambda)$, then the overall computational cost is $O(\Lambda \sum_{(l,r) \in \mathcal{G}} r)$ which is $O(\Lambda n^3)$ for the complete grid and $O(\Lambda n \log(n))$ for the dyadic grid. One can speed up the full procedure by computing the statistics $T_{l,r}$ and aggregating on the fly by checking whether $[l-r+1, l+r-1]$ intersects $\mathcal{S}$ before evaluating $T_{l,r} = 1$. Indeed, the connected components $C_k$ can be computed at each increasing scale $r$. Hence, at scale $r$, one only needs to compute the tests $T_{l,r}$ at locations $l$ such that $[l-r+1, l+r-1]$ does not intersect the connected components detected at scales $r' < r$.

## 2.3 General analysis

In this subsection, we provide an abstract theorem translating error controls of the multiple test procedure $T$ in terms of properties of $\hat{\tau}$. As explained in the introduction, the time series $(y_t)$ may contain change-points that are too small to be detected. Having this in mind, we define a subset $\mathcal{K}^* \subset [K]$ of indices corresponding to so-called significant change-points. As our purpose is to provide deterministic condition so that the change-points in $\mathcal{K}^*$, we need to introduce, for each $k \in \mathcal{K}^*$, an element of the grid $(\bar{\tau}_k, \bar{r}_k) \in \mathcal{G}$ at which the statistic $T$ is expected to detect $\tau_k$. One could think of $\bar{\tau}_k$ as some position close to $\tau_k$ and to $\bar{r}_k$ as some radius which is large enough to convey information on the change-point. Recall that the length $r_k$ of the change-point $\tau_k$ is defined by $r_k = \min(\tau_{k+1} - \tau_k, \tau_k - \tau_{k-1})$. We assume that the scales $\bar{r}_k$ and the location $\bar{\tau}_k$ of detection satisfy the two following conditions:

$$4(\bar{r}_k - 1) < r_k \quad \text{and} \quad |\bar{\tau}_k - \tau_k| \le \bar{r}_k - 1. \qquad (6)$$

The first condition ensures that the scale $\bar{r}_k < r_k/4 + 1$ is small enough compared to the length $r_k$. The second condition is always satisfied if $\bar{\tau}_k$ is the best approximation of $\tau_k$ in $\mathcal{D}_{\bar{r}_k}$ and if the grid $\mathcal{G}$ satisfies the following approximation property

(App): For all $r \in \mathcal{R}$ and all $l \in [r+1, n-r+1]$, there exists $l' \in \mathcal{D}_r$ such that $|l' - l| \le r - 1$.

This property entails that any point $l$ can be approximated at distance $r-1$ by some location in $\mathcal{D}_r$. This also implies that each point $l \in [r+1, n-r]$ belongs to at least one segment $(l' - r, l' + r)$ where $l_1$ lies in $\mathcal{D}_r$. In practice, the a-adic grids $\mathcal{G}_a$ and the complete grid satisfy (App).
<figure><img src="image_3.png" /><figcaption>FIG 2. The dyadic grid is represented as follows: for each $r = 2^i$ and $l \in D_r$, we draw the interval $[l - r + 1, l + r - 1]$ at position $(l, \log_2(r))$.</figcaption></figure>

and a location $l$ is an integer between $r+1$ and $n-r$. This couple $(l, r)$ refers to the segment $[l-r, l+r)$ centered at $l$ and with radius $r$. Formally, $\mathcal{G}$ is therefore a subset of $J_n = \{(l, r) : r = 1, \dots, \lfloor \frac{n}{2} \rfloor \text{ and } l = r+1, \dots, n-r+1\}$. Given a grid $\mathcal{G}$, we call $\mathcal{R}$ its collection of scales, that is $\mathcal{R} = \{r : \exists l \text{ s.t. } (l, r) \in \mathcal{G}\}$. Finally, for a scale $r \in \mathcal{R}$, $D_r$ stands for the corresponding collection of locations, that is $D_r = \{l : (l, r) \in \mathcal{G}\}$. Although we do not make any assumption on the grid $\mathcal{G}$ for the time being, we will mainly consider two specific grids in this section: the **complete** grid $\mathcal{G}_F = J_n$ and the **dyadic** grid $\mathcal{G}_D$ defined by $\mathcal{R} = \{1, 2, 4, \dots, 2^{\lfloor \log_2(n) \rfloor - 1}\}$, $D_1 = [2, n]$, and for $r \in \mathcal{R} \setminus \{1\}$,

$$D_r = \left\{ r+1, 3\lfloor r/2 \rfloor + 1, 4\lfloor r/2 \rfloor + 1, \dots, \left( \lfloor \frac{n}{r/2} \rfloor - 2 \right) \lfloor \frac{r}{2} \rfloor + 1, n-r+1 \right\}. \quad (5)$$

See Figure 2 for a visual representation of the dyadic grid. At some points, we shall also mention $a$-adic grids $\mathcal{G}_a$. For any $a \in (0, 1)$, $\mathcal{G}_a$ is defined by $\mathcal{R} = \{1, \lfloor a^{-1} \rfloor, \lfloor a^{-2} \rfloor, \dots, \lfloor a^{1-\lfloor \log(n)/\log(a) \rfloor} \rfloor\}$ and $D_r$ as in (5). Interestingly, the cardinality of the dyadic grid or more generally of the $a$-adic grid is order $O(n)$, whereas the complete grid $\mathcal{G}_D$ is quadratic.

Grids are reminiscent of the $c$-normal systems of intervals introduced by Nemirovsky [30] (see also [27] for a definition) although our definition allows for non-necessarily normal intervals.

Given a fixed grid $\mathcal{G}$, a multiscale test is simply a collection of test $T = (T_{l,r})_{(l,r) \in \mathcal{G}}$ indexed by the elements of $\mathcal{G}$, which amounts to testing at all scales $r \in \mathcal{R}$ and all locations $l \in D_r$ whether the functional $\Gamma(\mathbb{P}_l)$ is constant over the segment $[l-r, l+r)$. Equivalently, $T_{l,r}$ tests whether there exists a change-point in $[l-r+1, l+r-1]$.
<figure><img src="image_4.png" /><figcaption>Figure 3: Example of our change-point detection procedure with three change-points. The first two change-points have large heights and are detected at a small scale r (in magenta) while the third one is detected at a larger scale r.</figcaption></figure>

Next, we introduce an event on the tests $(T_{l,r})$ under which the change-point estimator $\hat{\tau}$ of Algorithm 1 performs well. In the following, we write $\mathcal{H}_0$, the collection of all possible $(l,r) \in J_n$ such that there is no change in $[l-r+1, l+r-1]$, i.e. $\Gamma(\mathbb{P}_t)$ is constant on $[l-r, l+r)$. Equivalently, we have

$$ (l, r) \in \mathcal{H}_0 \quad \text{iff} \quad (l-r, l+r) \cap \{\tau_k, k = 1, \dots, K\} = \emptyset. \qquad (7) $$

For a collection $\mathcal{K}^*$ and some elements of the grid $(\bar{\tau}_k, \bar{r}_k)$ satisfying (6), the Event $\mathcal{A}(T, \mathcal{K}^*, (\bar{\tau}_k, \bar{r}_k)_{k \in \mathcal{K}^*})$ is defined as the conjunction of the two following properties: (i) (No false positive) $T_{l,r} = 0$ for all $(l,r) \in \mathcal{H}_0 \cap \mathcal{G}$ (ii) (Detection of significant change-points) for every $k \in \mathcal{K}^*$, we have $T_{\bar{\tau}_k, \bar{r}_k} = 1$.

The first property states that T performs no type I errors on the event $\mathcal{A}(T, \mathcal{K}^*, (\bar{\tau}_k, \bar{r}_k)_{k \in \mathcal{K}^*})$, whereas the second property enforces that all the significant change-points are detected by the specific tests $T_{\bar{\tau}_k, \bar{r}_k}$.

**Theorem 1.** The following holds for any grid $\mathcal{G}$, any local test statistic T, any non-negative integer K, any distribution with K change-points, any $\mathcal{K}^* \subset [K]$ and scales and locations $(\bar{\tau}_k, \bar{r}_k)_{k \in \mathcal{K}^*}$ in $\mathcal{G}$ satisfying Assumption (6). Under the event $\mathcal{A}(T, \mathcal{K}^*, (\bar{\tau}_k, \bar{r}_k)_{k \in \mathcal{K}^*})$, the estimated change-point vector $\hat{\tau}$ returned by Algorithm 1 satisfies the two following properties

* **Significant change-points are detected:** for all $k \in \mathcal{K}^*$, there exists $k' \le \hat{K}$ such that $|\hat{\tau}_{k'} - \tau_k| \le \bar{r}_k - 1 < \frac{r_k}{4}$.

* (NoSp): No Spurious change-point is detected (1).

The first property states that so-called significant change-points $(\tau_k)_{k \in \mathcal{K}^*}$ are detected by the generic algorithm at the right scale. The no-spurious property (1) guarantees that, around any true change-point $\tau_k$, the procedure estimates at most one single change-point $\hat{\tau}_k$. Importantly, the theorem does not make any assumption on the non-significant change-points. In fact, change-points
## 2.2. From a multiscale test to a change-point detection procedure

Our purpose is to introduce a generic procedure to translate a multiscale procedure into a vector of change-points. Intuitively, if, for some $(l, r) \in \mathcal{G}$, we have $T_{l,r} = 1$, then the functional $\Gamma(\mathbb{P}_t)$ is certainly not constant over $[l-r, l+r)$ which entails that there is possibly at least one change-point in $[l-r+1, l+r-1]$. As a consequence, the multiscale test gives a collection $\mathcal{I}(T) = \{[l-r+1, l+r-1]\}$ s.t. $T_{l,r} = 1$} of intervals that tentatively contain at least one change-point.

If all these intervals were disjoint, then one simply would take $\hat{\tau}$ as the sequence of centers of these intervals. Unfortunately, when two intervals $[l_1 - r_1 + 1, l_1 + r_1 - 1]$ and $[l_2 - r_2 + 1, l_2 + r_2 - 1]$ in $\mathcal{I}(T)$ have a non-empty intersection, one cannot necessarily decipher whether there is only one change-point in the intersection of both intervals or if each interval contains a specific change-point. Hence, our general objective is to transform the collection $\mathcal{I}(T)$ into a collection of non-intersecting intervals by either discarding or merging some of them.

We propose the following bottom-up iterative procedure for building a collection of non-intersecting intervals. Start with $\mathcal{T}_0 = \mathcal{S}_0 = \emptyset$. For any scale $r \in \mathcal{R}$, we compute the collections $\mathcal{S}_r$ of intervals of scale $r$ and the collection $\mathcal{T}_r$ of locations based on the following

$$ \begin{aligned} \mathcal{T}_r &= \left\{ l \in \mathcal{D}_r, \ T_{l,r} = 1 \quad \text{and} \quad [l-r+1, l+r-1] \cap \left( \bigcup_{r' < r, r' \in \mathcal{R}} \mathcal{S}_{r'} \right) = \emptyset \right\} \\ \mathcal{S}_r &= \bigcup_{l \in \mathcal{T}_r} [l-r+1, l+r-1]. \end{aligned} $$

The sets $\mathcal{T}_1$ and $\mathcal{S}_1$ are made of all positions $l$ such that $T_{l,1} = 1$. More generally, $\mathcal{T}_r$ contains all locations $l$ such that $T_{l,r} = 1$ and the corresponding interval $[l-r+1, l+r-1]$ does not intersect with any of the detected intervals at a smaller scale $r' < r$. The set $\mathcal{S}_r$ contains all intervals associated to $\mathcal{T}_r$.

One can easily check that $\mathcal{S} = \bigcup_r \mathcal{S}_r$ is a union of closed non-intersecting intervals. Denote $\mathcal{C} = \{C_1, \dots, C_{\hat{K}}\}$ the partition of $\mathcal{S}$ into connected components such that, for all $1 \le i < j \le \hat{K}$, $\max C_i < \min C_j$. Finally, we estimate the vector of change-points $\hat{\tau}$ by taking the center of each segment $C_k$. In other words, we take $\hat{\tau}_k := \frac{1}{2}(\min C_k + \max C_k)$ for any $1 \le k \le \hat{K}$. This bottom-up aggregation procedure is summarized in Algorithm 1 and illustrated in Figure 3 below.

**Remark.** If, for some $r \in \mathcal{R}$ and some $l_1 < l_2 \in \mathcal{D}_r$, we have $T_{l_1,r} = 1$, $T_{l_2,r} = 1$, and $l_1 + r - 1 \ge l_2 - r + 1$, then $\mathcal{S}_r$ contains the segment $[l_1 - r + 1, l_2 + r - 1]$. In other words, our aggregation procedure merges two intervals if and only if they correspond to the same scales. In Section A, we also introduce a variant of the algorithm where, instead of merging these two intersecting with identical scale, we discard one of them.

**Computational Cost.** A naive implementation of Algorithm 1 – and also of Algorithm 2 defined in Appendix – requires to compute all tests $T_{l,r}$ on the grid,
$\tau_k$ with $k \in [K] \setminus K^*$ may or may not be detected. In general, we can only conclude from Theorem 1 that $|K^*| \le \hat{K} \le K$ on the event $\mathcal{A}(T, K^*, (\bar{\tau}_k, \bar{r}_k)_{k \in K^*})$.

Theorem 1 is abstract but its main virtue is to translate multiple testing properties into change-point detection properties. For a specific problem such as multivariate mean change-point detection considered in the next section, the construction of a near optimal procedure boils down to introducing a collection of local test statistics, such that (a) change-points $\tau_k$ belong to $K^*$ under minimal conditions, (b) the scale $\bar{r}_k$ is the smallest possible, and (c) the event $\mathcal{A}(T, K^*, (\bar{\tau}_k, \bar{r}_k)_{k \in K^*})$ holds with high probability.

In the case where all the change-points are significant, the result of Theorem 1 can be reformulated as follows:

**Corollary 1.** The following holds for any grid $\mathcal{G}$, any local test statistic $T$, any non-negative integer $K$, any distribution with $K$ change-points, any $(\bar{\tau}_k, \bar{r}_k)_{k=1,...,K}$ in $\mathcal{G}$ satisfying Assumption (6). Under the event $\mathcal{A}(T, [K], (\bar{\tau}_k, \bar{r}_k)_{k=1,...,K}$, the estimated change-point vector $\hat{\tau}$ returned by Algorithm 1 satisfies $\hat{K} = K$ and,

$$|\hat{\tau}_k - \tau_k| < \bar{r}_k - 1 \le \frac{r_k}{4} \quad \text{for all } k = 1, \dots, K.$$

Let us respectively define the Hausdorff distance and the Wasserstein distance of two vectors $(u_1, \dots, u_K)$ and $(v_1, \dots, v_K)$ in $\mathbb{R}^K$ by $d_H(u, v) = \max_{k=1,...,K} |u_k - v_k|$ and $d_W(u, v) = \sum_{k=1,...,K} |u_k - v_k|$. Then, Corollary 1 straightforwardly implies that, if $K^* = [K]$, then these two losses are bounded as follows

$$d_H(\hat{\tau}, \tau) \le \max_{k=1,...,K} (\bar{r}_k - 1) \quad \text{and} \quad d_W(\hat{\tau}, \tau) \le \sum_{k=1,...,K} (\bar{r}_k - 1).$$

As an alternative of Algorithm 1, one could use other bottom-up aggregating procedures. For instance, Algorithm 2 defined in Appendix A also satisfies Theorem 1. Although these two algorithms are closely related, Algorithm 1 is slightly more conservative than Algorithm 2 since it merges all detection intervals at a given resolution while Algorithm 2 only keeps one interval at a given resolution when multiple intervals intersect - the one with smallest index $t$. While the min-max properties of both methods are comparable - at least up to a multiple constant - the choice of aggregation method will have an influence in practice on the outcome: Algorithm 1 will be slightly more stable, detect less change-points, and provide wider confidence interval around them, while Algorithm 2 will be slightly more sensitive to smaller changes, i.e. detect smaller change-points, will be more precise, and somewhat less stable.

Theorem 1 ensures that, if $T_{\bar{\tau}_k, \bar{r}_k} = 1$ with $(\bar{\tau}_k, \bar{r}_k)$ satisfying Assumption (6), then the change-point $\tau_k$ is detected. Inspecting the proof of Theorem 1, one easily checks that Assumption (6) is minimal for Algorithm 1 (and also for Algorithm 2). Still, one may wonder whether any generic algorithm has to require that $4(\bar{r}_k - 1) < r_k$ to detect the change-points or if there exists a generic algorithm where the constant 4 in the above condition can be improved.

**Comparison with narrowest over threshold methods.** As mentioned in the introduction, other aggregation procedures have been proposed in the literature. In particular, the narrowest over threshold scheme proposed by [2] and later used in [24] is also closely related to the local segmentation algorithm of Chan and Chen [5]. A simple extension of these procedures for generic change-point problems and for a general collection of tests ($T_{l,r}$) would amount to modifying Algorithm 1 by selecting locations $l$ in $\mathcal{D}_r$ such that $T_{l,r} = 1$ and $[l - r + 1, l + r - 1]$ does not intersect
**Data:** $y_t, t = 1...n$ and local test statistics $(T_{l,r})_{(l,r) \in \mathcal{G}}$

**Result:** $(\hat{\tau}_k)_{k \le \hat{K}}$

$T_r, S_r = \emptyset$ for all $r \in \mathcal{R}$ and $S = \emptyset$;

for increasing $r \in \mathcal{R}$ do

for $l \in \mathcal{D}_r$ s.t. $T_{l,r} = 1$ do
    if $[l-r+1, l+r-1] \cap S = \emptyset$ then
      $T_r \leftarrow T_r \cup \{l\}$;
      $S_r \leftarrow S_r \cup [l-r+1, l+r-1]$;
    end
  end
  $S = S \cup S_r$;

end

Let $(C_k)_{k=1,...,\hat{K}}$ be the connected components of $S$ sorted in increasing order;

return $(\hat{\tau}_k = \frac{1}{2}(\min C_k + \max C_k))_{k=1,...,\hat{K}}$

<figure><img src="image_4.png" /><figcaption>FIG 3. Example of our change-point detection procedure with three change-points. The first two change-points have large heights and are detected at a small scale $r$ (in magenta) while the third one is detected at a larger scale $r$.</figcaption></figure>

whereas the aggregation procedure only needs to compute a number of tests $T_{l,r}$ proportional to the size of the grid. More precisely, if the computational cost of $T_{l,r}$ is $\Lambda_{l,r}$ for each $(l,r)$ in the grid $\mathcal{G}$, then the aggregation procedure requires $O(\sum_{(l,r) \in \mathcal{G}} \Lambda_{l,r})$ computations. If for all $(l,r)$, the cost $\Lambda_{l,r}$ is proportional to $r$, that is $\Lambda_{l,r} = O(r\Lambda)$, then the overall computational cost is $O(\Lambda \sum_{(l,r) \in \mathcal{G}} r)$ which is $O(\Lambda n^3)$ for the complete grid and $O(\Lambda n \log(n))$ for the dyadic grid. One can speed up the full procedure by computing the statistics $T_{l,r}$ and aggregating on the fly by checking whether $[l-r+1, l+r-1]$ intersects $S$ before evaluating $T_{l,r} = 1$. Indeed, the connected components $C_k$ can be computed at each increasing scale $r$. Hence, at scale $r$, one only needs to compute the tests
previously detected change-points, whereas we require in Algorithms 1 and 2, that $[l-r+1, l+r-1]$ does not intersect previously detected confidence intervals. In some way, the narrowest-over threshold scheme is therefore less conservative. Unfortunately, there is no generic result in the form of Theorem 1 for such procedures and, from informal arguments, we doubt that the corresponding procedure provably achieves (NoSp) under a control of the FWER of the tests. Inspecting the proof of Theorem 1 in [2] and Theorem 3 in [24] for univariate mean change-point problems, one observes that the chosen threshold is much larger than what is needed to control the FWER so that the theoretical threshold is certainly over-conservative – see step 5 of the proof of Theorem 1 in [2]. In contrast, Theorem 1 in [5] for univariate change-point problems is based on the minimal threshold, but the proof relies on the important assumption that the number $K$ of change-point remains bounded while $n$ goes to infinity. Besides, it is not clear how one could extend the arguments to more general settings.

## 3 Multivariate Gaussian change-point detection

We now turn to the multivariate change-point model introduced in Section 1.2. Throughout this section, we assume that the random vectors $\epsilon_t$ are independently and identically distributed with $\epsilon_t \sim N(0, \sigma^2 I_p)$. Since we shall apply the general aggregation procedures introduced in the previous section, our main job here is to introduce a near-optimal testing procedure.

Fix some quantity $\delta \in (0, 1)$. At the end of the section, $1 - \delta$ will correspond to the probability of the event $A(T, K^*, (\bar{\tau}_k, \bar{r}_k)_{k \in K^*})$ introduced in the previous section. Alternatively, one may interpret $\delta$ as an upper bound of the desired probability that the change-point detection procedure detects a spurious change-points. Recall that, for a change-point $\tau_k$, $s_k$ stands for the sparsity of the difference $\mu_{k+1} - \mu_k$. The energy of a given change-point $\tau_k$ is $c_0$-high if

$$r_k \Delta_k^2 \ge c_0 \sigma^2 \left[ s_k \log \left( 1 + \frac{\sqrt{p}}{s_k} \sqrt{\log \left( \frac{n}{r_k \delta} \right)} \right) + \log \left( \frac{n}{r_k \delta} \right) \right], \quad (8)$$

for some universal constant $c_0$ to be defined later. We show in this section that when $c_0$ is large enough, all high-energy change-points can be detected. Conversely, it is established in Section 5 that Condition (8) is (up to a multiplicative constant) optimal for detecting change-points and cannot be weakened.

Let us now discuss the different regimes contained in Equation (8). In what follows, define

$$\psi_{n,r,s}^{(g)} := s \log \left( 1 + \frac{\sqrt{p}}{s} \sqrt{\gamma_r} \right) + \gamma_r ; \quad \gamma_r := \log \left( \frac{n}{r\delta} \right),$$

in order to alleviate notations. If $\gamma_r \ge p/2$, then $\psi_{n,r,s}^{(g)} \asymp \gamma_r$ where $u \asymp v$ means that for two positive numerical constants $c_1$ and $c_2$, one has $c_1 v \le u \le c_2 v$. This corresponds to the minimal energy condition for detection in the univariate case, i.e. when $p = 1$; see [38]. The condition $\gamma_r \ge p/2$ occurs when $p$ is rather small and the scale $r$ is much smaller than $n$. If $\gamma_r \le p/2$, then

$$\psi_{n,r,s}^{(g)} \asymp \begin{cases} \gamma_r & \text{if } s \le \frac{\gamma_r}{\log(p)-\log(\gamma_r)} \\ s \log \left( 2 \frac{p}{s^2} \gamma_r \right) / \sqrt{p \gamma_r} & \text{if } \frac{\gamma_r}{\log(p)-\log(\gamma_r)} < s < \sqrt{p \gamma_r} \\ \gamma_r & \text{if } s \ge \sqrt{p \gamma_r} . \end{cases}$$

We define $K^* \subset [K]$ as the subset of indices such that $\tau_k$ satisfies (8). For any $k \in K^*$, we define
$T_{l,r}$ at locations $l$ such that $[l-r+1, l+r-1]$ does not intersect the connected components detected at scales $r' < r$.

## 2.3. General analysis

In this subsection, we provide an abstract theorem translating error controls of the multiple test procedure $T$ in terms of properties of $\hat{\tau}$. As explained in the introduction, the time series $(y_t)$ may contain change-points that are too small to be detected. Having this in mind, we define a subset $K^* \subset [K]$ of indices corresponding to so-called significant change-points. As our purpose is to provide deterministic condition so that the change-points in $K^*$, we need to introduce, for each $k \in K^*$, an element of the grid $(\bar{\tau}_k, \bar{r}_k) \in G$ at which the statistic $T$ is expected to detect $\tau_k$. One could think of $\bar{\tau}_k$ as some position close to $\tau_k$ and to $\bar{r}_k$ as some radius which is large enough to convey information on the change-point. Recall that the length $r_k$ of the change-point $\tau_k$ is defined by $r_k = \min(\tau_{k+1} - \tau_k, \tau_k - \tau_{k-1})$. We assume that the scales $\bar{r}_k$ and the location $\bar{\tau}_k$ of detection satisfy the two following conditions:

$$4(\bar{r}_k - 1) < r_k \quad \text{and} \quad |\bar{\tau}_k - \tau_k| \le \bar{r}_k - 1. \qquad (6)$$

The first condition ensures that the scale $\bar{r}_k < r_k/4+1$ is small enough compared to the length $r_k$. The second condition is always satisfied if $\bar{\tau}_k$ is the best approximation of $\tau_k$ in $D_{\bar{r}_k}$ and if the grid $G$ satisfies the following approximation property

(App): For all $r \in \mathbb{R}$ and all $l \in [r+1, n-r+1]$, there exists $l' \in D_r$ such that $|l' - l| \le r - 1$.

This property entails that any point $l$ can be approximated at distance $r-1$ by some location in $D_r$. This also implies that each point $l \in [r+1, n-r]$ belongs to at least one segment $(l' - r, l' + r)$ where $l_1$ lies in $D_r$. In practice, the a-adic grids $G_a$ and the complete grid satisfy (App).

Next, we introduce an event on the tests $(T_{l,r})$ under which the change-point estimator $\hat{\tau}$ of Algorithm 1 performs well. In the following, we write $H_0$, the collection of all possible $(l, r) \in J_n$ such that there is no change in $[l-r+1, l+r-1]$, i.e. $\Gamma(\mathbb{P}_t)$ is constant on $[l-r, l+r)$. Equivalently, we have

$$(l, r) \in H_0 \quad \text{iff} \quad (l-r, l+r) \cap \{\tau_k, k = 1, \dots, K\} = \emptyset. \qquad (7)$$

For a collection $K^*$ and some elements of the grid $(\bar{\tau}_k, \bar{r}_k)$ satisfying (6), the Event $A(T, K^*, (\bar{\tau}_k, \bar{r}_k)_{k \in K^*})$ is defined as the conjunction of the two following properties: (i) (No false positive) $T_{l,r} = 0$ for all $(l, r) \in H_0 \cap G$ (ii) (Detection of significant change-points) for every $k \in K^*$, we have $T_{\bar{\tau}_k, \bar{r}_k} = 1$.

The first property states that $T$ performs no type I errors on the event $A(T, K^*, (\bar{\tau}_k, \bar{r}_k)_{k \in K^*})$, whereas the second property enforces that all the significant change-points are detected by the specific tests $T_{\bar{\tau}_k, \bar{r}_k}$.
$r_k^*$ as the minimum radius $r$ such that an inequality similar to (8) is satisfied for $r\Delta_k^2$, namely

$$r_k^* = \min \left\{ r \in \mathbb{R}^+ : r\Delta_k^2 \ge c_0\sigma^2 \left[ s_k \log \left( 1 + \frac{\sqrt{p}}{s_k} \sqrt{\log \left( \frac{n}{r\delta} \right)} \right) + \log \left( \frac{n}{r\delta} \right) \right] \right\} . \quad (9)$$

In the following, we introduce multi-scale tests for respectively dense and sparse change-points. For simplicity, we restrict our attention to the dyadic grid $\mathcal{G}_D = (\mathcal{R}, \mathcal{D})$ introduced in the previous section (see Equation (5)), the complete grid being used in the next section.

To apply Theorem 1, we will consider an event $\mathcal{A}(T, \mathcal{K}^*, (\bar{\tau}_k, \bar{r}_k)_{k \in \mathcal{K}^*})$ in the proof of Corollary 2 where the scale $\bar{r}_k \in \mathcal{R}$ is of the same order as $r_k^* \in \mathbb{R}^+$.

## 3.1 Dense change-points

We focus here on dense change-points for which $s_k$ is possibly as large as $p$. Given $\kappa > 0$, $\tau_k$ is a $\kappa$-dense high-energy change-point if

$$r_k \Delta_k^2 \ge \kappa \sigma^2 \left( \sqrt{p \log \left( \frac{n}{r_k \delta} \right)} + \log \left( \frac{n}{r_k \delta} \right) \right) . \quad (10)$$

The requirement (10) is analogous to (8) when $s_k \ge [p \log(n/(r_k \delta))]^{1/2}$. For any $\kappa$-dense high-energy change-point, we define $\bar{r}_k^{(d)} \in \mathcal{R}$ as the minimum radius $r \in \mathcal{R}$ such that an inequality of the same type as (10) is satisfied for $r\Delta_k^2$,

$$\bar{r}_k^{(d)} = \min \left\{ r \in \mathcal{R} : 8r\Delta_k^2 \ge \kappa\sigma^2 \left( \sqrt{p \log \left( \frac{n}{r\delta} \right)} + \log \left( \frac{n}{r\delta} \right) \right) \right\} .$$

Intuitively, $\bar{r}_k^{(d)}$ corresponds to the smallest scale such that $\tau_k$ is guaranteed to be detected. By definition, we have $4(\bar{r}_k^{(d)} - 1) \le r_k$. Let $\bar{\tau}_k^{(d)}$ be the best approximation of $\tau_k$ in the grid with scale $\bar{r}_k^{(d)}$. By definition of the dyadic grid, we have $|\bar{\tau}_k^{(d)} - \tau_k| \le \bar{r}_k^{(d)}/4$.

For any positive integers $r \in [1; n]$ and $l \in [r + 1, n + 1 - r]$, we define the statistic $\Psi_{l,r}^{(d)} := \|C_{l,r}\|^2 - p$. If $\theta$ is constant over $[l - r, l + r)$, then the expectation of $\Psi_{l,r}^{(d)}$ is zero. Recall that the rescaled CUSUM statistic $C_{l,r}$ depends on the noise level $\sigma$, and the statistic $\Psi_{l,r}^{(d)}$ therefore requires the knowledge of $\sigma$. To calibrate the corresponding test $T_{l,r}^{(d)}$ rejecting for large values of $\Psi_{l,r}^{(d)}$ we introduce

$$T_{l,r}^{(d)} := \mathbf{1} \left\{ \Psi_{l,r}^{(d)} > x_r^{(d)} \right\} ; \quad x_r^{(d)} := 4 \left( \sqrt{p \log \left( \frac{2n}{r\delta} \right)} + \log \left( \frac{2n}{r\delta} \right) \right) .$$

**Proposition 1.** There exists a universal constant $\kappa_d > 0$ and an event $\xi^{(d)}$ of probability larger than $1 - 2\delta$ such that (i) $T_{l,r}^{(d)} = 0$ for all $(l, r) \in \mathcal{H}_0 \cap \mathcal{G}_D$ and (ii) $T_{\bar{\tau}_k^{(d)}, \bar{r}_k^{(d)}}^{(d)} = 1$ for all $\kappa_d$-dense high-energy change-point $\tau_k$.

The above proposition ensures that, on the event $\xi^{(d)}$, the collection of tests $T_{l,r}^{(d)}$ detects all dense high-energy change-points at the scale $\bar{r}_k^{(d)}$ and makes no false positives on the dyadic grid $\mathcal{G}_D$. If we plugged this collection of tests into the general multiple change-point procedure, then Theorem 1 would entail that all $\kappa_d$-dense high-energy change-points are discovered and localized and that $\hat{\tau}$ does not detect any spurious change-point. In the next subsection, we introduce alternative tests that are tailored to sparse change-points and thereby allow to detect change-points that are not $\kappa_d$-dense high-energy but still satisfy the energy condition (8).
**Theorem 1.** The following holds for any grid $\mathcal{G}$, any local test statistic $T$, any non-negative integer $K$, any distribution with $K$ change-points, any $K^* \subset [K]$ and scales and locations $(\bar{\tau}_k, \bar{r}_k)_{k \in K^*}$ in $\mathcal{G}$ satisfying Assumption (6). Under the event $\mathcal{A}(T, K^*, (\bar{\tau}_k, \bar{r}_k)_{k \in K^*})$, the estimated change-point vector $\hat{\tau}$ returned by Algorithm 1 satisfies the two following properties

* **Significant change-points are detected:** for all $k \in K^*$, there exists $k' \le K$ such that $|\hat{\tau}_{k'} - \tau_k| \le \bar{r}_k - 1 < \frac{r_k}{4}$.

* **(NoSp):** No Spurious change-point is detected (1).

The first property states that so-called significant change-points $(\tau_k)_{k \in K^*}$ are detected by the generic algorithm at the right scale. The no-spurious property (1) guarantees that, around any true change-point $\tau_k$, the procedure estimates at most one single change-point $\hat{\tau}_l$. Importantly, the theorem does not make any assumption on the non-significant change-points.

In fact, change-points $\tau_k$ with $k \in [K] \setminus K^*$ may or may not be detected. In general, we can only conclude from Theorem 1 that $|K^*| \le \hat{K} \le K$ on the event $\mathcal{A}(T, K^*, (\bar{\tau}_k, \bar{r}_k)_{k \in K^*})$.

Theorem 1 is abstract but its main virtue is to translate multiple testing properties into change-point detection properties. For a specific problem such as multivariate mean change-point detection considered in the next section, the construction of a near optimal procedure boils down to introducing a collection of local test statistics, such that (a) change-points $\tau_k$ belong to $K^*$ under minimal conditions, (b) the scale $\bar{r}_k$ is the smallest possible, and (c) the event $\mathcal{A}(T, K^*, (\bar{\tau}_k, \bar{r}_k)_{k \in K^*})$ holds with high probability.

In the case where all the change-points are significant, the result of Theorem 1 can be reformulated as follows:

**Corollary 1.** The following holds for any grid $\mathcal{G}$, any local test statistic $T$, any non-negative integer $K$, any distribution with $K$ change-points, any sequence $(\bar{\tau}_k, \bar{r}_k)_{k=1,...,K}$ in $\mathcal{G}$ satisfying Assumption (6).

Under the event $\mathcal{A}(T, [K], (\bar{\tau}_k, \bar{r}_k))$, the estimated change-point vector $\hat{\tau}$ returned by Algorithm 1 satisfies $\hat{K} = K$ and,

$$|\hat{\tau}_k - \tau_k| < \bar{r}_k - 1 \le \frac{r_k}{4} \quad \text{for all } k = 1, \dots, K.$$

Let us respectively define the Hausdorff distance and the Wasserstein distance of two vectors $(u_1, \dots, u_K)$ and $(v_1, \dots, v_K)$ in $\mathbb{R}^K$ by

$d_H(u, v) = \max_{k=1,...,K} |u_k - v_k|$ and $d_W(u, v) = \sum_{k=1,...,K} |u_k - v_k|$. Then, Corollary 1 straightforwardly implies that, if $K^* = [K]$, then these two losses are bounded as follows

$$d_H(\hat{\tau}, \tau) \le \max_{k=1,...,K} (\bar{r}_k - 1) \quad \text{and} \quad d_W(\hat{\tau}, \tau) \le \sum_{k=1,...,K} (\bar{r}_k - 1).$$

As an alternative of Algorithm 1, one could use other bottom-up aggregating procedures. For instance, Algorithm 2 defined in Appendix A also satisfies Theorem 1. Although these two algorithms are closely related, Algorithm 1 is slightly
## 3.2 Sparse change-points

### 3.2.1 Energy condition

For a given $1 \le k \le K$, the change-point $\tau_k$ is a $\kappa$-sparse high-energy change-point if $s_k \le [p \log(n/(r_k \delta))]^{1/2}$ and

$$r_k \Delta_k^2 \ge \kappa \sigma^2 \left( s_k \log \left( \frac{p}{s_k^2} \log \left( \frac{n}{r_k \delta} \right) \right) + \log \left( \frac{n}{r_k \delta} \right) \right) . \quad (11)$$

If $\tau_k$ is a $\kappa$-sparse high-energy change-point, we define $\bar{r}_k^{(s)}$ as the minimum scale such that an inequality similar to (11) is satisfied:

$$\bar{r}_k^{(s)} = \min \left\{ r \in \mathcal{R} : 8r \Delta_k^2 \ge \kappa \sigma^2 \left( s_k \log \left( \frac{p}{s_k^2} \log \left( \frac{n}{r_k \delta} \right) \right) + \log \left( \frac{n}{r_k \delta} \right) \right) \right\} .$$

As in the dense case, we have $4(\bar{r}_k^{(s)} - 1) \le r_k$. Set $\bar{\tau}_k^{(s)}$ as the best approximation of $\tau_k$ in the grid $\mathcal{D}_{\bar{\tau}_k^{(s)}}$ at scale $\tau_k$. By definition of the dyadic grid, we have $|\bar{\tau}_k^{(s)} - \tau_k| \le \bar{r}_k^{(s)}/4$. We introduce below two statistics for handling this problem.

### 3.2.2 Berk-Jones Test

The Berk-Jones test [29] is a variation of the Higher-Criticism test originally introduced in [11] for signal detection. It has been previously studied in [6] for sparse segment detection. We decided to use the Berk-Jones test in this paper because of its intrinsic formulation in terms of the quantiles of a Bernoulli distribution, but the Higher-Criticism test would reach the same rates of detection within a constant factor. We use the notation $\mathbb{N}^*$ to denote the set of positive integers. Given $(l, r)$ in the grid $\mathcal{G}_D$, we first introduce $N_{x,l,r}$ as the number of coordinates of $\mathbf{C}_{l,r}$ that are larger than $x$ in absolute value.

$$N_{x,l,r} = \sum_{i=1}^{p} \mathbf{1}_{|\mathcal{C}_{l,r,i}| > x} \quad (12)$$

If $(l, r) \in \mathcal{H}_0$, then the rescaled CUSUM statistic follows a standard normal distribution and $N_{x,l,r}$ therefore follows a Binomial distribution with parameters $p$ and $2\bar{\Phi}(x)$. The Berk-Jones test amounts to rejecting the null, when at least one of the statistics $N_{x,l,r}$, for $x \in \mathbb{N}^*$, is significantly large. Next, we formalize what we mean by 'large'.

For any $u > 0$, any $q_0 \in [0, 1]$, and positive integer $p_0$, denote $\bar{Q}(u, p_0, q_0) = \mathbb{P}[\mathcal{B}(p_0, q_0) > u]$ the tail distribution function of a Binomial distribution with parameters $p_0$ and $q_0$. Given $\delta \in [0, 1]$, we then write $\bar{Q}^{-1}(\delta, p_0, q_0)$ for the corresponding quantile function,

$$\bar{Q}^{-1}(\delta, p_0, q_0) = \inf_{u} [\mathbb{P}[\mathcal{B}(p_0, q_0) > u] \le \delta] .$$

Given a scale $r \in \mathcal{R}$ and a positive integer $x$, we define the weights

$$\delta_{x,r}^{(\text{BJ})} = \frac{6\delta r}{\pi^2 x^2 |\mathcal{D}_r| n} . \quad (13)$$

This allows us to define the Berk-Jones statistic over $[l-r, l+r)$ as the test rejecting the null when at least one $N_{x,l,r}$ is large.

$$T_{l,r}^{(\text{BJ})} = \max_{x \in \mathbb{N}^*} \mathbf{1} \left\{ N_{x,l,r} > \bar{Q}^{-1}(\delta_{x,r}^{(\text{BJ})}, p, 2\bar{\Phi}(x)) \right\} . \quad (14)$$
more conservative than Algorithm 2 since it merges all detection intervals at a given resolution while Algorithm 2 only keeps one interval at a given resolution when multiple intervals intersect – the one with smallest index $t$. While the minimax properties of both methods are comparable – at least up to a multiple constant – the choice of aggregation method will have an influence in practice on the outcome: Algorithm 1 will be slightly more stable, detect less change-points, and provide wider confidence interval around them, while Algorithm 2 will be slightly more sensitive to smaller changes, i.e. detect smaller change-points, will be more precise, and somewhat less stable.

Theorem 1 ensures that, if $T_{\bar{\tau}_k, \bar{r}_k} = 1$ with $(\bar{\tau}_k, \bar{r}_k)$ satisfying Assumption (6), then the change-point $\tau_k$ is detected. Inspecting the proof of Theorem 1, one easily checks that Assumption (6) is minimal for Algorithm 1 (and also for Algorithm 2). Still, one may wonder whether any generic algorithm has to require that $4(\bar{r}_k - 1) < r_k$ to detect the change-points or if there exists a generic algorithm where the constant 4 in the above condition can be improved.

**Comparison with narrowest over threshold methods.** As mentioned in the introduction, other aggregation procedures have been proposed in the literature. In particular, the narrowest over threshold scheme proposed by [2] and later used in [24] is also closely related to the local segmentation algorithm of Chan and Chen [5]. A simple extension of these procedures for generic change-point problems and for a general collection of tests ($T_{l,r}$) would amount to modifying Algorithm 1 by selecting locations $l$ in $D_r$ such that $T_{l,r} = 1$ and $[l-r+1, l+r-1]$ does not intersect previously detected change-points, whereas we require in Algorithms 1 and 2, that $[l-r+1, l+r-1]$ does not intersect previously detected confidence intervals. In some way, the narrowest-over threshold scheme is therefore less conservative. Unfortunately, there is no generic result in the form of Theorem 1 for such procedures and, from informal arguments, we doubt that the corresponding procedure provably achieves (NoSp) under a control of the FWER of the tests. Inspecting the proof of Theorem 1 in [2] and Theorem 3 in [24] for univariate mean change-point problems, one observes that the chosen threshold is much larger than what is needed to control the FWER so that the theoretical threshold is certainly over-conservative – see step 5 of the proof of Theorem 1 in [2]. In contrast, Theorem 1 in [5] for univariate change-point problems is based on the minimal threshold, but the proof relies on the important assumption that the number $K$ of change-point remains bounded while $n$ goes to infinity. Besides, it is not clear how one could extend the arguments to more general settings.

### 3. Multivariate Gaussian change-point detection

We now turn to the multivariate change-point model introduced in Section 1.2. Throughout this section, we assume that the random vectors $\varepsilon_t$ are independently and identically distributed with $\varepsilon_t \sim N(0, \sigma^2 \mathbf{I}_p)$. Since we shall apply the general aggregation procedures introduced in the previous section, our main job here is to introduce a near-optimal testing procedure.
Equivalently, $T_{l,r}^{(\text{BJ})}$ is an aggregated test based on the statistics $N_{x,l,r}$ with weights $\delta_{x,r}^{(\text{BJ})}$. From the above remark and a union bound, we deduce that the probability that the collection of tests $\{T_{l,r}^{(\text{BJ})}, (l,r) \in \mathcal{G}_D\}$ rejects a least one false positive is at most $\delta$:

$$ \mathbb{P}\left[\max_{(l,r) \in \mathcal{H}_0 \cap \mathcal{G}_D} T_{l,r}^{(\text{BJ})} = 1\right] \le \sum_{r \in \mathcal{R}} \sum_{l \in \mathcal{D}_r} \sum_{x \in \mathbb{N}^*} \delta_{x,r}^{(\text{BJ})} \le \sum_{r \in \mathcal{R}} \sum_{l \in \mathcal{D}_r} \frac{\delta r}{|\mathcal{D}_r|n} \le \sum_{r \in \mathcal{R}} \frac{\delta r}{n} \le \delta, $$

where we recall that $(l,r) \in \mathcal{H}_0$ if and only if $\Theta$ is constant on $[l-r, l+r)$. Although one may think from the definition (14) that $T_{l,r}^{(\text{BJ})}$ involves an infinite number of $N_{x,l,r}$, this is not the case. Indeed, $N_{x,l,r}$ is a non-increasing function of $x$ whereas for all $x$ such that $2p\bar{\Phi}(x) \le \delta_{x,r}^{(\text{BJ})}$, we have $\bar{Q}^{-1}(\delta_{x,r}^{(\text{BJ})}, p, 2\bar{\Phi}(x)) = 0$. Writing $x_{0,r}$ the smallest $x$ such that $2p\bar{\Phi}(x) \le \delta_{x,r}^{(\text{BJ})}$ we derive

$$ T_{l,r}^{(\text{BJ})} = \max_{x=1,\dots,x_{0,r}} \mathbf{1}\left\{N_{x,l,r} > \bar{Q}^{-1}(\delta_{x,r}^{(\text{BJ})}, p, 2\bar{\Phi}(x))\right\}. $$

Since, for any $x > 0$, we have $\bar{\Phi}(x) \le e^{-x^2/2}$, one can deduce that $x_{0,r} \le c[\log(np/(r\delta))]^{1/2}$, for some numerical constant $c > 0$.

### 3.2.3 Partial norm statistics

The Berk-Jones test is able to detect change-points $\tau_k$ for which there exists $s$ such that the $s$ largest squared coordinates of $\mu_k - \mu_{k-1}$ are larger than $C(\log(ep/s^2) + \log(n/r_k)/s)$ with a large enough constant $C$. However, it may happen that $\tau_k$ satisfies the energy condition (8) and that the $s$ largest coordinates of $\mu_k - \mu_{k-1}$ are negligible compared to $\log(n/r_k)/s$, mainly because $s \mapsto 1/s$ is not summable. To solve this issue, we introduce a second sparse statistic based on the partial sums. Let

$$ \mathcal{Z} = \{1, 2, 2^2, \dots, 2^{\lceil \log_2(p) \rceil}\} $$

denote the dyadic set. Only the sparsities $s \in \mathcal{Z}$ will be analysed by the partial norm statistic. For any $(l,r)$ in the grid $\mathcal{G}_D$, we respectively write $\mathbf{C}_{l,r,(1)}$, $\mathbf{C}_{l,r,(2)}$, \dots the reordered entries of $\mathbf{C}_{l,r}$ by decreasing absolute value, that is $|\mathbf{C}_{l,r,(1)}| \ge \dots \ge |\mathbf{C}_{l,r,(p)}|$. Then, for $s \in \mathcal{Z}$, we define the partial CUSUM norm by

$$ \Psi_{l,r,s}^{(\text{p})} = \sum_{i=1}^{s} (\mathbf{C}_{l,r,(i)})^2 . \qquad (15) $$

Then, we define the test $T_{l,r}^{(\text{p})}$ rejecting the null when at least one of the partial norms is large

$$ x_{r,s}^{(\text{p})} := x_{r,s}^{(\text{p})}(\delta) = 4s \log\left(\frac{2ep}{s}\right) + 4 \log\left(\frac{n}{r\delta}\right); \quad T_{l,r}^{(\text{p})} = \max_{s \in \mathcal{Z}} \mathbf{1}\left\{\Psi_{l,r,s}^{(\text{p})} > x_{r,s}^{(\text{p})}\right\}. $$

Finally, we define the sparse test by aggregating both the Berk-Jones test and the partial norm test. For any $(l,r) \in \mathcal{G}_D$, let $T_{l,r}^{(\text{s})} = T_{l,r}^{(\text{p})} \lor T_{l,r}^{(\text{BJ})}$. The next proposition controls the error of this collection of tests.

**Proposition 2.** There exists a universal constant $\kappa_s > 0$ and an event $\xi^{(s)}$ of probability larger than $1 - 4\delta$ such that (i) $T_{l,r}^{(\text{s})} = 0$ for all $(l,r) \in \mathcal{H}_0 \cap \mathcal{G}_D$ and (ii) $T_{\bar{\tau}_k^{(s)}, \bar{\tau}_k^{(s)}}^{(\text{s})} = 1$ for all $\kappa_s$-sparse high-energy change-point $\tau_k$.
Fix some quantity $\delta \in (0, 1)$. At the end of the section, $1-\delta$ will correspond to the probability of the event $\mathcal{A}(T, \mathcal{K}^*, (\bar{\tau}_k, \bar{r}_k)_{k \in \mathcal{K}^*})$ introduced in the previous section. Alternatively, one may interpret $\delta$ as an upper bound of the desired probability that the change-point detection procedure detects a spurious change-points. Recall that, for a change-point $\tau_k$, $s_k$ stands for the sparsity of the difference $\mu_{k+1} - \mu_k$. The energy of a given change-point $\tau_k$ is $c_0$-high if

$$r_k \Delta_k^2 \geq c_0 \sigma^2 \left[ s_k \log \left( 1 + \frac{\sqrt{p}}{s_k} \sqrt{\log \left( \frac{n}{r_k \delta} \right)} \right) + \log \left( \frac{n}{r_k \delta} \right) \right], \quad (8)$$

for some universal constant $c_0$ to be defined later. We show in this section that when $c_0$ is large enough, all high-energy change-points can be detected. Conversely, it is established in Section 5 that Condition (8) is (up to a multiplicative constant) optimal for detecting change-points and cannot be weakened.

Let us now discuss the different regimes contained in Equation (8). In what follows, define

$$\psi_{n,r,s}^{(g)} := s \log \left( 1 + \frac{\sqrt{p}}{s} \sqrt{\gamma_r} \right) + \gamma_r; \quad \gamma_r := \log \left( \frac{n}{r\delta} \right),$$

in order to alleviate notations. If $\gamma_r \ge p/2$, then $\psi_{n,r,s}^{(g)} \asymp \gamma_r$ where $u \asymp v$ means that for two positive numerical constants $c_1$ and $c_2$, one has $c_1v \le u \le c_2v$. This corresponds to the minimal energy condition for detection in the univariate case, i.e. when $p=1$; see [38]. The condition $\gamma_r \ge p/2$ occurs when $p$ is rather small and the scale $r$ is much smaller than $n$. If $\gamma_r \le p/2$, then

$$\psi_{n,r,s}^{(g)} \asymp \begin{cases} \gamma_r & \text{if } s \le \frac{\gamma_r}{\log(p)-\log(\gamma_r)} \\ s \log\left(2 \frac{p}{s^2} \gamma_r\right) & \text{if } \frac{\gamma_r}{\log(p)-\log(\gamma_r)} < s < \sqrt{p\gamma_r} \\ \sqrt{p\gamma_r} & \text{if } s \ge \sqrt{p\gamma_r}. \end{cases}$$

We define $\mathcal{K}^* \subset [K]$ as the subset of indices such that $\tau_k$ satisfies (8). For any $k \in \mathcal{K}^*$, we define $r_k^*$ as the minimum radius $r$ such that an inequality similar to (8) is satisfied for $r\Delta_k^2$, namely

$$r_k^* = \min \left\{ r \in \mathbb{R}^+ : r \Delta_k^2 \ge c_0 \sigma^2 \left[ s_k \log \left( 1 + \frac{\sqrt{p}}{s_k} \sqrt{\log \left( \frac{n}{r\delta} \right)} \right) + \log \left( \frac{n}{r\delta} \right) \right] \right\}. \quad (9)$$

In the following, we introduce multi-scale tests for respectively dense and sparse change-points. For simplicity, we restrict our attention to the dyadic grid $\mathcal{G}_D = (\mathcal{R}, \mathcal{D})$ introduced in the previous section (see Equation (5)), the complete grid being used in the next section.

To apply Theorem 1, we will consider an event $\mathcal{A}(T, \mathcal{K}^*, (\bar{\tau}_k, \bar{r}_k)_{k \in \mathcal{K}^*})$ in the proof of Corollary 2 where the scale $\bar{r}_k \in \mathcal{R}$ is of the same order as $r_k^* \in \mathbb{R}^+$.
Here we introduced two different statistics for the same sparse regime $s_k \le [p \log(n/(r_k \delta))]^{1/2}$ - the Berk-Jones statistic and the partial sums statistic - mainly to solve a problem of integrability. We made this choice for the sake of simplicity, but we could have used a single test, as presented in [28]

$$ \Psi_{x,l,r}^{(\text{LGS})} = \sum_{i=1}^{p} (C_{l,r,i}^2 - \mathbb{E}[Z|Z \ge x]) \mathbf{1}\{C_{l,r,i}^2 \ge x\}, $$

where $Z$ follows a standard normal distribution $N(0, 1)$. This statistic leads to the same type of result as the Berk-Jones statistic when enough coordinates $\mu_k - \mu_{k-1}$ are large in absolute value, and it is comparable to the partial sums statistic when its threshold $x$ becomes low enough.

### 3.3 Consequences

To conclude this section, it suffices to observe that, for $c_0$ in (8), any $c_0$-high-energy change-point $\tau_k$ in the sense of (8) is either a $\frac{c_0}{2}$-dense or a $\frac{c_0}{2}$-sparse high-energy change-point. Hence, upon defining the test $T_{l,r} = T_{l,r}^{(d)} \lor T_{l,r}^{(s)}$ for $(l, r) \in \mathcal{G}_D$, we consider the change-point procedure $\hat{\tau}$ defined in Algorithm 1. Gathering Theorem 1 with Proposition 1 and Proposition 2, we obtain the following.

**Corollary 2.** There exists a universal constant $c_0 > 0$ such that, with probability higher than $1-6\delta$, the estimator $\hat{\tau}$ satisfies (NoSp) and detects all $c_0$-high-energy change-points (as defined in (8)) $\tau_k$ in the sense

$$ d_{H,1}(\hat{\tau}, \tau_k) < \frac{r_k^*}{2} \le \frac{r_k}{2}, $$

where $r_k^*$ is defined in (9).

If the change-points are of high-energy, that is $K^* = [K]$, then Corollary 2 can be reformulated as follows:

**Corollary 3.** Assume that for all $k = 1, \dots, K$, $\tau_k$ is a $c_0$-high-energy change-point (see (8)) where $c_0$ is the same as in Corollary 2. Then, with probability higher than $1-6\delta$, the estimator $\hat{\tau}$ satisfies $\hat{K} = K$ and

$$ |\hat{\tau}_k - \tau_k| < \frac{r_k^*}{2} \le \frac{r_k}{2}, \quad \text{for all } k = 1, \dots, K. $$

In particular, one can respectively bound the Hausdorff and the Wasserstein losses, with probability higher than $1-6\delta$ by

$$ d_H(\hat{\tau}, \tau) \le \max_{k=1,\dots,K} \frac{r_k^*}{2} \quad \text{and} \quad d_W(\hat{\tau}, \tau) \le \sum_{k=1,\dots,K} \frac{r_k^*}{2}. \qquad (16) $$

In Section 5, we establish that the Condition (8) is (up to a multiplicative constant) unimprovable and corresponds to the detection threshold for multivariate change-points.

Corollary 3 can be compared to the result of [46] on multivariate change-point detection in the multiple change-point setting. Using a method based on the CUSUM statistic and assuming that there are only high-energy change-points, the authors also obtain an upper bound on the energy necessary to detect the change-points. However, this result does not adapt to $r_k, \Delta_k, s_k$, and the detection rate is suboptimal in many regimes. Writing $r = \min_{k=1,\dots,K} r_k$, $\Delta = \min_{k=1,\dots,K} \Delta_k$ and $s = \max_{k=1,\dots,K} s_k$, Theorem 5 of [46] requires two conditions of the type $r\Delta^2 \ge c(\frac{n}{r})^4 \log(np)$ and $r\Delta^2 \ge c s \frac{n}{r} \log(np)$. This detection rate is therefore suboptimal by a polynomial factor in $n/r$ when
### 3.1. Dense change-points

We focus here on dense change-points for which $s_k$ is possibly as large as $p$. Given $\kappa > 0$, $\tau_k$ is a $\kappa$-dense high-energy change-point if

$$r_k \Delta_k^2 \geq \kappa \sigma^2 \left( \sqrt{p \log \left( \frac{n}{r_k \delta} \right)} + \log \left( \frac{n}{r_k \delta} \right) \right) . \quad (10)$$

The requirement (10) is analogous to (8) when $s_k \geq [p \log(n/(r_k \delta))]^{1/2}$. For any $\kappa$-dense high-energy change-point, we define $\bar{r}_k^{(d)} \in \mathcal{R}$ as the minimum radius $r \in \mathcal{R}$ such that an inequality of the same type as (10) is satisfied for $r \Delta_k^2$,

$$\bar{r}_k^{(d)} = \min \left\{ r \in \mathcal{R} : 8r \Delta_k^2 \geq \kappa \sigma^2 \left( \sqrt{p \log \left( \frac{n}{r \delta} \right)} + \log \left( \frac{n}{r \delta} \right) \right) \right\} .$$

Intuitively, $\bar{r}_k^{(d)}$ corresponds to the smallest scale such that $\tau_k$ is guaranteed to be detected. By definition, we have $4(\bar{r}_k^{(d)} - 1) \leq r_k$. Let $\bar{\tau}_k^{(d)}$ be the best approximation of $\tau_k$ in the grid with scale $\bar{r}_k^{(d)}$. By definition of the dyadic grid, we have $|\bar{\tau}_k^{(d)} - \tau_k| \leq \bar{r}_k^{(d)}/4$.

For any positive integers $r \in [1; n]$ and $l \in [r + 1, n + 1 - r]$, we define the statistic $\Psi_{l,r}^{(d)} := \|C_{l,r}\|^2 - p$. If $\theta$ is constant over $[l - r, l + r)$, then the expectation of $\Psi_{l,r}^{(d)}$ is zero. Recall that the rescaled CUSUM statistic $C_{l,r}$ depends on the noise level $\sigma$, and the statistic $\Psi_{l,r}^{(d)}$ therefore requires the knowledge of $\sigma$. To calibrate the corresponding test $T_{l,r}^{(d)}$ rejecting for large values of $\Psi_{l,r}^{(d)}$ we introduce

$$T_{l,r}^{(d)} := \mathbf{1} \{ \Psi_{l,r}^{(d)} > x_r^{(d)} \} ; \quad x_r^{(d)} := 4 \left( \sqrt{p \log \left( \frac{2n}{r\delta} \right)} + \log \left( \frac{2n}{r\delta} \right) \right) .$$

**Proposition 1.** There exists a universal constant $\kappa_d > 0$ and an event $\xi^{(d)}$ of probability larger than $1 - 2\delta$ such that (i) $T_{l,r}^{(d)} = 0$ for all $(l, r) \in \mathcal{H}_0 \cap \mathcal{G}_D$ and (ii) $T_{\bar{\tau}_k^{(d)}, \bar{r}_k^{(d)}}^{(d)} = 1$ for all $\kappa_d$-dense high-energy change-point $\tau_k$.

The above proposition ensures that, on the event $\xi^{(d)}$, the collection of tests $T_{l,r}^{(d)}$ detects all dense high-energy change-points at the scale $\bar{r}_k^{(d)}$ and makes no false positives on the dyadic grid $\mathcal{G}_D$. If we plugged this collection of tests into the general multiple change-point procedure, then Theorem 1 would entail that all $\kappa_d$-dense high-energy change-points are discovered and localized and that $\hat{\tau}$ does not detect any spurious change-point. In the next subsection, we introduce alternative tests that are tailored to sparse change-points and thereby allow to detect change-points that are not $\kappa_d$-dense high-energy but still satisfy the energy condition (8).
r is of smaller order than n, and by a logarithmic factor $log(np)$ instead of $log(1 + \sqrt{p}/s \log(n/r)) + \frac{1}{s} \log(n/r)$ when r is of order n. Closer to our results, [5] have introduced another bottom-up procedure in the very specific asymptotic setting $n = e^{p\zeta}$ for $\zeta \in (0, 1)$ with a fixed K number of change-points. Assuming that, for each change-point, at least s coordinates of $\mu_{k+1} - \mu_{k+1}$ are larger than $\zeta$ in absolute value, [5] establish that their procedure provably detects the change-points as long as

$$rs\zeta^2 \ge c \begin{cases} \sqrt{p \log(n)} & \text{if } s \ge 0.5\sqrt{p \log(n)} \\ s \log\left(\frac{p}{s^2} \log(n)\right) & \text{if } s \le 0.5\sqrt{p \log(n)} \end{cases}.$$

In their specific asymptotic regime and when all non-zero coordinates are of the same order, and all the change-points have a similar length $r_k$, their result is similar to ours up to the logarithmic terms. Indeed, for equispaced change-points, our logarithmic term $log(n/r_k) = \log(K)$ is much smaller than $log(n)$. Besides, their result does not handle the presence of low-energy change-points and does not hold beyond the asymptotic regime $n = e^{p\zeta}$. In contrast, our condition (8) for high-energy change-points entails that the detection conditions are qualitatively different for other scalings in n and p. On the technical side, our condition (8) is of $l_2$ type whereas that in [5] is of minimal non-zero type. Recovering the tight $l_2$ conditions turns out to be much more challenging as we need to handle situations where some coordinates have different orders of magnitude. This is the main reason why we need to resort to a combination of the Berk-Jones and the partial-norm statistics.

**Comparison to one change-point problem.** When one knows that $K \le 1$ (at most one change-point), then [28] proved that it is possible to detect $\tau_1$ if and only if $r_1\Delta_1^2 \ge c\sigma^2[s_1 \log(1 + \frac{1}{s_1}\sqrt{p \log \log 8n}) + \log \log 8n]$. As in the univariate setting, the problem with only one change-point is simpler than for general $K \ge 2$. As for our procedure, Liu et al. [28] rely on statistics based on the CUSUM - a chi square statistics in the dense case and a thresholded sum of squared coordinates in the sparse case - to detect and localize $\tau_1$. It turns out that the detection procedure of [28] adapts to distance $r_1 = \max(\tau_1 - 1, n + 1 - \tau_1)$ the boundary, and one could refine their result by stating that $\tau_1$ is detectable if and only if $r_1\Delta_1^2 \ge c\sigma^2[s_1 \log(1 + \frac{1}{s_1}\sqrt{p \log \log(2n/r_1)}) + \log \log(2n/r_1)]$ which is more smaller when $r_1$ is of the order of n. This refined result is in the same spirit as our bounds for multiple change-point, but the rate is faster because one obtains $\log \log(n/r_1)$ - instead of $\log(n/r_k)$ in our case. The reason for this faster rate is due to the relative simplicity of the problem with only one change-point. Indeed, in single change-point detection, there is no need to look for change-points at all positions and scale at the same time, since scale and positions are related. This implies that it is possible to attain faster rates than in multiple change-point detection. The comparison between single and multiple change-point detection is thoroughly done in [38] for univariate models.

**Computational Cost.** The cost of the tests $T_{l,r}^{(d)}$ in the dense regime is $O(rp)$. The computation of the partial norm statistic requires to sort the coordinates $C_{l,r,i}$ of the CUSUM statistic, which takes $O(p(r+\log(p)))$ operations. Since only the thresholds $x \le c \log(np/(r\delta))^{1/2}$ are needed to compute the Berk-Jones statistic, it holds that, for $\delta \ge (np)^{-c}$ with a numerical constant $c > 0$, the computational cost of the Berk-Jones statistic is $O(p(r+\log(np)))$. Thus, for each $(l, r)$, the overall computational cost of the test $T_{l,r} = T_{l,r}^{(d)} \lor T_{l,r}^{(s)}$ is $\Lambda = O(p(r+\log(np)))$, and the computational cost of the whole change-point detection procedure on the dyadic grid is $O(np \log(np))$.
## 3.2. Sparse change-points

### 3.2.1. Energy condition

For a given $1 \le k \le K$, the change-point $\tau_k$ is a $\kappa$-sparse high-energy change-point if $s_k \le [p \log(n/(r_k \delta))]^{1/2}$ and

$$r_k \Delta_k^2 \ge \kappa \sigma^2 \left( s_k \log \left( \frac{p}{s_k^2} \log \left( \frac{n}{r_k \delta} \right) \right) + \log \left( \frac{n}{r_k \delta} \right) \right) . \quad (11)$$

If $\tau_k$ is a $\kappa$-sparse high-energy change-point, we define $\bar{r}_k^{(s)}$ as the minimum scale such that an inequality similar to (11) is satisfied:

$$\bar{r}_k^{(s)} = \min \left\{ r \in \mathcal{R} : 8r \Delta_k^2 \ge \kappa \sigma^2 \left( s_k \log \left( \frac{p}{s_k^2} \log \left( \frac{n}{r_k \delta} \right) \right) + \log \left( \frac{n}{r_k \delta} \right) \right) \right\} .$$

As in the dense case, we have $4(\bar{r}_k^{(s)} - 1) \le r_k$. Set $\bar{\tau}_k^{(s)}$ as the best approximation of $\tau_k$ in the grid $\mathcal{D}_{\bar{\tau}_k^{(s)}}$ at scale $\tau_k$. By definition of the dyadic grid, we have $|\bar{\tau}_k^{(s)} - \tau_k| \le \bar{r}_k^{(s)}/4$. We introduce below two statistics for handling this problem.

### 3.2.2. Berk-Jones test

The Berk-Jones test [29] is a variation of the Higher-Criticism test originally introduced in [11] for signal detection. It has been previously studied in [6] for sparse segment detection. We decided to use the Berk-Jones test in this paper because of its intrinsic formulation in terms of the quantiles of a Bernoulli distribution, but the Higher-Criticism test would reach the same rates of detection within a constant factor. We use the notation $\mathbb{N}^*$ to denote the set of positive integers. Given $(l, r)$ in the grid $\mathcal{G}_D$, we first introduce $N_{x,l,r}$ as the number of coordinates of $\mathbf{C}_{l,r}$ that are larger than $x$ in absolute value.

$$N_{x,l,r} = \sum_{i=1}^{p} \mathbf{1}_{|\mathbf{C}_{l,r,i}| > x} \quad (12)$$

If $(l, r) \in \mathcal{H}_0$, then the rescaled CUSUM statistic follows a standard normal distribution and $N_{x,l,r}$ therefore follows a Binomial distribution with parameters $p$ and $2\bar{\Phi}(x)$. The Berk-Jones test amounts to rejecting the null, when at least one of the statistics $N_{x,l,r}$, for $x \in \mathbb{N}^*$, is significantly large. Next, we formalize what we mean by 'large'.

For any $u > 0$, any $q_0 \in [0, 1]$, and positive integer $p_0$, denote $\bar{Q}(u, p_0, q_0) = \mathbb{P}[\mathcal{B}(p_0, q_0) > u]$ the tail distribution function of a Binomial distribution with parameters $p_0$ and $q_0$. Given $\delta \in [0, 1]$, we then write $\bar{Q}^{-1}(\delta, p_0, q_0)$ for the corresponding quantile function,

$$\bar{Q}^{-1}(\delta, p_0, q_0) = \inf_u [\mathbb{P}[\mathcal{B}(p_0, q_0) > u] \le \delta] .$$
# 4 Multi-scale change-point detection with sub-Gaussian noise

We now turn to the more general case of sub-Gaussian distributions [37]. Given a random variable $Z$, define its $\psi_2$-norm by $\|Z\|_{\psi_2} = \inf\{x > 0, \mathbb{E}[\exp(Z^2/x^2)] \le 2\}$. Given $L > 0$, a mean zero real random variable is said to be $L$-sub-Gaussian if $\|Z\|_{\psi_2} \le L$. This implies in particular that, for all $x \ge 0$, one has $\mathbb{P}(|Z| \ge x) \le 2\exp(-x^2/L^2)$. Throughout this section, we assume that, for $t = 1, \dots, n$, the random vectors $\varepsilon_t$ are independent, have independent $L$-sub-Gaussian components $\varepsilon_{t,i}$, for $i = 1, \dots, p$ with variance $\sigma^2$. As in the previous section, we apply the general aggregation procedures introduced in Section 2. As a consequence, our main task boils down to introducing a near-optimal multiple testing procedure indexed by a grid for detecting the existence of a change-point. Here, we shall rely on the complete grid $\mathcal{G}_F = J_n = \{(l, r) : r = 1, \dots, \lfloor \frac{n}{2} \rfloor \text{ and } l = r + 1, \dots, n - r\}$ whose size is quadratic with respect to $n$. All the results presented in this section are still valid (but with different numerical constants) if we keep the dyadic grid $\mathcal{G}_D$ as in the previous section. Here, we use the complete grid as a proof of concept that one can rely on the full collection of possible segments without deteriorating the rates. Still, controlling the behavior of the procedure on the complete grid is technically more involved and requires chaining arguments. A detailed comparison between the complete and dyadic grids is made in Section 7.

In order to emphasize the common points with the previous section, we use the same notation $\mathcal{K}^*$ for the collection of high-energy change-points<fnref n="1" />, $\bar{r}_k$ for the scales associated to the $k$-th change-points<fnref n="2" />, $\Psi$ for the statistics, $T$ for the test and $x$ for the thresholds although these quantities are slightly changed to cope with the sub-Gaussian tail distribution. We follow the same scheme as for the Gaussian case and first introduce multi-scale tests for dense change-points before turning to sparse change-points. As in the previous section, we consider some $\delta \in (0, 1)$ corresponding to the type I error probability.

## 4.1 Dense change-points with sub-Gaussian noise

Recall that, for a change-point $\tau_k$, $s_k$ stands for the sparsity of the difference $\mu_{k+1} - \mu_k$. We focus here on dense change-points for which $s_k$ is possibly as large as $p$. Given $\kappa > 0$, $\tau_k$ is a $\kappa$-dense high-energy change-point if

$$r_k \Delta_k^2 \ge \kappa L^2 \left( \sqrt{p \log \left( \frac{n}{r_k \delta} \right)} + \log \left( \frac{n}{r_k \delta} \right) \right) . \quad (17)$$

This condition is very similar to its counterpart (10) for Gaussian noise. Still, we introduce it here for the sake of completeness. For $k \in [K]$ such that $\tau_k$ is a $\kappa$-dense high-energy change-point, we define $\bar{r}_k^{(d)}$ as the minimum length such that an inequality similar to (17) is satisfied :

$$\bar{r}_k^{(d)} = \min \left\{ r \in \mathbb{N}^* : 4r \Delta_k^2 \ge \kappa L^2 \left( \sqrt{p \log \left( \frac{n}{r\delta} \right)} + \log \left( \frac{n}{r\delta} \right) \right) \right\} .$$

As in the Gaussian case in Section 3, $\bar{r}_k^{(d)}$ corresponds to the smallest scale such that $\tau_k$ is guaranteed to be detected. For any $\kappa$-dense high-energy change-point, it holds that $4(\bar{r}_k^{(d)} - 1) < r_k$. For any positive integers $(l, r) \in \mathcal{G}_F$, we consider the same CUSUM-based statistic $\Psi_{l,r}^{(d)} := \|\mathbf{C}_{l,r}\|^2 - p$ as

<footnote n="1">¹See Equation (20) as the energy condition is slightly different in the sub-Gaussian setting.</footnote>

<footnote n="2">²Re-defined in Equation (21).</footnote>
Given a scale $r \in \mathbb{R}$ and a positive integer $x$, we define the weights

$$ \delta_{x,r}^{(\text{BJ})} = \frac{6\delta r}{\pi^2 x^2 |\mathcal{D}_r| n} . \qquad (13) $$

This allows us to define the Berk-Jones statistic over $[l - r, l + r)$ as the test rejecting the null when at least one $N_{x,l,r}$ is large.

$$ T_{l,r}^{(\text{BJ})} = \max_{x \in \mathbb{N}^*} \mathbf{1}\left\{N_{x,l,r} > \overline{Q}^{-1}(\delta_{x,r}^{(\text{BJ})}, p, 2\overline{\Phi}(x))\right\}. \qquad (14) $$

Equivalently, $T_{l,r}^{(\text{BJ})}$ is an aggregated test based on the statistics $N_{x,l,r}$ with weights $\delta_{x,r}^{(\text{BJ})}$. From the above remark and a union bound, we deduce that the probability that the collection of tests $\{T_{l,r}^{(\text{BJ})}, (l,r) \in \mathcal{G}_D\}$ rejects a least one false positive is at most $\delta$:

$$ \mathbb{P}\left[\max_{(l,r) \in \mathcal{H}_0 \cap \mathcal{G}_D} T_{l,r}^{(\text{BJ})} = 1\right] \le \sum_{r \in \mathbb{R}} \sum_{l \in \mathcal{D}_r} \sum_{x \in \mathbb{N}^*} \delta_{x,r}^{(\text{BJ})} \le \sum_{r \in \mathbb{R}} \sum_{l \in \mathcal{D}_r} \frac{\delta r}{|\mathcal{D}_r|n} \le \sum_{r \in \mathbb{R}} \frac{\delta r}{n} \le \delta, $$

where we recall that $(l,r) \in \mathcal{H}_0$ if and only if $\Theta$ is constant on $[l - r, l + r)$. Although one may think from the definition (14) that $T_{l,r}^{(\text{BJ})}$ involves an infinite number of $N_{x,l,r}$, this is not the case. Indeed, $N_{x,l,r}$ is a non-increasing function of $x$ whereas for all $x$ such that $2p\overline{\Phi}(x) \le \delta_{x,r}^{(\text{BJ})}$, we have $\overline{Q}^{-1}(\delta_{x,r}^{(\text{BJ})}, p, 2\overline{\Phi}(x)) = 0$. Writing $x_{0,r}$ the smallest $x$ such that $2p\overline{\Phi}(x) \le \delta_{x,r}^{(\text{BJ})}$ we derive

$$ T_{l,r}^{(\text{BJ})} = \max_{x=1, \dots, x_{0,r}} \mathbf{1}\left\{N_{x,l,r} > \overline{Q}^{-1}(\delta_{x,r}^{(\text{BJ})}, p, 2\overline{\Phi}(x))\right\}. $$

Since, for any $x > 0$, we have $\overline{\Phi}(x) \le e^{-x^2/2}$, one can deduce that $x_{0,r} \le c[\log(np/(r\delta))]^{1/2}$, for some numerical constant $c > 0$.

### 3.2.3. Partial norm statistics

The Berk-Jones test is able to detect change-points $\tau_k$ for which there exists $s$ such that the $s$ largest squared coordinates of $\mu_k - \mu_{k-1}$ are larger than $C(\log(ep/s^2) + \log(n/r_k)/s)$ with a large enough constant $C$. However, it may happen that $\tau_k$ satisfies the energy condition (8) and that the $s$ largest coordinates of $\mu_k - \mu_{k-1}$ are negligible compared to $\log(n/r_k)/s$, mainly because $s \mapsto 1/s$ is not summable. To solve this issue, we introduce a second sparse statistic based on the partial sums. Let

$$ \mathcal{Z} = \{1, 2, 2^2, \dots, 2^{\lceil \log_2(p) \rceil}\} $$

denote the dyadic set. Only the sparsities $s \in \mathcal{Z}$ will be analysed by the partial norm statistic. For any $(l,r)$ in the grid $\mathcal{G}_D$, we respectively write $\mathbf{C}_{l,r,(1)}$, $\mathbf{C}_{l,r,(2)}$, . . . the reordered entries of $\mathbf{C}_{l,r}$ by decreasing absolute value, that is
for Gaussian noise. Let $\bar{c}_{\text{thresh}}^{(d)} > 0$ be a tuning parameter to be discussed later. To calibrate the corresponding multiple test procedures ($T_{l,r}^{(d)}$) with $(l,r) \in \mathcal{G}_F$ rejecting for large values of $\Psi_{l,r}^{(d)}$ we introduce

$$T_{l,r}^{(d)} := \mathbf{1} \{ \Psi_{l,r}^{(d)} > x_r^{(d)} \} ; \quad x_r^{(d)} = \bar{c}_{\text{thresh}}^{(d)} \frac{L^2}{\sigma^2} \left( \sqrt{p \log \left( \frac{n}{r\delta} \right)} + \log \left( \frac{n}{r\delta} \right) \right) .$$

**Proposition 3.** There exists a numerical constant $\bar{c}_{\text{thresh}}^{(d)} > 0$ such that the following holds for any $\kappa_d > 32\bar{c}_{\text{thresh}}^{(d)}$. With probability higher than $1 - \delta$, one has (i) $T_{l,r}^{(d)} = 0$ for all $(l,r) \in \mathcal{G}_F \cap \mathcal{H}_0$ and (ii) $T_{\tau_k, \bar{r}_k^{(d)}}^{(d)} = 1$ for all $\kappa_d$-dense high-energy change-points $\tau_k$.

In comparison to Proposition 1 in the previous section, there are two differences. First, we need to cope with sub-Gaussian distribution by applying the Hanson-Wright inequality. Most importantly, the grid $\mathcal{G}_F$ is much larger than $\mathcal{G}_D$ so that we cannot simply consider each test $T_{l,r}$ separately and simply apply a union bound as in the previous section. To handle the dependencies between the statistics $\Psi_{l,r}^{(d)}$, we have to apply a chaining argument. In fact, the thresholds $x_r^{(d)}$ are similar to their counterpart in the previous section, whereas the number $|\mathcal{G}_F|$ of tests is now proportional to $n^2$. In principle, the benefit of using the full grid $\mathcal{G}_F$ is that $(\tau_k, \bar{r}_k^{(d)})$ belongs to $\mathcal{G}_F$ so that we can consider the CUSUM statistic based on a segment $[\tau_k - \bar{r}_k^{(d)}, \tau_k + \bar{r}_k^{(d)}]$ centered around the change-point $\tau_k$. In contrast, $(\tau_k, \bar{r}_k^{(d)})$ does not necessarily belong to the dyadic grid $\mathcal{G}_D$ and we needed to consider its best approximation $(\bar{\tau}_k^{(d)}, \bar{r}_k^{(d)})$. The segment $[\bar{\tau}_k^{(d)} - \bar{r}_k^{(d)}, \bar{\tau}_k^{(d)} + \bar{r}_k^{(d)}]$ is therefore not centered on $\tau_k$ and the corresponding statistic $\Psi_{\bar{\tau}_k^{(d)}, \bar{r}_k^{(d)}}^{(d)}$ is in expectation smaller than $\Psi_{\tau_k, \bar{r}_k^{(d)}}^{(d)}$. In summary, both the collections of dense tests $\Psi_{l,r}^{(d)}$ on $\mathcal{G}_D$ and $\mathcal{G}_F$ are able to detect change-points whose energy is, up to some multiplicative constants, higher than $L^2[[p \log(\frac{n}{r_k \delta})]^{1/2} + \log(\frac{n}{r_k \delta})]]$.

## 4.2 Sparse change-points with sub-Gaussian noise

Unlike in the Gaussian case, we do not know the exact distribution of the noise. As a consequence, the Berk-Jones test and more generally higher-criticism type tests cannot be applied to this setting. This is why we only rely on the partial norm statistic. Recall that $\mathcal{Z} = \{1, 2, 2^2, \dots, 2^{\log_2(p)}\}$ stands for a dyadic set of sparsities. For $(l,r) \in \mathcal{G}_F$ and $s \in \mathcal{Z}$, we also recall that the partial CUSUM norm is defined as $\Psi_{l,r,s}^{(p)} = \sum_{i=1}^s (C_{l,r,(i)})^2$. Then, for any $(l,r) \in \mathcal{G}_F$, the test $T_{l,r}^{(p)}$ rejects the null when at least one of the partial norms is large

$$x_{r,s}^{(p)} = s + \bar{c}_{\text{thresh}}^{(p)} \frac{L^2}{\sigma^2} \left[ s \log \left( \frac{2ep}{s} \right) + \log \left( \frac{n}{r\delta} \right) \right]; \quad T_{l,r}^{(p)} = \max_{s \in \mathcal{Z}} \mathbf{1} \{ \Psi_{l,r,s}^{(p)} > x_{r,s}^{(p)} \},$$

where $\bar{c}_{\text{thresh}}^{(p)}$ is a tuning parameter in Proposition 4 below. The partial norm test alone is not able to detect sparse high-energy change-points in the sense of (11) and we need to introduce a stronger condition on the energy. Given $\kappa > 0$, a change-point $\tau_k$ is a $\kappa$-sparse high-energy change-point in the sub-Gaussian setting if $s_k \le [p \log(\frac{n}{r_k \delta})]^{1/2}$ and

$$r_k \Delta_k^2 \ge \kappa L^2 \left[ s_k \log \left( \frac{ep}{s_k} \right) + \log \left( \frac{n}{r_k \delta} \right) \right]. \quad (18)$$

Both Conditions (11) and (18) are compared at the end of the subsection. For a $\kappa$-sparse high-energy change-point $\tau_k$, we define its scale $\bar{r}_k^{(s)}$ by

$$\bar{r}_k^{(s)} = \min \left\{ r \in \mathbb{N}^* : 4r \Delta_k^2 \ge \kappa L^2 \left[ s_k \log \left( \frac{ep}{s_k} \right) + \log \left( \frac{n}{r\delta} \right) \right] \right\}. \quad (19)$$
$|C_{l,r,(1)}| \ge \dots \ge |C_{l,r,(p)}|$. Then, for $s \in \mathcal{Z}$, we define the partial CUSUM norm by

$$ \Psi_{l,r,s}^{(p)} = \sum_{i=1}^{s} (C_{l,r,(i)})^2 . \qquad (15) $$

Then, we define the test $T_{l,r}^{(p)}$ rejecting the null when at least one of the partial norms is large

$$ x_{r,s}^{(p)} := x_{r,s}^{(p)}(\delta) = 4s \log \left( \frac{2ep}{s} \right) + 4 \log \left( \frac{n}{r\delta} \right); \quad T_{l,r}^{(p)} = \max_{s \in \mathcal{Z}} \mathbf{1} \left\{ \Psi_{l,r,s}^{(p)} > x_{r,s}^{(p)} \right\}. $$

Finally, we define the sparse test by aggregating both the Berk-Jones test and the partial norm test. For any $(l,r) \in \mathcal{G}_D$, let $T_{l,r}^{(s)} = T_{l,r}^{(p)} \lor T_{l,r}^{(\text{BJ})}$. The next proposition controls the error of this collection of tests.

**Proposition 2.** There exists a universal constant $\kappa_s > 0$ and an event $\xi^{(s)}$ of probability larger than $1 - 4\delta$ such that (i) $T_{l,r}^{(s)} = 0$ for all $(l,r) \in \mathcal{H}_0 \cap \mathcal{G}_D$ and (ii) $T_{\bar{\tau}_k^{(s)}, \bar{\tau}_k^{(s)}}^{(s)} = 1$ for all $\kappa_s$-sparse high-energy change-point $\tau_k$.

Here we introduced two different statistics for the same sparse regime $s_k \le [p \log(n/(r_k \delta))]^{1/2}$ – the Berk-Jones statistic and the partial sums statistic – mainly to solve a problem of integrability. We made this choice for the sake of simplicity, but we could have used a single test, as presented in [28]

$$ \Psi_{x,l,r}^{(\text{LGS})} = \sum_{i=1}^{p} (C_{l,r,i}^2 - \mathbb{E}[Z|Z \ge x]) \mathbf{1}\{C_{l,r,i}^2 \ge x\}, $$

where $Z$ follows a standard normal distribution $N(0, 1)$. This statistic leads to the same type of result as the Berk-Jones statistic when enough coordinates $\mu_k - \mu_{k-1}$ are large in absolute value, and it is comparable to the partial sums statistic when its threshold $x$ becomes low enough.

### 3.3. Consequences

To conclude this section, it suffices to observe that, for $c_0$ in (8), any $c_0$-high-energy change-point $\tau_k$ in the sense of (8) is either a $c_0/2$-dense or a $c_0/2$-sparse high-energy change-point. Hence, upon defining the test $T_{l,r} = T_{l,r}^{(d)} \lor T_{l,r}^{(s)}$ for $(l,r) \in \mathcal{G}_D$, we consider the change-point procedure $\hat{\tau}$ defined in Algorithm 1. Gathering Theorem 1 with Proposition 1 and Proposition 2, we obtain the following.

**Corollary 2.** There exists a universal constant $c_0 > 0$ such that, with probability higher than $1 - 6\delta$, the estimator $\hat{\tau}$ satisfies (NoSp) and detects all $c_0$-high-energy change-points (as defined in (8)) $\tau_k$ in the sense

$$ d_{H,1}(\hat{\tau}, \tau_k) < \frac{r_k^*}{2} \le \frac{r_k}{2}, $$
For any $\kappa$-sparse high-energy change-point, it holds that $4(\bar{r}_k^{(s)} - 1) \le r_k$.

**Proposition 4.** There exists a numerical constant $\bar{c}_{\text{thresh}}^{(p)} > 0$ such that the following holds for any $\kappa_s > 32\bar{c}_{\text{thresh}}^{(p)}$. With probability higher than $1 - \delta$, one has (i) $T_{l,r}^{(p)} = 0$ for all $(l,r) \in \mathcal{G}_F \cap \mathcal{H}_0$ and (ii) $T_{\tau_k, \bar{r}_k^{(s)}}^{(p)} = 1$ for all $\kappa_s$-sparse high-energy change-point $\tau_k$ in the sense of (18).

As for Proposition 3, the proof relies on a careful analysis of the joint distributions of the statistics $\Psi_{l,r,s}^{(p)}$ to handle the multiplicity of $\mathcal{G}_F$.

## 4.3 Consequences

Let $c_0 > 0$ be some constant that we will discuss later. A change-point $\tau_k$ is then said to be a $c_0$-high-energy change-points –in the sub-Gaussian setting– if

$$r_k \Delta_k^2 \ge c_0 L^2 \left[ \left( \sqrt{p \log \left( \frac{n}{r_k \delta} \right)} \wedge \left( s_k \log \left( \frac{ep}{s_k} \right) \right) \right) + \log \left( \frac{n}{r_k \delta} \right) \right] . \quad (20)$$

We here re-introduce $\mathcal{K}^* \subset [K]$ as the subset of indices such that $\tau_k$ satisfies (20).

We gather both tests by considering, for any $(l,r) \in \mathcal{G}_F$, the test $T_{l,r} = T_{l,r}^{(d)} \vee T_{l,r}^{(p)}$ with tuning parameters $\bar{c}_{\text{thresh}}^{(d)}$ and $\bar{c}_{\text{thresh}}^{(p)}$ as in Propositions 3 and 4. Consider any $c_0 > 32(\bar{c}_{\text{thresh}}^{(d)} \vee \bar{c}_{\text{thresh}}^{(p)})$ and any $c_0$-high-energy change-point $\tau_k$, which is either a $c_0$-sparse or a $c_0$-dense high-energy change-point. Defining

$$\bar{r}_k = \bar{r}_k^{(d)} \wedge \bar{r}_k^{(s)}, \quad (21)$$

we straightforwardly derive from Proposition 3 and Proposition 4 the following result.

**Corollary 4.** There exists two numerical constants $\bar{c}_{\text{thresh}}^{(p)} > 0$ and $\bar{c}_{\text{thresh}}^{(d)} > 0$ such that the following holds. With probability higher than $1 - \delta$, it holds that (i) $T_{l,r} = 0$ for all $(l,r) \in \mathcal{G}_F \cap \mathcal{H}_0$ and (ii) $T_{\tau_k, \bar{r}_k} = 1$ for any $c_0$-high-energy change-point $\tau_k$ in the sense of (20).

Then, it suffices to combine this multiple testing procedure with Algorithm 1 to get the change-point procedure $\hat{\tau}$. Since, for a high-energy change-point in the sense of (20), we have $4(\bar{r}_k - 1) < r_k$, we are in position to apply Theorem 1.

**Corollary 5.** There exist two numerical constant $\bar{c}_{\text{thresh}}^{(p)} > 0$ and $\bar{c}_{\text{thresh}}^{(d)} > 0$ such that the following holds. With probability higher than $1 - \delta$, the estimator $\hat{\tau}$ satisfies (NoSp) and detects $c_0$-high-energy change-point $\tau_k$ (as defined in(20)), that is

$$d_{H,1}(\hat{\tau}, \tau_k) \le \bar{r}_k - 1 \le \frac{r_k}{4},$$

where $\bar{r}_k$ is defined in (21).

In the case where all change-points are $c_0$-high-energy change-points in the sense of (20), all of them are detected, and a result similar to Corollary 3 holds here, replacing $r_k^*/2$ by $\bar{r}_k - 1$. Also, both the Hausdorff distance and the Wasserstein distance, can be bounded as in Equation (16) if we replace $r_k^*/2$ by $\bar{r}_k - 1$.

As already stated, we could have obtained a similar result (but with different constants) using the dyadic grid $\mathcal{G}_D$ instead of $\mathcal{G}_F$. To conclude this section, let us compare the conditions (20) and (8) for high-energy. Define

$$\psi_{n,r,s}^{(sg)} = \sqrt{p\gamma_r} \wedge \left( s \log \left( \frac{ep}{s} \right) \right) + \gamma_r,$$
where $r_k^*$ is defined in (9).

If the change-points are of high-energy, that is $K^* = [K]$, then Corollary 2 can be reformulated as follows:

**Corollary 3.** Assume that for all $k = 1, \dots, K$, $\tau_k$ is a $c_0$-high-energy change-point (see (8)) where $c_0$ is the same as in Corollary 2. Then, with probability higher than $1 - 6\delta$, the estimator $\hat{\tau}$ satisfies $\hat{K} = K$ and

$$|\hat{\tau}_k - \tau_k| < \frac{r_k^*}{2} \le \frac{r_k}{2}, \quad \text{for all } k = 1, \dots, K.$$

In particular, one can respectively bound the Hausdorff and the Wasserstein losses, with probability higher than $1 - 6\delta$ by

$$d_H(\hat{\tau}, \tau) \le \max_{k=1,\dots,K} \frac{r_k^*}{2} \quad \text{and} \quad d_W(\hat{\tau}, \tau) \le \sum_{k=1,\dots,K} \frac{r_k^*}{2}. \qquad (16)$$

In Section 5, we establish that the Condition (8) is (up to a multiplicative constant) unimprovable and corresponds to the detection threshold for multivariate change-points.

Corollary 3 can be compared to the result of [46] on multivariate change-point detection in the multiple change-point setting. Using a method based on the CUSUM statistic and assuming that there are only high-energy change-points, the authors also obtain an upper bound on the energy necessary to detect the change-points. However, this result does not adapt to $r_k$, $\Delta_k$, $s_k$, and the detection rate is suboptimal in many regimes. Writing $r = \min_{k=1,\dots,K} r_k$, $\Delta = \min_{k=1,\dots,K} \Delta_k$ and $s = \max_{k=1,\dots,K} s_k$, Theorem 5 of [46] requires two conditions of the type $r\Delta^2 \ge c(\frac{n}{r})^4 \log(np)$ and $r\Delta^2 \ge c s \frac{n}{r} \log(np)$. This detection rate is therefore suboptimal by a polynomial factor in $n/r$ when $r$ is of smaller order than $n$, and by a logarithmic factor $\log(np)$ instead of $\log(1 + \sqrt{p/s}\log(n/r)) + \frac{1}{s}\log(n/r)$ when $r$ is of order $n$. Closer to our results, [5] have introduced another bottom-up procedure in the very specific asymptotic setting $n = e^{p\zeta}$ for $\zeta \in (0,1)$ with a fixed $K$ number of change-points. Assuming that, for each change-point, at least $s$ coordinates of $\mu_{k+1} - \mu_{k+1}$ are larger than $\zeta$ in absolute value, [5] establish that their procedure provably detects the change-points as long as

$$rs\zeta^2 \ge c \begin{cases} \sqrt{p \log(n)} & \text{if } s \ge 0.5 \sqrt{p \log(n)} \\ s \log\left(\frac{p}{s^2} \log(n)\right) & \text{if } s \le 0.5 \sqrt{p \log(n)} \end{cases}.$$

In their specific asymptotic regime and when all non-zero coordinates are of the same order, and all the change-points have a similar length $r_k$, their result is similar to ours up to the logarithmic terms. Indeed, for equispaced change-points, our logarithmic term $\log(n/r_k) = \log(K)$ is much smaller than $\log(n)$. Besides, their result does not handle the presence of low-energy change-points and does not hold beyond the asymptotic regime $n = e^{p\zeta}$. In contrast, our condition (8) for high-energy change-points entails that the detection conditions
where we recall that $\gamma_r = \log(\frac{n}{r\delta})$. If $\gamma_r \ge p/2$, then $\psi_{n,r,s}^{(sg)} \asymp \gamma_r$. In low dimension, the energy threshold for multivariate change-point detection is the same as in the univariate setting, see [38]. If $\gamma_r \le p/2$, then

$$ \psi_{n,r,s}^{(sg)} \asymp \begin{cases} \gamma_r & \text{if } s \le \frac{\gamma_r}{\log(p)-\log(\gamma_r)} \\ s \log\left(e^{\frac{p}{s}}\right) & \text{if } \frac{\gamma_r}{\log(p)-\log(\gamma_r)} < s < \frac{\sqrt{p\gamma_r}}{\log(p)-\log(\gamma_r)} \\ \sqrt{p\gamma_r} & \text{if } s \ge \frac{\sqrt{p\gamma_r}}{\log(p)-\log(\gamma_r)} \end{cases} $$

As a consequence, $\psi_{n,r,s}^{(sg)}$ and $\psi_{n,r,s}^{(g)}$ are of the same order of magnitude for all $s$ when $\gamma_r \ge p/2$. When $\log(n/r\delta) < p$, they are also of the same order of magnitude except when $s$ is close but smaller than $\sqrt{p\gamma_r}$, for which the ratio $\psi_{n,r,s}^{(sg)}/\psi_{n,r,s}^{(g)}$ between these two quantities can be as large as $\log(p) - \log(\gamma_r)$. This gap corresponds to the regime where the test based on the Berk-Jones statistic defined in Equation (14), used in the Gaussian case, outperforms the test based on the partial CUSUM norm statistic defined in Equation (15).

In the definitions of the tests, the tuning constants $\bar{c}_{\text{thresh}}^{(p)}$ and $\bar{c}_{\text{thresh}}^{(d)}$ are left implicit, although one can find suitable values by following the proofs of Propositions 3 and 4. In practice, the practitioner can calibrate them by a Monte-Carlo method by simulating a Gaussian multivariate times series without any change-points. Then, $\bar{c}_{\text{thresh}}^{(p)}$ and $\bar{c}_{\text{thresh}}^{(d)}$ are chosen so that the Family-wise error rate (FWER) of the two collections $(T_{l,r}^{(d)})$ and $T_{l,r}^{(p)}$ is equal to $\delta$.

**Computational Cost.** The computational cost of the statistic $T_{l,r} = T_{l,r}^{(d)} \lor T_{l,r}^{(p)}$ is $O(p(r + \log(p)))$. Thus, a naive computation of all the tests $T_{l,r}$ for $(l,r)$ in the complete grid $\mathcal{G}_F$ requires $O(p \log(p) \sum_{(l,r) \in \mathcal{G}_F} r) = O(np(n^2 + \log(p)))$ operations. Nevertheless, using the fact that $\sum_{i=l+1}^{l+r} Y_i = (\sum_{i=l}^{l+r-1} Y_i) + Y_{l+r} - Y_l$, it is possible to compute all the tests at scale $r$ with cost $O(np \log(p))$. Since there are $n$ possible scales $r$ on the complete grid, the whole procedure cost is $O(n^2p \log(p))$. Using a grid $\mathcal{G} = \{(l,r) \in \mathcal{G}_F : r \in \mathcal{R}\}$ that contains dyadic scales and all possible locations $l$ for each scale, the whole change-point detection would then require only $O(np \log(n) \log(p))$ computations, since there are only $\log(n)$ possible scales $r$ for such grids.

# 5 Minimax lower bound

In this section, we write for any $\Theta \in \mathbb{R}^{p \times n}$, the distribution of the time series $Y = (y_1, \dots, y_n)$ in the model (2) with Gaussian noise $\varepsilon_t \sim \mathcal{N}(0, \sigma^2 \mathbf{I}_p)$. In Section 3, we have established that any change-point satisfying the condition (8), that is

$$ r_k \Delta_k^2 \ge c_0 \sigma^2 \left[ s_k \log \left( 1 + \frac{\sqrt{p}}{s_k} \sqrt{\log \left( \frac{n}{r_k \delta} \right)} \right) + \log \left( \frac{n}{r_k \delta} \right) \right], $$

is detected by our change-point procedure. We now show that this energy condition is unimprovable from a minimax point of view. More precisely, let us define, for any $u > 0$, the class $\tilde{P}(u)$ of mean parameters $\Theta$ with arbitrary $K \ge 0$ number of change points and such that any change-point $\tau_k$ for $1 \le k \le K$ satisfies

$$ r_k \Delta_k^2 \ge \frac{1}{2} \sigma^2 \left[ s_k \log \left( 1 + u \frac{\sqrt{p}}{s_k} \sqrt{\log \left( \frac{n}{r_k} \right)} \right) + u \log \left( \frac{n}{r_k} \right) \right]. \quad (22) $$

For $u$ small enough, it turns out no change-point estimator is able to detect all change-points without estimating any spurious change-point with high probability on the full class $\tilde{P}(u)$. Still,
are qualitatively different for other scalings in $n$ and $p$. On the technical side, our condition (8) is of $l_2$ type whereas that in [5] is of minimal non-zero type. Recovering the tight $l_2$ conditions turns out to be much more challenging as we need to handle situations where some coordinates have different orders of magnitude. This is the main reason why we need to resort to a combination of the Berk-Jones and the partial-norm statistics.

**Comparison to one change-point problem.** When one knows that $K \le 1$ (at most one change-point), then [28] proved that it is possible to detect $\tau_1$ if and only if $r_1\Delta_1^2 \ge c\sigma^2[s_1 \log(1 + \frac{1}{s_1}\sqrt{p \log \log 8n}) + \log \log 8n]$. As in the univariate setting, the problem with only one change-point is simpler than for general $K \ge 2$. As for our procedure, Liu et al. [28] rely on statistics based on the CUSUM – a chi square statistics in the dense case and a thresholded sum of squared coordinates in the sparse case – to detect and localize $\tau_1$. It turns out that the detection procedure of [28] adapts to distance $r_1 = \max(\tau_1-1, n+1-\tau_1)$ the boundary, and one could refine their result by stating that $\tau_1$ is detectable if and only if $r_1\Delta_1^2 \ge c\sigma^2[s_1 \log(1 + \frac{1}{s_1}\sqrt{p \log \log(2n/r_1)}) + \log \log(2n/r_1)]$ which is more smaller when $r_1$ is of the order of $n$. This refined result is in the same spirit as our bounds for multiple change-point, but the rate is faster because one obtains $\log \log(n/r_1)$ – instead of $\log(n/r_k)$ in our case. The reason for this faster rate is due to the relative simplicity of the problem with only one change-point. Indeed, in single change-point detection, there is no need to look for change-points at all positions and scale at the same time, since scale and positions are related. This implies that it is possible to attain faster rates than in multiple change-point detection. The comparison between single and multiple change-point detection is thoroughly done in [38] for univariate models.

**Computational Cost.** The cost of the tests $T_{l,r}^{(d)}$ in the dense regime is $O(rp)$. The computation of the partial norm statistic requires to sort the coordinates $C_{l,r,i}$ of the CUSUM statistic, which takes $O(p(r + \log(p)))$ operations. Since only the thresholds $x \le c \log(np/(r\delta))^{1/2}$ are needed to compute the Berk-Jones statistic, it holds that, for $\delta \ge (np)^{-c}$ with a numerical constant $c > 0$, the computational cost of the Berk-Jones statistic is $O(p(r + \log(np)))$. Thus, for each $(l,r)$, the overall computational cost of the test $T_{l,r} = T_{l,r}^{(d)} \cup T_{l,r}^{(s)}$ is $\Lambda = O(p(r + \log(np)))$, and the computational cost of the whole change-point detection procedure on the dyadic grid is $O(np \log(np))$.

## 4. Multi-scale change-point detection with sub-Gaussian noise

We now turn to the more general case of sub-Gaussian distributions [37]. Given a random variable $Z$, define its $\psi_2$-norm by $\|Z\|_{\psi_2} = \inf\{x > 0, \mathbb{E}[\exp(Z^2/x^2)] \le 2\}$. Given $L > 0$, a mean zero real random variable is said to be $L$-sub-Gaussian if $\|Z\|_{\psi_2} \le L$. This implies in particular that, for all $x \ge 0$, one has $\mathbb{P}(|Z| \ge x) \le 2\exp(-x^2/L^2)$. Throughout this section, we assume that, for $t = 1, \dots, n$, the random vectors $\epsilon_t$ are independent, have independent
using this large class provides somewhat pessimistic bounds. For instance, the most challenging distributions in $\bar{P}(u)$ for the purpose of change-point detection satisfy $s_k = p$ and $r_k = 1$ (very close change-points). As a consequence, relying on the full collection $\bar{P}(u)$ turns too pessimistic. To establish that our bounds are adaptive with respect to the sparsity $s_k$ and the length $r_k$, we define, for any positive integers $1 \le r \le \lfloor n/2 \rfloor$ and any $1 \le s \le p$ the collection

$$ \bar{P}(u, r, s) = \{ \Theta \in \bar{P}(u) : \min_k r_k \ge r \text{ and } \max_k s_k \le s \}. $$

By convention, constant means $\Theta$ with no change-points ($K = 0$) also belong to $\bar{P}(u, r, s)$. In the class $\bar{P}(u, r, s)$, all change-points have a sparsity at most $s$ and a length at least $r$. Hence, $\bar{P}(u, r, s)$ becomes larger when $s$ increases or when $r$ increases.

**Theorem 2.** Fix any $u \in (0, 1/8)$. For any $\sigma > 0$, $n \ge 2$, $p \ge 1$, any length $1 \le r \le n/4$, and any sparsity $1 \le s \le p$, we have

$$ \inf_{\hat{\tau}} \sup_{\Theta \in \bar{P}(u,r,s)} \mathbb{P}_{\Theta}(\hat{K} \neq K) \ge \frac{1}{4}, $$

where the infimum is taken over all estimators $\hat{\tau}$ of the change-point vector $\tau$ and $\hat{K} = |\hat{\tau}|$.

Thus, in the Gaussian setting, if all the change-points have a high-energy in the sense of (8) but with a smaller multiplicative constant factor, no change-point estimator can consistently estimate the true number of change-points. The next corollary restates this negative results in the same lines as Corollary 3.

**Corollary 6.** Fix any $u \in (0, 1/8)$. For any $\sigma > 0$, $n \ge 2$, $p > 1$, any length $1 \le r \le n/4$, any sparsity $1 \le s \le p$, and any estimator $\hat{\tau}$, there exists some $\Theta \in \bar{P}(u, r, s)$ such that with $\mathbb{P}_{\Theta}$-probability larger than $1/4$, at least one of the two following properties is satisfied

* $\hat{\tau}$ contains at least one spurious change-point

* at least a change-point $\tau_k$ with $1 \le k \le K$ is not detected, i.e. there is no change-point estimated in the interval $[(\tau_{k-1} + \tau_k)/2, (\tau_k + \tau_{k+1})/2]$.

This corollary is to be compared to Corollary 3 - indeed, the energy condition in Equation (22) differs from Equation (8) only by a numerical multiplicative constant. As a consequence, the energy condition (22) is minimal for detection by a change-point estimator that achieves (NoSp).

# 6 Application to other change-point problems

In this section, we apply the general methodology of Section 2 to two other problems, namely detection of covariance and nonparametric change-points. This allows us to obtain the first tight minimax detection conditions for these problems.

## 6.1 Covariance change-point detection

Following Wang et al. [40], we consider the covariance change-point model where the covariance matrices $\Sigma_t$ of the centered random vectors $y_t \in \mathbb{R}^p$ are piece-wise constant. Then, the goal is to estimate the times $0 < \tau_1 < \dots < \tau_K < \tau_{K+1} = n+1$ such that $\Sigma_t$ is varying. See [40] for motivations. As in that work, we assume that the random vectors $y_t$ are independent and are sub-Gaussian with a uniformly bounded Orlicz norm, that is $\max_{t=1,\dots,n} \|y_t\|_{\psi_2} \le B$ for some known fixed $B$. The Orlicz norm of a random vector $y$ is the supremum of the Orlicz norm of any
L-sub-Gaussian components $\varepsilon_{t,i}$, for $i = 1, \dots, p$ with variance $\sigma^2$. As in the previous section, we apply the general aggregation procedures introduced in Section 2. As a consequence, our main task boils down to introducing a near-optimal multiple testing procedure indexed by a grid for detecting the existence of a change-point. Here, we shall rely on the complete grid $\mathcal{G}_F = J_n = \{(l,r) : r = 1, \dots, \lfloor \frac{n}{2} \rfloor \text{ and } l = r + 1, \dots, n - r\}$ whose size is quadratic with respect to $n$. All the results presented in this section are still valid (but with different numerical constants) if we keep the dyadic grid $\mathcal{G}_D$ as in the previous section. Here, we use the complete grid as a proof of concept that one can rely on the full collection of possible segments without deteriorating the rates. Still, controlling the behavior of the procedure on the complete grid is technically more involved and requires chaining arguments. A detailed comparison between the complete and dyadic grids is made in Section 7.

In order to emphasize the common points with the previous section, we use the same notation $\mathcal{K}^*$ for the collection of high-energy change-points<fnref n="1" />, $\bar{r}_k$ for the scales associated to the k-th change-points<fnref n="2" />, $\Psi$ for the statistics, $T$ for the test and $x$ for the thresholds although these quantities are slightly changed to cope with the sub-Gaussian tail distribution. We follow the same scheme as for the Gaussian case and first introduce multi-scale tests for dense change-points before turning to sparse change-points. As in the previous section, we consider some $\delta \in (0, 1)$ corresponding to the type I error probability.

## 4.1. Dense change-points with sub-Gaussian noise

Recall that, for a change-point $\tau_k$, $s_k$ stands for the sparsity of the difference $\mu_{k+1} - \mu_k$. We focus here on dense change-points for which $s_k$ is possibly as large as $p$. Given $\kappa > 0$, $\tau_k$ is a $\kappa$-dense high-energy change-point if

$$r_k \Delta_k^2 \geq \kappa L^2 \left( \sqrt{p \log \left( \frac{n}{r_k \delta} \right)} + \log \left( \frac{n}{r_k \delta} \right) \right). \quad (17)$$

This condition is very similar to its counterpart (10) for Gaussian noise. Still, we introduce it here for the sake of completeness. For $k \in [K]$ such that $\tau_k$ is a $\kappa$-dense high-energy change-point, we define $\bar{r}_k^{(d)}$ as the minimum length such that an inequality similar to (17) is satisfied:

$$\bar{r}_k^{(d)} = \min \left\{ r \in \mathbb{N}^* : 4r \Delta_k^2 \geq \kappa L^2 \left( \sqrt{p \log \left( \frac{n}{r\delta} \right)} + \log \left( \frac{n}{r\delta} \right) \right) \right\}.$$

As in the Gaussian case in Section 3, $\bar{r}_k^{(d)}$ corresponds to the smallest scale such that $\tau_k$ is guaranteed to be detected. For any $\kappa$-dense high-energy change-point, it holds that $4(\bar{r}_k^{(d)} - 1) < r_k$. For any positive integers $(l, r) \in \mathcal{G}_F$, we consider the same CUSUM-based statistic $\Psi_{l,r}^{(d)} := \|C_{l,r}\|^2 - p$ as for Gaussian noise.

<footnote n="1">¹See Equation (20) as the energy condition is slightly different in the sub-Gaussian setting.</footnote>

<footnote n="2">²Re-defined in Equation (21).</footnote>
uni-dimensional projection of $y$ - see e.g. [37]. If the $y_t$'s follow a normal distribution, this amounts to assuming that $\max_{t=1,...,n} \|\Sigma_t\|_{op} \le 2B^2$ where $\|\cdot\|_{op}$ is for the operator norm. The purpose of Wang et al. was to detect small changes in operator norm, that is detecting instants $\tau_k$ such that $\Sigma_{\tau_k} \ne \Sigma_{\tau_{k-1}}$ with $\|\Sigma_{\tau_k} - \Sigma_{\tau_{k-1}}\|_{op}$ possibly small. Apart from the operator norm, other norms have also been considered e.g. in [10]. Here, we focus on the operator norm as in [40].

Recalling the generic procedure introduced in Section 2, we consider the dyadic grid $\mathcal{G}_D$ and some $\delta \in (0, 1)$. For any $(l, r) \in \mathcal{G}$, we respectively write $\hat{\Sigma}_{l,-r}$ and $\hat{\Sigma}_{l,r}$ for the empirical covariance matrices

$$ \hat{\Sigma}_{l,-r} = r^{-1} \sum_{t=l-r}^{l-1} y_t y_t^T ; \quad \hat{\Sigma}_{l,r} = r^{-1} \sum_{t=l}^{l+r-1} y_t y_t^T . $$

Then, we consider the test $T_{l,r}$ rejecting for large values of $\|\hat{\Sigma}_{l,r} - \hat{\Sigma}_{l,-r}\|_{op}$.

$$ T_{l,r} = \mathbf{1} \left\{ \|\hat{\Sigma}_{l,r} - \hat{\Sigma}_{l,-r}\|_{op} \ge c_0 B^2 \left[ \sqrt{\frac{p}{r}} + \frac{p}{r} + \sqrt{\frac{\log(\frac{2n}{\delta r})}{r}} + \frac{\log(\frac{2n}{\delta r})}{r} \right] \right\}, \quad (23) $$

where the numerical tuning constant $c_0$ is set in the proof of the following proposition. Relying on concentration bounds [23] for the empirical covariance matrix of sub-Gaussian random vectors, we easily prove that the FWER of the multiple testing procedure $(T_{l,r})$ with $(l, r) \in \mathcal{G}_D$ is small. Then, we can analyze the type II error probability and plug it into the generic result (Theorem 1) to control the behavior of the change-point estimator $\hat{\tau}$. This leads us to the following result. In the sequel, a change-point $\tau_k$ is said to have a high-energy if

$$ r_k \|\Sigma_{\tau_k} - \Sigma_{\tau_{k-1}}\|_{op}^2 \ge c_1 B^4 \left[ \left( p + \log \left( \frac{2n}{r_k \delta} \right) \right) \wedge r_k \right], \quad (24) $$

where the numerical constant $c_1$ is introduced in the proof of the following proposition. We recall that, by definition of the model, we have $\|\Sigma_{\tau_k} - \Sigma_{\tau_{k-1}}\|_{op} \le 4B^2$.

**Proposition 5.** There exist positive numerical constants $c_0$, $c_1$, and $c_2$ such that the following holds for any $B > 0$ and any sequence of independent centered random vectors $(y_t)$ satisfying $\max_t \|y_t\|_{\psi_2} \le B$. With probability higher than $1-\delta$, the change-point estimator $\hat{\tau}$ satisfies (NoSp) and detects all high-energy change-points in the sense of (24). Besides, any such high-energy change-point $\tau_k$ satisfies

$$ d_{H,1}(\hat{\tau}, \tau_k) \le c_2 B^4 \frac{p + \log(2\delta^{-1}B^{-4}n \|\Sigma_{\tau_k} - \Sigma_{\tau_{k-1}}\|_{op}^2)}{\|\Sigma_{\tau_k} - \Sigma_{\tau_{k-1}}\|_{op}^2} \le \frac{r_k}{4}, \quad (25) $$

under the same event of probability than $1 - \delta$.

Let us compare our condition (24) for detection with Theorem 2 in Wang et al. [40]. The authors assume that all the change-points satisfy

$$ \min_k r_k \min_k \|\Sigma_{\tau_k} - \Sigma_{\tau_{k-1}}\|_{op}^2 \ge c_1' B^4 p \log(n) . $$

In addition to the fact that we allow some change-points to have an arbitrarily low energy, our requirement for detection scales like $\sqrt{p} + \sqrt{\log(n/r_k)}$ instead of $\sqrt{p} \log(n)$.

The next proposition establishes that the latter condition is minimal. By homogeneity, we can only consider the case where $B = 3/2$. We focus our attention on Gaussian distributions so
Let $\bar{c}_{\text{thresh}}^{(d)} > 0$ be a tuning parameter to be discussed later. To calibrate the corresponding multiple test procedures ($T_{l,r}^{(d)}$) with $(l, r) \in \mathcal{G}_F$ rejecting for large values of $\Psi_{l,r}^{(d)}$ we introduce

$$T_{l,r}^{(d)} := \mathbf{1}\{\Psi_{l,r}^{(d)} > x_r^{(d)}\}; \quad x_r^{(d)} = \bar{c}_{\text{thresh}}^{(d)} \frac{L^2}{\sigma^2} \left(\sqrt{p \log\left(\frac{n}{r\delta}\right)} + \log\left(\frac{n}{r\delta}\right)\right).$$

**Proposition 3.** There exists a numerical constant $\bar{c}_{\text{thresh}}^{(d)} > 0$ such that the following holds for any $\kappa_d > 32\bar{c}_{\text{thresh}}^{(d)}$. With probability higher than $1 - \delta$, one has (i) $T_{l,r}^{(d)} = 0$ for all $(l, r) \in \mathcal{G}_F \cap \mathcal{H}_0$ and (ii) $T_{\tau_k, \bar{r}_k^{(d)}}^{(d)} = 1$ for all $\kappa_d$-dense high-energy change-points $\tau_k$.

In comparison to Proposition 1 in the previous section, there are two differences. First, we need to cope with sub-Gaussian distribution by applying the Hanson-Wright inequality. Most importantly, the grid $\mathcal{G}_F$ is much larger than $\mathcal{G}_D$ so that we cannot simply consider each test $T_{l,r}$ separately and simply apply a union bound as in the previous section. To handle the dependencies between the statistics $\Psi_{l,r}^{(d)}$, we have to apply a chaining argument. In fact, the thresholds $x_r^{(d)}$ are similar to their counterpart in the previous section, whereas the number $|\mathcal{G}_F|$ of tests is now proportional to $n^2$. In principle, the benefit of using the full grid $\mathcal{G}_F$ is that $(\tau_k, \bar{r}_k^{(d)})$ belongs to $\mathcal{G}_F$ so that we can consider the CUSUM statistic based on a segment $[\tau_k - \bar{r}_k^{(d)}, \tau_k + \bar{r}_k^{(d)}]$ centered around the change-point $\tau_k$. In contrast, $(\tau_k, \bar{r}_k^{(d)})$ does not necessarily belong to the dyadic grid $\mathcal{G}_D$ and we needed to consider its best approximation $(\bar{\tau}_k^{(d)}, \bar{r}_k^{(d)})$. The segment $[\bar{\tau}_k^{(d)} - \bar{r}_k^{(d)}, \bar{\tau}_k^{(d)} + \bar{r}_k^{(d)}]$ is therefore not centered on $\tau_k$ and the corresponding statistic $\Psi_{\bar{\tau}_k^{(d)}, \bar{r}_k^{(d)}}^{(d)}$ is in expectation smaller than $\Psi_{\tau_k, \bar{r}_k^{(d)}}^{(d)}$. In summary, both the collections of dense tests $\Psi_{l,r}^{(d)}$ on $\mathcal{G}_D$ and $\mathcal{G}_F$ are able to detect change-points whose energy is, up to some multiplicative constants, higher than $L^2\left[[p \log(\frac{n}{r_k \delta})]^{1/2} + \log(\frac{n}{r_k \delta})\right]$.

## 4.2. Sparse change-points with sub-Gaussian noise

Unlike in the Gaussian case, we do not know the exact distribution of the noise. As a consequence, the Berk-Jones test and more generally higher-criticism type tests cannot be applied to this setting. This is why we only rely on the partial norm statistic. Recall that $\mathcal{Z} = \{1, 2, 2^2, \dots, 2^{\lfloor \log_2(p) \rfloor}\}$ stands for a dyadic set of sparsities. For $(l, r) \in \mathcal{G}_F$ and $s \in \mathcal{Z}$, we also recall that the partial CUSUM norm is defined as $\Psi_{l,r,s}^{(p)} = \sum_{i=1}^s (C_{l,r,(i)})^2$. Then, for any $(l, r) \in \mathcal{G}_F$, the test $T_{l,r}^{(p)}$ rejects the null when at least one of the partial norms is large

$$x_{r,s}^{(p)} = s + \bar{c}_{\text{thresh}}^{(p)} \frac{L^2}{\sigma^2} \left[ s \log \left( \frac{2ep}{s} \right) + \log \left( \frac{n}{r\delta} \right) \right]; T_{l,r}^{(p)} = \max_{s \in \mathcal{Z}} \mathbf{1}\{\Psi_{l,r,s}^{(p)} > x_{r,s}^{(p)}\},$$
that the distribution of the sequence $(y_1, \dots, y_n)$ is uniquely defined by the sequence $(\Sigma_1, \dots, \Sigma_n)$ of covariance matrices. Given an integer $1 \le r \le n/4$ and $\zeta \in (0, 1/\sqrt{2})$, we define $\tilde{P}(r, \zeta)$ the collection of sequences $\eta = (\Sigma_1, \dots, \Sigma_n)$ of covariance matrices that satisfy either $\Sigma_t = I_p$ or $\|\Sigma_t\|_{op} = 1 + \zeta$. Besides, the corresponding change-points $(\tau_1, \dots, \tau_K)$ of $\eta$ must satisfy $\min_k r_k \ge r$ and $\min_k \|\Sigma_{\tau_k} - \Sigma_{\tau_{k-1}}\|_{op} \ge \zeta$. For $\eta \in \tilde{P}(r, \zeta)$, we write $\mathbb{P}_\eta$ for the corresponding distribution of $(y_1, \dots, y_n)$.

**Proposition 6.** There exists a positive numerical constant $c$ such that, for any $n, p$ and any length $1 \le r \le n/4$ the following holds. Provided that $r\zeta^2 \le c(p + \log(n/r)) \wedge \frac{r}{2}$, we have

$$ \inf_{\hat{\tau}} \sup_{\eta \in \tilde{P}(r, \zeta)} \mathbb{P}_{\eta}(\hat{K} \neq K) \ge \frac{1}{4}. $$

As a consequence, our procedure $\hat{\tau}$ achieves the minimal separation condition (24) for change-point detection. In their work, [40] obtain faster localization errors than (25) to the price of stronger separation conditions. Our focus in this work is to provide optimal detection conditions and we did not try to optimize (24).

## 6.2 Univariate nonparametric change-point detection

We now turn to the univariate nonparametric change-point model considered in [33]. Let $m \ge 1$ be any positive integer. At each time $t = 1, \dots, n$, the random vector $y_t$ is an $m$-sample of a univariate distribution with cumulative distribution function $F_t$. Then, we aim at detecting a vector $\tau = (\tau_1, \dots, \tau_K)$ of change-points such that $F_{\tau_k} \neq F_{\tau_{k-1}}$. As in [33], we quantify the distance between two distributions by the Kolmogorov distance $\|F_1 - F_2\|_{\infty} = \sup_{z \in \mathbb{R}} |F_1(z) - F_2(z)|$.

As in the previous subsection, we build a procedure $\hat{\tau}$ with our generic algorithm on the dyadic grid. Regarding the collection of tests $(T_{l,r})$, we consider two-sample Kolmogorov-Smirnov tests. More precisely, we denote $\hat{F}_t$ the empirical distribution function associated with the sample $y_t$ and we define the test

$$ T_{l,r} = 1 \left\{ \left\| r^{-1} \left( \sum_{t=l}^{l+r-1} \hat{F}_t - \sum_{t=l-r}^{l-1} \hat{F}_t \right) \right\|_{\infty} \ge \sqrt{2 \frac{\log(4n/(\delta r))}{mr}} \right\}. $$

In the following, a change-point $\tau_k$ is said to have a high-energy if

$$ r_k \|F_{\tau_k} - F_{\tau_{k-1}}\|_{\infty}^2 \ge \frac{c_1}{m} \log \left( \frac{n}{r_k \delta} \right), \qquad (26) $$

where the numerical constant $c_1$ is introduced in the proof of the next proposition. As in Subsection 6.1, it is straightforward to prove, based on Dvoretzky-Kiefer-Wolfowitz inequality, that the FWER of the multiple testing procedures $(T_{l,r})$ with $(l, r) \in \mathcal{G}_D$ is small. Then, we analyze the type II error probability of this test and plug it into the generic result (Theorem 1) to control the behavior of the change-point estimator $\hat{\tau}$.

**Proposition 7.** There exist positive numerical constants $c_1$ and $c_2$ such that the following holds. With probability higher than $1 - \delta$, the change-point estimator $\hat{\tau}$ satisfies (NoSp) and detects all high-energy change-points $\tau_k$ in the sense of (26). Besides, any such high-energy change-points $\tau_k$ satisfies

$$ d_{H,1}(\hat{\tau}_{k'}, \tau_k) \le c_2 \frac{\log(\delta^{-1}nm \|F_{\tau_k} - F_{\tau_{k-1}}\|_{\infty}^2)}{m \|F_{\tau_k} - F_{\tau_{k-1}}\|_{\infty}^2} \le \frac{r_k}{4}, \qquad (27) $$

under the same event of probability than $1 - \delta$.
where $\bar{c}_{\text{thresh}}^{(p)}$ is a tuning parameter in Proposition 4 below. The partial norm test alone is not able to detect sparse high-energy change-points in the sense of (11) and we need to introduce a stronger condition on the energy. Given $\kappa > 0$, a change-point $\tau_k$ is a $\kappa$-sparse high-energy change-point in the sub-Gaussian setting if $s_k \le [p \log(\frac{n}{r_k \delta})]^{1/2}$ and

$$r_k \Delta_k^2 \ge \kappa L^2 \left[ s_k \log \left( \frac{ep}{s_k} \right) + \log \left( \frac{n}{r_k \delta} \right) \right] . \quad (18)$$

Both Conditions (11) and (18) are compared at the end of the subsection. For a $\kappa$-sparse high-energy change-point $\tau_k$, we define its scale $\bar{r}_k^{(s)}$ by

$$\bar{r}_k^{(s)} = \min \left\{ r \in \mathbb{N}^* : 4r \Delta_k^2 \ge \kappa L^2 \left[ s_k \log \left( \frac{ep}{s_k} \right) + \log \left( \frac{n}{r_k \delta} \right) \right] \right\} . \quad (19)$$

For any $\kappa$-sparse high-energy change-point, it holds that $4(\bar{r}_k^{(s)} - 1) \le r_k$.

**Proposition 4.** There exists a numerical constant $\bar{c}_{\text{thresh}}^{(p)} > 0$ such that the following holds for any $\kappa_s > 32\bar{c}_{\text{thresh}}^{(p)}$. With probability higher than $1 - \delta$, one has (i) $T_{l,r}^{(p)} = 0$ for all $(l,r) \in \mathcal{G}_F \cap \mathcal{H}_0$ and (ii) $T_{\tau_k, \bar{r}_k^{(s)}}^{(p)} = 1$ for all $\kappa_s$-sparse high-energy change-point $\tau_k$ in the sense of (18).

As for Proposition 3, the proof relies on a careful analysis of the joint distributions of the statistics $\Psi_{l,r,s}^{(p)}$ to handle the multiplicity of $\mathcal{G}_F$.

### 4.3. Consequences

Let $c_0 > 0$ be some constant that we will discuss later. A change-point $\tau_k$ is then said to be a $c_0$-high-energy change-points –in the sub-Gaussian setting– if

$$r_k \Delta_k^2 \ge c_0 L^2 \left[ \left( \sqrt{p \log \left( \frac{n}{r_k \delta} \right)} \wedge \left( s_k \log \left( \frac{ep}{s_k} \right) \right) \right) + \log \left( \frac{n}{r_k \delta} \right) \right] . \quad (20)$$

We here re-introduce $\mathcal{K}^* \subset [K]$ as the subset of indices such that $\tau_k$ satisfies (20).

We gather both tests by considering, for any $(l,r) \in \mathcal{G}_F$, the test $T_{l,r} = T_{l,r}^{(d)} \vee T_{l,r}^{(p)}$ with tuning parameters $\bar{c}_{\text{thresh}}^{(d)}$ and $\bar{c}_{\text{thresh}}^{(p)}$ as in Propositions 3 and 4. Consider any $c_0 > 32(\bar{c}_{\text{thresh}}^{(d)} \vee \bar{c}_{\text{thresh}}^{(p)})$ and any $c_0$-high-energy change-point $\tau_k$, which is either a $c_0$-sparse or a $c_0$-dense high-energy change-point. Defining

$$\bar{r}_k = \bar{r}_k^{(d)} \wedge \bar{r}_k^{(s)}, \quad (21)$$

we straightforwardly derive from Proposition 3 and Proposition 4 the following result.

**Corollary 4.** There exists two numerical constants $\bar{c}_{\text{thresh}}^{(p)} > 0$ and $\bar{c}_{\text{thresh}}^{(d)} > 0$ such that the following holds. With probability higher than $1 - \delta$, it holds that (i) $T_{l,r} = 0$ for all $(l,r) \in \mathcal{G}_F \cap \mathcal{H}_0$ and (ii) $T_{\tau_k, \bar{r}_k} = 1$ for any $c_0$-high-energy change-point $\tau_k$ in the sense of (20).
In [33], the authors introduce a procedure detecting all the change-points provided that

$$ \min_k r_k \min_k \|F_{\tau_k} - F_{\tau_{k-1}}\|_{\infty}^2 \geq c_1 \frac{\log(n)}{m} . $$

Comparing this last condition with (26), we observe that our logarithmic term is tighter and that we allow arbitrarily low-energy change-points.

The next proposition establishes that the condition (26) is unimprovable. Given an integer $1 \le r \le n/4$ and $\zeta \in (0, 1/4)$, we focus our attention on the collection $\tilde{P}(r, \zeta)$ of sequences $(F_1, \dots, F_n)$ of distributions such that the corresponding change-points $(\tau_1, \dots, \tau_K)$ satisfy $\min_k r_k \ge r$ and $\min_k \|F_{\tau_k} - F_{\tau_{k-1}}\|_{\infty} \ge \zeta$. For $\eta \in \tilde{P}(r, \zeta)$, we write $\mathbb{P}_{\eta}$ for the corresponding distribution of the sequence $(y_1, \dots, y_n)$.

**Proposition 8.** There exists a positive numerical constant $c$ such that, for any $n, p$ and any length $1 \le r \le n/4$ the following holds. Provided that $r\zeta^2 \le c'\log(n/r)/m$, we have

$$ \inf_{\hat{\tau}} \sup_{\eta \in \tilde{P}(r, \zeta)} \mathbb{P}_{\eta}(\hat{K} \neq K) \geq \frac{1}{4}. $$

# 7 Discussion

## 7.1 Noise distribution for multivariate change-point detection

**Comparison between Gaussian and sub-Gaussian rates.** In this work, we have studied two types of noise distribution: Gaussian (Section 3) and general sub-Gaussian distributions (Section 4) without further knowledge on the distribution functions. Since the Gaussian setting is a specific instance of the sub-Gaussian setting, it is clear that the minimax lower bounds from Section 5 apply in both settings. As described in the previous subsection, the performances in the sub-Gaussian case almost match those in the Gaussian setting except for $s_k$ slightly lower but close to $\sqrt{p \log(en/r_k)}$. Indeed, in that regime, Berk-Jones or Higher-Criticism type statistics heavily rely on the probability distribution function of the noise, which is not available in the general sub-Gaussian case. Still, we could slightly improve the sub-Gaussian rates if we further assume that the noise components are identically distributed with common CDF $F$.

* If $F$ is known (know noise distribution), then one may adapt Berk-Jones test by replacing $\bar{\Phi}(x)$ in Equation (14) by $F(-x) + (1 - F(x))$. This would allow us to recover the exact same detection condition as in the Gaussian setting.

* If $F$ is unknown and if there are not too many change-points, one could hope to estimate the quantiles of the CUSUM statistic at each scale $r$ and plug them into a Berk-Jones statistics. This goes however beyond the scope of this paper.

**Unknown variance or more general variance matrix.** We assumed in the sparse multivariate sections that the variance $\sigma^2$ is known. Whereas the partial norm test only requires the knowledge of an upper bound on $\sigma$, the dense statistic $\Psi_{l,r}^{(d)}$ requires the exact knowledge of the variance. As soon as there are not too many change-points, it is possible to roughly estimate $\sigma$ and therefore accommodate the partial norm test with an unknown variance. In contrast, the dense statistics needs to be replaced by a $U$-statistics. Consider any even positive integer $r$ and define

$$ \tilde{\mathbf{C}}_{l,r}(Y) = \frac{\sqrt{r}}{2} \left( \frac{2}{r} \sum_{t=1}^{r/2} Y_{l-2(t-1)-1} - \frac{2}{r} \sum_{t=1}^{r/2} Y_{l+2(t-1)} \right), \quad \tilde{\mathbf{C}}'_{l,r}(Y) = \frac{\sqrt{r}}{2} \left( \frac{2}{r} \sum_{t=1}^{r/2} Y_{l-2t} - \frac{2}{r} \sum_{t=1}^{r/2} Y_{l+2(t-1)+1} \right), $$
Then, it suffices to combine this multiple testing procedure with Algorithm 1 to get the change-point procedure $\hat{\tau}$. Since, for a high-energy change-point in the sense of (20), we have $4(\bar{r}_k - 1) < r_k$, we are in position to apply Theorem 1.

**Corollary 5.** There exist two numerical constant $\bar{c}_{\text{thresh}}^{(p)} > 0$ and $\bar{c}_{\text{thresh}}^{(d)} > 0$ such that the following holds. With probability higher than $1 - \delta$, the estimator $\hat{\tau}$ satisfies (NoSp) and detects $c_0$-high-energy change-point $\tau_k$ (as defined in (20)), that is

$$d_{H,1}(\hat{\tau}, \tau_k) \le \bar{r}_k - 1 \le \frac{r_k}{4},$$

where $\bar{r}_k$ is defined in (21).

In the case where all change-points are $c_0$-high-energy change-points in the sense of (20), all of them are detected, and a result similar to Corollary 3 holds here, replacing $r_k^*/2$ by $\bar{r}_k - 1$. Also, both the Hausdorff distance and the Wasserstein distance, can be bounded as in Equation (16) if we replace $r_k^*/2$ by $\bar{r}_k - 1$.

As already stated, we could have obtained a similar result (but with different constants) using the dyadic grid $\mathcal{G}_D$ instead of $\mathcal{G}_F$. To conclude this section, let us compare the conditions (20) and (8) for high-energy. Define

$$\psi_{n,r,s}^{(sg)} = \sqrt{p\gamma_r} \wedge \left(s \log \left(\frac{ep}{s}\right)\right) + \gamma_r,$$

where we recall that $\gamma_r = \log(\frac{n}{r\delta})$. If $\gamma_r \ge p/2$, then $\psi_{n,r,s}^{(sg)} \asymp \gamma_r$. In low dimension, the energy threshold for multivariate change-point detection is the same as in the univariate setting, see [38]. If $\gamma_r \le p/2$, then

$$\psi_{n,r,s}^{(sg)} \asymp \begin{cases} \gamma_r & \text{if } s \le \frac{\gamma_r}{\log(p)-\log(\gamma_r)} \\ s \log\left(\frac{e^p}{s}\right) & \text{if } \frac{\gamma_r}{\log(p)-\log(\gamma_r)} < s < \frac{\sqrt{p\gamma_r}}{\log(p)-\log(\gamma_r)} \\ \sqrt{p\gamma_r} & \text{if } s \ge \frac{\sqrt{p\gamma_r}}{\log(p)-\log(\gamma_r)} \end{cases}$$

As a consequence, $\psi_{n,r,s}^{(sg)}$ and $\psi_{n,r,s}^{(g)}$ are of the same order of magnitude for all $s$ when $\gamma_r \ge p/2$. When $\log(n/r\delta) < p$, they are also of the same order of magnitude except when $s$ is close but smaller than $\sqrt{p\gamma_r}$, for which the ratio $\psi_{n,r,s}^{(sg)}/\psi_{n,r,s}^{(g)}$ between these two quantities can be as large as $\log(p) - \log(\gamma_r)$. This gap corresponds to the regime where the test based on the Berk-Jones statistic defined in Equation (14), used in the Gaussian case, outperforms the test based on the partial CUSUM norm statistic defined in Equation (15).

In the definitions of the tests, the tuning constants $\bar{c}_{\text{thresh}}^{(p)}$ and $\bar{c}_{\text{thresh}}^{(d)}$ are left implicit, although one can find suitable values by following the proofs of Propositions 3 and 4. In practice, the practitioner can calibrate them by a Monte-Carlo method by simulating a Gaussian multivariate times series without any change-points. Then, $\bar{c}_{\text{thresh}}^{(p)}$ and $\bar{c}_{\text{thresh}}^{(d)}$ are chosen so that the Family-wise error rate (FWER) of the two collections ($T_{l,r}^{(d)}$ and $T_{l,r}^{(p)}$ is equal to $\delta$.
where $\tilde{\mathbf{C}}_{l,r}(Y)$ and $\tilde{\mathbf{C}}'_{l,r}(Y)$ are independent. If there is one change-point at position $l$ and no other change-points in $(l-r, l+r)$, then these statistics are identically distributed and we consider $\tilde{\Psi}_{l,r}^{(d)} = \langle \tilde{\mathbf{C}}_{l,r}(Y), \tilde{\mathbf{C}}'_{l,r}(Y) \rangle$ whose expectation is null when there are no change-points in the segment. As a consequence, $\tilde{\Psi}_{l,r}^{(d)}$ does not require the knowledge of $\sigma$; only an upper bound of $\sigma$ is required to calibrate the corresponding test. Such a U-statistics has already been introduced in [45] and analyzed in an asymptotic setting. Unfortunately, since we can only consider even $r$, this precludes us to detecting change-points that are very close together with $r_k = 1$.

In the general case where there is spatial covariance in the noise, that is $\text{var}(\epsilon_t) = \Sigma$ for an unknown but general $\Sigma$, we can still use the same U-statistic described in the previous paragraph for the dense case. For the sparse case, one could use the supremum norm of the CUSUM statistics as in Jirak [21] and Yu and Chen [48]. To calibrate those tests, we need to estimate both the Frobenius and the operator norm of $\Sigma$, which seems to be doable as soon as there are not too many change-points. If the spatial covariance matrix $\text{var}(\epsilon_t)$ is unknown and even allowed to change with time, we suspect that the problem becomes intrinsically more involved.

## 7.2 Optimal Localization rates

In this work, we mainly considered the problem of **detecting** change-points in the mean of a random vector. We provided tight conditions on the energy so that a change-point is detectable. When such a change-point $\tau_k$ is detected, Corollary 2 states that its position is estimated up to an error of $r_k^*$, which is also of the order of $\sigma^2 \Psi_{n, \tau_k, s_k}^{(g)} \Delta_k^{-2}$ – see the definition (9). It is not clear whether this error is optimal or not.

In the univariate setting ($p=1$), [38] has established that, above the detection threshold, a specific change-point position $\tau_k$ can be localized at the rate $\sigma^2 \Delta_k^{-2}$. In the multivariate setting, the situation is more tricky and there are certainly several localization regimes beyond the detection threshold. It is an interesting direction of research to pinpoint the exact localization rate between $\sigma^2 \Delta_k^{-2}$ and $\sigma^2 \Psi_{n, \tau_k, s_k}^{(g)} \Delta_k^{-2}$. We leave this for future work.

## 7.3 On the choice of the grid in the generic algorithm

Our general procedure is defined for almost any arbitrary grid. Optimal procedures with the dyadic grid are introduced in Sections 3 and 6, whereas we use a near-optimal procedure on the complete grid in Section 4.

From a computational perspective, the procedure’s worst-case complexity is proportional to the size $|\mathcal{G}|$ of the grid $\mathcal{G}$. In that respect, the dyadic grid and more generally the *a*-adic grids benefit from a linear size whereas the size of the complete grid is quadratic.

From a mathematical perspective, it is much easier to control the behaviour of the procedure for an *a*-adic grid by a simple Bonferroni correction on all the statistics as it turns out that this correction is sufficient for our purpose – see the proofs of Section 3. In constrast, controlling larger collections of tests turns out to be much more challenging as one needs to carefully take into account the dependences between the test statistics, which becomes all the more challenging for complex models. As an example, we introduced in Section 3 Berk-Jones statistics to achieve the tight minimax condition for change-point detection. Unfortunately, we did not manage to apply a suitable chaining argument to these statistics and were therefore unable to control the behavior of the corresponding change-point detection procedure on the complete grid.

From a purely statistical perspective, it is difficult to appreciate the respective benefits of denser or sparser grids. On the one hand, for denser grids, the approximation $\bar{\tau}_k$ of $\tau_k$ at scale $r$ will be
**Computational Cost.** The computational cost of the statistic $T_{l,r} = T_{l,r}^{(d)} \lor T_{l,r}^{(p)}$ is $O(p(r + \log(p)))$. Thus, a naive computation of all the tests $T_{l,r}$ for $(l, r)$ in the complete grid $\mathcal{G}_F$ requires $O(p \log(p) \sum_{(l,r) \in \mathcal{G}_F} r) = O(pn(n^2 + \log(p)))$ operations. Nevertheless, using the fact that $\sum_{i=l+1}^{l+r} Y_i = (\sum_{i=1}^{l+r-1} Y_i) + Y_{l+r} - Y_l$, it is possible to compute all the tests at scale $r$ with cost $O(np \log(p))$. Since there are $n$ possible scales $r$ on the complete grid, the whole procedure cost is $O(n^2p \log(p))$. Using a grid $\mathcal{G} = \{(l, r) \in \mathcal{G}_F : r \in \mathbb{R}\}$ that contains dyadic scales and all possible locations $l$ for each scale, the whole change-point detection would then require only $O(np \log(n) \log(p))$ computations, since there are only $\log(n)$ possible scales $r$ for such grids.

## 5. Minimax lower bound

In this section, we write for any $\Theta \in \mathbb{R}^{p \times n}$, the distribution of the time series $Y = (y_1, \dots, y_n)$ in the model (2) with Gaussian noise $\varepsilon_t \sim \mathcal{N}(0, \sigma^2 \mathbf{I}_p)$. In Section 3, we have established that any change-point satisfying the condition (8), that is

$$r_k \Delta_k^2 \geq c_0 \sigma^2 \left[ s_k \log \left( 1 + \frac{\sqrt{p}}{s_k} \sqrt{\log \left( \frac{n}{r_k \delta} \right)} \right) + \log \left( \frac{n}{r_k \delta} \right) \right],$$

is detected by our change-point procedure. We now show that this energy condition is unimprovable from a minimax point of view. More precisely, let us define, for any $u > 0$, the class $\tilde{P}(u)$ of mean parameters $\Theta$ with arbitrary $K \geq 0$ number of change points and such that any change-point $\tau_k$ for $1 \leq k \leq K$ satisfies

$$r_k \Delta_k^2 \geq \frac{1}{2} \sigma^2 \left[ s_k \log \left( 1 + u \frac{\sqrt{p}}{s_k} \sqrt{\log \left( \frac{n}{r_k} \right)} \right) + u \log \left( \frac{n}{r_k} \right) \right]. \quad (22)$$

For $u$ small enough, it turns out no change-point estimator is able to detect all change-points without estimating any spurious change-point with high probability on the full class $\tilde{P}(u)$. Still, using this large class provides somewhat pessimistic bounds. For instance, the most challenging distributions in $\tilde{P}(u)$ for the purpose of change-point detection satisfy $s_k = p$ and $r_k = 1$ (very close change-points). As a consequence, relying on the full collection $\tilde{P}(u)$ turns too pessimistic. To establish that our bounds are adaptive with respect to the sparsity $s_k$ and the length $r_k$, we define, for any positive integers $1 \leq r \leq \lfloor n/2 \rfloor$ and any $1 \leq s \leq p$ the collection

$$\tilde{P}(u, r, s) = \{ \Theta \in \tilde{P}(u) : \min_k r_k \geq r \text{ and } \max_k s_k \leq s \}.$$

By convention, constant means $\Theta$ with no change-points ($K = 0$) also belong to $\tilde{P}(u, r, s)$. In the class $\tilde{P}(u, r, s)$, all change-points have a sparsity at most $s$ and a length at least $r$. Hence, $\tilde{P}(u, r, s)$ becomes larger when $s$ increases or when $r$ increases.
closer to $\tau_k$ so that the corresponding test $T_{\bar{\tau}_k, r}$ may be more powerful. On the other hand, for a denser grid, the tests possibly suffer from a higher price for multiplicity. This price can be mild if one takes into account the dependences between the tests. Still, except perhaps in the univariate Gaussian change-point model for which delicate controls of the CUSUM process exist, it is challenging to provide theoretical guidance towards the best choice of the grid.

## 7.4 Optimality of the generic algorithm in a broader context.

Algorithm 1 aggregates homogeneity tests and provides theoretical guarantees on the event $A(T, K^*, (\bar{\tau}_k, \bar{r}_k)_{k \in K^*})$ - i.e. the event where the outcomes of the tests are consistent - as stated in Theorem 1. In the possibly sparse high-dimensional mean change-point model, we introduced a suitable multiple testing procedure which, when combined with Algorithm 1, leads to a minimax optimal change-point detection procedure.

We described in Section 2 how to adapt this approach to other change-point problems and this was already illustrated in Section 6 with covariance and nonparametric problems. One may then wonder whether this roadmap still leads to minimax optimal procedures for general problems. Consider the general setting from Section 1 where we are interested in detecting change-points in $(\Gamma(\mathbb{P}_t))_{t \in [n]}$. Upon endowing the space $\mathcal{V}$ with some distance $d$, we define, for any $k$,

$$ \bar{\Delta}_k = d(\Gamma(\mathbb{P}_{\tau_k}), \Gamma(\mathbb{P}_{\tau_{k-1}})), $$

which corresponds to the change-point height. Then, one may wonder how large $\bar{\Delta}_k$ has to be - as a function of $r_k$ - so that a change-point detection procedure achieving the no-spurious property (NoSp) with high probability is able to detect $\tau_k$. In this discussion, we restrict our attention to independent observations, that is the random variables $y_t$ are assumed to be independent and we consider the dyadic grid $\mathcal{G}_D$.

Fix $\delta \in (0, 1)$. At each scale $r \in \{1, 2, \dots, 2^{\lfloor \log_2(n) \rfloor - 1}\}$ and for each $l \in \mathcal{D}_r$, with $\mathcal{D}_r$ defined in (5), we consider the testing problem $H_{0,l,r} : \{\mathbb{P} : \Gamma(\mathbb{P}_{l-r}) = \dots = \Gamma(\mathbb{P}_{l+r-1})\}$ versus

$$ H_{\rho,l,r} : \left\{ \mathbb{P} : \begin{array}{l} \Gamma(\mathbb{P}_{l-r}) = \dots = \Gamma(\mathbb{P}_{l-m-1}) \\ \Gamma(\mathbb{P}_{l-m}) = \dots = \Gamma(\mathbb{P}_{l+r-1}) \quad \text{for some integer } m \in [-r/2, r/2] \\ d(\Gamma(\mathbb{P}_{l-m-1}), \Gamma(\mathbb{P}_{l-m}) \ge \rho) \end{array} \right\} $$

This amounts to testing whether there is a single change-point near $l$ of height at least $\rho$ in the segment $(l-r, l+r)$. Given $\delta \in (0, 1)$ and a test $T$ we define the $\delta$-separation distance of $T$ by

$$ \rho_{l,r}^*(T, \delta) = \inf \left\{ \rho : \sup_{\mathbb{P} \in H_{0,l,r}} \mathbb{P}(T=1) \lor \sup_{\mathbb{P} \in H_{\rho,l,r}} \mathbb{P}(T=0) \le \delta \right\}. $$

This corresponds to the minimal change-point height that is detected by the test $T$. Then, the minimax separation distance $\rho_{l,r}^*(\delta)$ is simply $\inf_T \rho_{l,r}(T, \delta)$, i.e. the infimum over all tests $T$ of the separation distance. By translation invariance of the testing problem, note that $\rho_{l,r}^*(\delta)$ does not depend on $l$ and is henceforth denoted $\rho_r^*(\delta)$.

For any $(l, r)$, take any test $T_{l,r}$ (nearly)<fnref n="3" /> achieving the minimax separation distance $\rho_r^*(\delta | \mathcal{D}_r |^{-1}\beta_r)$ with $\beta_r = 6 \log_2^{-2}(n/r)\pi^{-2}$. Then, it follows from a simple union bound on the

<footnote n="3">³Since the minimax separation distance is defined as an infimum, it is not necessarily achieved by a test. Still, we can build a test whose separation distance is arbitrarily close to the optimal one. We neglect the additive error term for the purpose of the discussion.</footnote>
**Theorem 2.** Fix any $u \in (0, 1/8)$. For any $\sigma > 0$, $n \ge 2$, $p \ge 1$, any length $1 \le r \le n/4$, and any sparsity $1 \le s \le p$, we have

$$ \inf_{\hat{\tau}} \sup_{\Theta \in \bar{P}(u,r,s)} \mathbb{P}_{\Theta}(\hat{K} \neq K) \geq \frac{1}{4}, $$

where the infimum is taken over all estimators $\hat{\tau}$ of the change-point vector $\tau$ and $\hat{K} = |\hat{\tau}|$.

Thus, in the Gaussian setting, if all the change-points have a high-energy in the sense of (8) but with a smaller multiplicative constant factor, no change-point estimator can consistently estimate the true number of change-points. The next corollary restates this negative results in the same lines as Corollary 3.

**Corollary 6.** Fix any $u \in (0, 1/8)$. For any $\sigma > 0$, $n \ge 2$, $p > 1$, any length $1 \le r \le n/4$, any sparsity $1 \le s \le p$, and any estimator $\hat{\tau}$, there exists some $\Theta \in \bar{P}(u, r, s)$ such that with $\mathbb{P}_{\Theta}$-probability larger than $1/4$, at least one of the two following properties is satisfied

* $\hat{\tau}$ contains at least one spurious change-point

* at least a change-point $\tau_k$ with $1 \le k \le K$ is not detected, i.e. there is no change-point estimated in the interval $[(\tau_{k-1} + \tau_k)/2, (\tau_k + \tau_{k+1})/2]$.

This corollary is to be compared to Corollary 3 – indeed, the energy condition in Equation (22) differs from Equation (8) only by a numerical multiplicative constant. As a consequence, the energy condition (22) is minimal for detection by a change-point estimator that achieves (NoSp).

# 6. Application to other change-point problems

In this section, we apply the general methodology of Section 2 to two other problems, namely detection of covariance and nonparametric change-points. This allows us to obtain the first tight minimax detection conditions for these problems.

## 6.1. Covariance change-point detection

Following Wang et al. [40], we consider the covariance change-point model where the covariance matrices $\Sigma_t$ of the centered random vectors $y_t \in \mathbb{R}^p$ are piecewise constant. Then, the goal is to estimate the times $0 < \tau_1 < \dots < \tau_K < \tau_{K+1} = n+1$ such that $\Sigma_t$ is varying. See [40] for motivations. As in that work, we assume that the random vectors $y_t$ are independent and are sub-Gaussian with a uniformly bounded Orlicz norm, that is $\max_{t=1,\dots,n} \|y_t\|_{\psi_2} \le B$ for some known fixed $B$. The Orlicz norm of a random vector $y$ is the supremum of the Orlicz norm of any uni-dimensional projection of $y$ – see e.g. [37]. If the $y_t$'s follow a normal distribution, this amounts to assuming that $\max_{t=1,\dots,n} \|\Sigma_t\|_{op} \le 2B^2$ where $\|\cdot\|_{op}$ is for the operator norm. The purpose of Wang et al. was to detect
dyadic grid that, with probability higher than $1 - \delta$, the collection of tests $T_{l,r}$, where $(l,r)$ belongs to the dyadic grid, does not detect any false positive and detects any change-point $\tau_k$ such that $\bar{\Delta}_k$ is higher than $\rho_{\tilde{r}_k}^*(\delta|\mathcal{D}_{\tilde{r}_k}|^{-1}\beta_{\tilde{r}_k})$, where $\tilde{r}_k$ is the largest scale in $\mathcal{R}$ such that $4(\tilde{r}_k - 1) \le r_k$. As a consequence of Theorem 1, the corresponding detection procedure achieves, with probability higher than $1 - \delta$, the property (NoSp) and detects any change-point satisfying the energy condition $\bar{\Delta}_k \ge \rho_{\tilde{r}_k}^*(r\delta\beta_r/2n)$.

Conversely, we believe that this energy condition is almost tight. Indeed, fix any even range $r \ge 2$. To simplify the discussion suppose that $n/(2r)$ is an integer. We consider a specific instance of the problem where the statistician knows that there are $n/(2r) - 1$ evenly-spaced change-points respectively at $2r+1, 4r+1, \dots, n-2r+1$ that allow to reduce the change-point detection problem to $n/(2r)$ change-point detection problem in intervals $(l-r, l+r]$ for $l = r+1, 3r+1, 5r+1, \dots$. Furthermore, it is known that, in each such segment, there exists at most one change-point that is situated in $[l-0.5r, l+0.5r]$, and if the change-point is present then its height is at least $\rho = \rho_r^*(\delta) - \zeta$ for $\zeta$ arbitrarily small. Since all $n/(2r) - 1$ evenly-spaced change-points $2r+1, 4r+1, \dots, n-2r+1$ are known to the statistician, detecting all remaining change-points is equivalent to building an $n/(2r)$ multiple test of the hypotheses $H_{0,l,r}$ versus $H_{\rho,l,r}$ for $l = r+1, 3r+1, 5r+1, \dots$. If a change-point procedure achieves (NoSp) and detects all change-points with radius at least $r/2$ and height at least $\rho$ with probability at least $1 - \delta$, then one is able, with probability uniformly higher than $1 - \delta$, to simultaneously perform without error $n/(2r)$ independent tests $H_{0,l,r}$ versus $H_{\rho,l,r}$. Since any single test must endure an error with probability at least $\delta$ in the worst case, no collection of independents tests is able to endure less than $1 - (1 - \delta)^{n/(2r)}$. When $n/r$ is large and $\delta < 2r/n$, the latter is of the order of $\delta 2r/n$. Based on this, we conjecture that no change-point procedure is able to achieve, with probability higher than $1 - \delta$ the property (NoSp), and also to detect all change-points with radius at least $r/2$ and height at least $\rho_r^*(2r\delta/n) - \zeta$ for $\zeta > 0$ arbitrarily small.

Comparing the performances of our procedure with the negative arguments that we just outlined, we see that aggregating optimal tests on a dyadic grid allows to detect change-points with (almost) uniform height higher $\rho_{\tilde{r}_k}^*(r_k\delta\beta_{\tilde{r}_k}/(2n))$ whereas, as explained above, we conjecture that a change-point $\tau_k$ can be detected only if $\bar{\Delta}_k \ge \rho_{\tilde{r}_k}^*(2r_k\delta/n)$. Since $\tilde{r}_k \ge (r_k/8) \vee 1-$ as we considered the dyadic grid when constructing $\tilde{r}_k$ - the difference between these two bounds is mostly due to the term $\beta_r$ which is of the order of $\log^2(n/r)$. Whereas it is possible to detect change-points at a given scale with a test of type I error probability $2r\delta/n$, our multi-scale procedure relies on a collection of single tests with type I error probability of the order of $r\delta/n/\log^2(n/r)$. This mild mismatch - that we introduce to deal with the multiplicity of scales - of order $\log^2(n/r)$ is harmless for the Gaussian mean-detection problem. Indeed, one may deduce from our analysis in Section 3 that $\rho_{\tilde{r}_k}^*(2r_k\delta/n)$ is of the same order as $\rho_{\tilde{r}_k}^*(\delta|\mathcal{D}_{\tilde{r}_k}|^{-1}\beta_{\tilde{r}_k})$.

In conclusion, one can build through Algorithm 1 an almost optimal change-point procedure in any model provided that we are given optimal homogeneity tests of the form $H_{0,l,r}$ versus $H_{\rho,l,r}$. This provides a universal reduction of the problem of change-point detection to the problem of homogeneity testing.

# 8 Numerical Experiments

In this section, we illustrate the behavior of our procedure to detect change-points in a sparse high-dimensional setting (2).
small changes in operator norm, that is detecting instants $\tau_k$ such that $\Sigma_{\tau_k} \neq \Sigma_{\tau_k-1}$ with $\|\Sigma_{\tau_k} - \Sigma_{\tau_{k-1}}\|_{op}$ possibly small. Apart from the operator norm, other norms have also been considered e.g. in [10]. Here, we focus on the operator norm as in [40].

Recalling the generic procedure introduced in Section 2, we consider the dyadic grid $\mathcal{G}_D$ and some $\delta \in (0, 1)$. For any $(l, r) \in \mathcal{G}$, we respectively write $\hat{\Sigma}_{l,-r}$ and $\hat{\Sigma}_{l,r}$ for the empirical covariance matrices

$$ \hat{\Sigma}_{l,-r} = r^{-1} \sum_{t=l-r}^{l-1} y_t y_t^T ; \quad \hat{\Sigma}_{l,r} = r^{-1} \sum_{t=l}^{l+r-1} y_t y_t^T . $$

Then, we consider the test $T_{l,r}$ rejecting for large values of $\|\hat{\Sigma}_{l,r} - \hat{\Sigma}_{l,-r}\|_{op}$.

$$ T_{l,r} = 1 \left\{ \|\hat{\Sigma}_{l,r} - \hat{\Sigma}_{l,-r}\|_{op} \ge c_0 B^2 \left[ \sqrt{\frac{p}{r}} + \frac{p}{r} + \sqrt{\frac{\log(\frac{2n}{\delta r})}{r}} + \frac{\log(\frac{2n}{\delta r})}{r} \right] \right\}, \quad (23) $$

where the numerical tuning constant $c_0$ is set in the proof of the following proposition. Relying on concentration bounds [23] for the empirical covariance matrix of sub-Gaussian random vectors, we easily prove that the FWER of the multiple testing procedure $(T_{l,r})$ with $(l, r) \in \mathcal{G}_D$ is small. Then, we can analyze the type II error probability and plug it into the generic result (Theorem 1) to control the behavior of the change-point estimator $\hat{\tau}$. This leads us to the following result. In the sequel, a change-point $\tau_k$ is said to have a high-energy if

$$ r_k \|\Sigma_{\tau_k} - \Sigma_{\tau_{k-1}}\|_{op}^2 \ge c_1 B^4 \left[ \left( p + \log \left( \frac{2n}{r_k \delta} \right) \right) \wedge r_k \right], \quad (24) $$

where the numerical constant $c_1$ is introduced in the proof of the following proposition. We recall that, by definition of the model, we have $\|\Sigma_{\tau_k} - \Sigma_{\tau_{k-1}}\|_{op} \le 4B^2$.

**Proposition 5.** There exist positive numerical constants $c_0$, $c_1$, and $c_2$ such that the following holds for any $B > 0$ and any sequence of independent centered random vectors $(y_t)$ satisfying $\max_t \|y_t\|_{\psi_2} \le B$. With probability higher than $1-\delta$, the change-point estimator $\hat{\tau}$ satisfies (NoSp) and detects all high-energy change-points in the sense of (24). Besides, any such high-energy change-point $\tau_k$ satisfies

$$ d_{H,1}(\hat{\tau}, \tau_k) \le c_2 B^4 \frac{p + \log (2\delta^{-1}B^{-4}n \|\Sigma_{\tau_k} - \Sigma_{\tau_{k-1}}\|_{op}^2)}{\|\Sigma_{\tau_k} - \Sigma_{\tau_{k-1}}\|_{op}^2} \le \frac{r_k}{4}, \quad (25) $$

under the same event of probability than $1-\delta$.

Let us compare our condition (24) for detection with Theorem 2 in Wang et al. [40]. The authors assume that all the change-points satisfy

$$ \min_k r_k \min_k \|\Sigma_{\tau_k} - \Sigma_{\tau_{k-1}}\|_{op}^2 \ge c_1 B^4 p \log(n) . $$
**Performance Measure.** To assess the quality of change-point estimator $\hat{\tau}$, we first measure whether the estimated number of change-points $\hat{K} = |\hat{\tau}|$ is equal to the true number $K$ of change-points. We also define the **SAND** loss as the proportion of **Spurious estimated change-points** And true change-points that are **Not Detected**:

$$ \text{SAND}((\tau_k), (\hat{\tau}_{k'})) = \frac{1}{K} \sum_{k=1}^{K} \left| \left| \left(\frac{\tau_k + \tau_{k-1}}{2}, \frac{\tau_k + \tau_{k+1}}{2}\right) \right| \cap \left\{ \hat{\tau}_k, k \in [\hat{K}] \right\} \right| - 1 \right| . $$

**Change-point Detection Methods.** In the experiments, we implemented the bottom-up aggregation procedure Algorithm 1 with partial norm tests $T^{(p)}$ and dense test $T^{(d)}$ corresponding to Section 4 on a semi-complete grid $\mathcal{G}_F = \{(l, r) : l \in \{r+1, \dots, n-r+1, r \in \mathbb{R}\} - \text{we take scales } r \text{ in the dyadic set for computational purposes. On a location } l \text{ and a scale } r, \text{ each test statistic can be seen as a partial norm test relying on the statistic } \Psi_{l,r,s}^{(p)} \text{ defined in Section 4.2 and a threshold Threshold}(r, s) \text{ which is either equal to } x_r^{(d)} \text{ when } s=d \text{ - see Section 4.1 - or to } x_{r,s}^{(p)} \text{ when } s \in \mathcal{Z}_r := \{1, 2, 4, \dots, 2^{\lceil \log_2(s_{\max}) \rceil}\} \text{ with } s_{\max} := \frac{\sqrt{p\gamma_r}}{\log(p)-\log(\gamma_r)} \text{ - see Section 4.3 for the definition of the boundary between sparse and dense regimes } s_{\max}. \text{ We actually do not use the definition of } x_r^{(d)} \text{ and } x_{r,s}^{(p)} \text{ for our thresholds Threshold}(r, s) \text{ since they rely on constants that are not necessarily tight, but we rather calibrate them by a Monte-Carlo method using 10.000 independent samples. For each sample consisting in a time series made of } n \text{ gaussian normal centered vector in } \mathbb{R}^p, \text{ and for each } r \in \mathbb{R}, s \in \mathcal{Z}_r \cup \{p\}, \text{ we compute the maximum over all } l \text{ of the statistics } \Psi_{l,r,s}^{(p)}. \text{ Considering the list of all the 10.000 maximums and taking } \delta = 5\%, \text{ Threshold}(r, s) \text{ is then defined as the } (1-\delta/(2|\mathcal{R}||\mathcal{Z}_r|))\text{-quantile if } s \in \mathcal{Z}_r \text{ and as the } (1-\delta/(2|\mathcal{R}|))\text{-quantile if } s = p, \text{ so that, by a union bound, the total probability of finding a false positive is less than } \delta. \text{ Note that this calibration step only depends on } n, p, \text{ and } \sigma \text{ and only needs to be performed once and for all.}$$

We compare our procedure with the inspect method of [46] which is available as an R package. The tuning parameters of inspect are computed with the automatic method defined in the same R package.

In all the following experiments, we fix the dimension $p = 100$ and the sample size $n = 200$. We generate a piecewise constant signal $(\eta_t)_{t=1}^n$ in $\mathbb{R}^p$ with possible change-points $(\tau_1, \dots, \tau_K)$ using one of the three following settings. We then add a scaling factor $\alpha > 0$ and apply our procedure to the data $y_t = \alpha\eta_t + \varepsilon_t$, which amounts to setting $\theta_t = \alpha\eta_t$ in model (2). We fix the variance of all the coordinates of $\varepsilon_t$ to be equal to one. Increasing $\alpha$ on a grid with step 0.1 allows us to experimentally identify a transition between the regime where we do not detect precisely the change-points - in which case the two losses tend to be close to one - and the regime where we do detect the change-points - in which cases the losses are smaller. We consider three simulation settings:

1. **Segment.** We generate a signal $\eta$ which is zero everywhere, except on $[80, 100]$ where we set it equal to a random vector $\Delta$ with $\|\Delta\| = 1$ and $\|\Delta\|_0 = s$, for $s = 1, 20, 100$. In each one of these cases, we choose the location of the $s$ non null coordinates of $\Delta$ uniformly at random and their value uniformly at random in the set $\{-1/\sqrt{s}, 1/\sqrt{s}\}$. Each time, $\eta$ has 2 true change-points, and we generate the noise $(\varepsilon_t)$ as independent centered and normalized gaussian vectors.

2. **Multiple Change points.** We generate 10 uniform random locations $\tau_1 < \tau_2 < \dots < \tau_{10}$ on $[1, 200]$. For each location $\tau_i$, we generate a uniform random integer $s_i \in [1, 100]$ and a vector $\Delta_i$ as in the segment setting with $\|\Delta_i\| = 1$ and $\|\Delta_i\|_0 = s_i$. We generate a uniform