Valid and Exact Statistical Inference for Multi-dimensional Multiple Change-Points by Selective Inference

Consider a dataset in the form of a multi-dimensional sequence with the length $N$ and the dimension $D$. We denote the dataset as $X\in\mathbb{R}^{D\times N}$ where each row is called a component and each column is called a location. The $i^{\rm th}$ component, denoted as $X_{i,:}\in\mathbb{R}^{N}$ for $i\in[D]$, is represented as one-dimensional sequence with the length $N$. The $j^{\rm th}$ location, denoted as $X_{:,j}\in\mathbb{R}^{D}$ for $j\in[N]$, is represented as a vector with the length $D$. Furthermore, we denote the $ND$-dimensional vector obtained by concatenating rows of $X$ as ${\rm vec}(X)=[X_{1,:}^{\top},\ldots,X_{D,:}^{\top}]^{\top}\in\mathbb{R}^{ND}$ where ${\rm vec}(\cdot)$ is an operator to transform a matrix into a vector with concatenated rows.

We interpret that the dataset $X$ is a random sample from a statistical model

where $M\in\mathbb{R}^{D\times N}$ is the matrix containing the mean values of each element of $X$, $\Xi\in\mathbb{R}^{D\times D}$ is the covariance matrix for representing correlations between components and $\Sigma\in\mathbb{R}^{N\times N}$ is the covariance matrix for representing correlations between locations. Here, we do not assume any specific structure for the mean matrix $M$ and assume it to be unknown. On the other hand, the covariance matrices $\Xi$ and $\Sigma$ are assumed to be estimable from an independent dataset where it is known that no change points exist.

We consider the problem of the identifying both the locations and the components of multiple CPs in the dataset $X$ when there exist changes in a subset of the $D$ components at each of the multiple CP locations. Let $K$ be the number of detected CPs. The locations of the $K$ detected CPs are denoted as $1\leq\tau_{1},\ldots,\tau_{K}<N$ and the subsets of the components where changes present are denoted as $\bm{\theta}^{1},\ldots,\bm{\theta}^{K}$ where $\bm{\theta}^{k}$, $k\in[K]$, is the set of the components in which changes exist at the location $\tau_{k}$. Elements of $\bm{\theta}^{k}$ is denoted as $\theta^{k}_{1},\ldots,\theta^{k}_{|\bm{\theta}^{k}|}$. For notational simplicity, we also define $\tau_{0}=0$ , $\tau_{K+1}=N$ and $\bm{\theta}^{0}=\emptyset$, $\bm{\theta}^{K+1}=\emptyset$. Given a multi-dimensional sequence dataset $X$, the task is to identify both the locations ${\mathcal{T}}=(\tau_{1},\ldots,\tau_{K})$ and the sets of the components $\Theta=(\bm{\theta}^{1},\ldots,\bm{\theta}^{K})$. of the detected multiple CPs (see Fig. 1).

As a CP detection method for both locations and components, we employ a method combining dynamic-programming (DP)-based approach (Bellman, 1966; Guthery, 1974; Bai and Perron, 2003; Lavielle, 2005; Rigaill, 2010) and scan-statistic (SS)-based approach (Enikeeva et al., 2019). This CP detection method is formulated as the following optimization problem:

where

Here, $L$ and $W$ are hyperparameters which are explained below, $\bm{1}_{D,\bm{\theta}}$ is $D$-dimensional binary vector in which the $i^{\rm th}$ element is $1$ if $i\in\bm{\theta}$ and $0$ otherwise, and ${\rm diag}(\bm{1}_{D,\bm{\theta}})\in\mathbb{R}^{D\times D}$ is a matrix in which vector $\bm{1}_{D,\bm{\theta}}$ is arranged diagonally. The CP detection method in (2) is motivated as follows. We employ the scan statistic proposed in (Enikeeva et al., 2019) as the building block, which was introduced for the problem of detecting a single CP from a multi-dimensional sequence¹¹1 Double CUMSUM (Cho, 2016) is known as another CP detection method that enables us to select not only the location but also the subset of the components of the change. We focus on the scan-statistic in this study, but conjecture that our SI framework can be applied also to double CUMSUM. . As a criterion for a single CP at the location $\tau\in[N]$ and the components $\bm{\theta}\subseteq[D]$, the scan statistic is defined as

We extend the scan statistic in various directions in (2). First, we extended the scan statistic to handle multiple CPs, i.e., multiple locations $\tau_{1},\ldots,\tau_{K}$ and the corresponding sets of the components $\bm{\theta}^{1},\cdots,\bm{\theta}^{K}$. Next, in each component $\theta^{k}_{h}\in\bm{\theta}^{k},h\in[|\bm{\theta}^{k}|]$ we consider the mean signal difference between the $L$ points before and after the CP location $\tau_{k}$, i.e., the mean difference between $X_{h,\tau_{k}-L+1:\tau_{k}}$ and $X_{h,\tau_{k}+1:\tau_{k}+L}$. This is to avoid the influence of neighboring CPs when evaluating the statistical reliability of each location $\tau_{k}$ and the components $\theta^{k}_{h}$. Here, $L$ is considered as a hyperparameter Furthermore, the reason why we use $Z_{(s,e)}^{\prime}(\tau)$ instead of $Z_{(s,e)}(\tau)$ in (3a) is that we consider changes to be common among components if it is within the range of $\pm W$ rather than defining common changes as the changes at exactly the same location. Here, $W$ is also considered as a hyperparameter.

The optimal solution of (2) can be efficiently obtained by using DP and the algorithm for efficient scan-statistic computation described in Enikeeva et al. (2019). The pseudo-code for solving the optimization problem (2) is presented in Algorithms 1 (SortComp and UpdateDP functions are presented in Appendix A).

The main goal of this paper is to develop a valid and exact (non-asymptotic) statistical inference method for detected locations and components of multiple CPs. For each of the detected location $\tau_{k}$ and component $\theta^{k}_{h}$, $k\in[K],h\in[|\bm{\theta}^{k}|]$, we consider the statistical test with the following null and alternative hypotheses:

This test simply concerns whether the averages of $L$ points before and after the CP location $\tau_{k}\pm W$ is equal or not in the component $\theta_{h}^{k}$. Note that, throughout the paper, we do NOT assume that the mean of each component $M_{:,i},i\in[D]$ is piecewise-constant function, i.e., the proposed statistical inference method in §3 is valid in the sense that the type I error is properly controlled at the specified significance level $\alpha$ (e.g., $\alpha=0.05$). While existing statistical inference methods for CPs often require the assumption on the true structure (such as piecewise constantness, piecewise linearity, etc.), the conditional SI framework in this paper enables us to be free from assumptions about the true structure on $M_{i,j},i\in[D],j\in[N]$ although the power of the method depend on the unknown structure.

For the statistical test in (4), a reasonable test statistic is the observed mean signal difference

This test statistic can be written as

by defining an $ND$-dimensional vector

where $\bm{1}_{\theta^{k}_{h}}$ is $D$-dimensional binary vector in which the $\theta_{h}^{k}$th element is 1 and 0 otherwise, and $\bm{\mathbf{1}}_{s:e}$ is $N$-dimensional binary vector in which the $j^{\rm th}$ element is 1 if $j\in[s:e]$ and 0 otherwise. In case both the CP location $\tau_{k}$ and the CP component $\theta^{k}_{h}$ are known before looking at data, since ${\rm vec}(X)$ follows the multivariate normal distribution in (1), the null distribution of the test statistic $\bm{\eta}_{k,h}^{\top}{\rm vec}(X)$ also follows a multivariate normal distribution under the null hypothesis (4a). Unfortunately, however, because the CP location $\tau_{k}$ and the CP component $\theta^{k}_{h}$ are actually selected by looking at the observed dataset $X^{\text{obs}}$, it is highly difficult to derive the sampling distribution of the test statistic. In the literature of traditional statistical CP detection, in order to mitigate this difficulty, only simple problem settings, e.g., with a single location or/and a single component, have been studied by relying on asymptotic approximation with $N\to\infty$ or/and $D\to\infty$. However, these asymptotic approaches are not available for complex problem settings such as (2), and are only applicable when certain regularity assumptions on sequential correlation are satisfied.

To conduct exact (non-asymptotic) statistical inference on the locations and the components of the detected CPs, we employ conditional SI. The basic idea of conditional SI is to make inferences based on the sampling distribution conditional on the selected hypothesis (Lee et al., 2016; Fithian et al., 2014). Let us denote the set of pairs of location and component of the detected CPs when the dataset $X$ is applied to Algorithm 1 as

To apply conditional SI framework into our problem, we consider the following sampling distribution of the test statistic $\bm{\eta}_{k,h}^{\top}{\rm vec}(X)$ conditional on the selection event that a change is detected at location $\tau_{k}$ in the component $\theta^{k}_{h}$, i.e.,

Then, for testing the statistical significance of the change, we introduce so-called selective $p$-value $p_{k,h}$, which satisfies the following sampling property

where $\alpha\in(0,1)$ is the significance level of the statistical test (e.g., $\alpha=0.05$)²²2 Note that a valid $p$-value in a statistical test should satisfies $\mathbb{P}(p\leq\alpha)=\alpha~{}\forall\alpha\in(0,1)$ under the null hypothesis so that the probability of the type I error (false positive) coincides with the significance level $\alpha$. . To define the selective $p$-value, with a slight abuse of notation, we make distinction between the random variables $X$ and the observed dataset $X^{\rm obs}$. The selective $p$-value is defined as

where $\bm{q}_{k,h}(X)$ is the nuisance parameter defined as

Note that the selective $p$-values depend on the nuisance component $\bm{q}_{k,h}(X)$, but the sampling property in (5) is kept satisfied without conditioning on the nuisance event $\bm{q}_{k,h}(X)=\bm{q}_{k,h}(X^{\rm obs})$ because $\bm{q}_{k,h}(X)$ is independent of the test statistic $\bm{\eta}_{k,h}^{\top}{\rm vec}(X)$. The selective $(1-\alpha)$ confidence interval (CI) for any $\alpha\in(0,1)$ which satisfies the $(1-\alpha)$ coverage property conditional on the selection event $(\tau_{k},\theta_{h}^{k})\in{\mathcal{M}}(X)$ can be also defined similarly. See Lee et al. (2016) and Fithian et al. (2014) for a more detailed discussion and the proof that the selective $p$-value (8) satisfies the sampling property (5) and how to define selective CIs.

The discussion so far indicates that, if we can compute the conditional probability in (8), we can conduct a valid exact test for each of the detected change at the location $\tau_{k}$ in the component $\theta^{k}_{h}$ for $k\in[K],h\in[|\bm{\theta}^{k}|]$. Most of the existing conditional SI studies consider the case where the hypothesis selection event is represented as an intersection of a set of linear or quadratic inequalities. For example, in the seminal work of conditional SI in Lee et al. (2016), the selective $p$-value and CI computations are conducted by exploiting the fact that the event that some features (and their signs) are selected by Lasso is represented as an intersection of a set of linear inequalities, i.e., the selection event is represented as a polyhedron of sample (data) space. Unfortunately, however, the selection event in (8) is highly complicated and cannot be simply represented an intersection of linear nor quadratic inequalities. We resolve this computational challenge in the next section. Our basic idea is to consider a family of over-conditioning problems such that the computation of the conditional probability of each of the over-conditioning problems can be tractable and then sum them up to obtain the selective $p$-value in (8). This idea is motivated by Liu et al. (2018) and Le Duy and Takeuchi (2021) in which the authors proposed methods to increase the power of conditional SI for Lasso in Lee et al. (2016).