Robust Inference for Change Points in High Dimension

High-dimensional data analysis often encounters testing and estimation of change-points in the mean, and it has attracted a lot of attention in statistics recently. See Horváth and Hušková, (2012), Jirak, (2015), Cho, (2016), Wang and Samworth, (2018), Liu et al., (2020), Wang et al., (2022), Zhang et al., (2021) and Yu and Chen, (2021, 2022) for some recent literature. Among the proposed tests and estimation methods, most of them require quite strong moment conditions (e.g., Gaussian or sub-Gaussian assumption, or sixth moment assumption) and some of them also require weak component-wise dependence assumption. There are only a few exceptions, such as Yu and Chen, (2022), where they used anti-symmetric and nonlinear kernels in a U-statistics framework to achieve robustness. However, the limiting distribution of their test statistic is non-pivotal and their procedure requires bootstrap calibration, which could be computationally demanding. In addition, their test statistic targets the sparse alternative only. As pointed out in athe review paper by Liu et al., (2022), the interest in the dense alternative can be well motivated by real data and is often the type of alternative the practitioners want to detect. For example, copy number variations in cancer cells are commonly manifested as change-points occurring at the same positions across many related data sequences corresponding to cancer samples and biologically related individuals; see Fan and Mackey, (2017).

In this article, we propose a new test for a change point in the mean of high-dimensional data that works for a broad class of data generating processes. In particular, our test targets the dense alternative, is robust to heavy-tailedness, and can accommodate both weak and strong coordinate-wise dependence. Our test is built on two recent advances in high-dimensional testing: spatial sign based two sample test developed in Chakraborty and Chaudhuri, (2017) and U-statistics based change-point test developed in Wang et al., (2022). Spatial sign based tests have been studied in the literature of multivariate data and they are usually used to handle heavy-tailedness, see Oja, (2010) for a book-length review. However, it was until recently that Wang et al., (2015) and Chakraborty and Chaudhuri, (2017) discovered that spatial sign could also help relax the restrictive moment conditions in high dimensional testing problems. In Wang et al., (2022), they advanced the high-dimensional two sample U-statistic pioneered by Chen and Qin, (2010) to the change-point setting by adopting the self-normalization (SN) (Shao,, 2010; Shao and Zhang,, 2010). Their test targets dense alternative, but requires sixth moment assumption and only allows for weak coordinate-wise dependence.

Building on these two recent advances, we shall propose a spatial signed SN-based test for a change point in the mean of high-dimensional data. Our contribution to the literature is threefold. Firstly, we derive the limiting null distribution of our test statistic under the so-called fixed-$n$ asymptotics, where the sample size $n$ is fixed and dimension $p$ grows to infinity. We discovered that the fixed-$n$ asymptotics provide a better approximation to the finite sample distribution when the sample size is small or moderate. We also let $n$ grows to infinity after we derive $n$-dependent asymptotic distribution, and obtain the limit under the sequential asymptotics (Phillips and Moon,, 1999). This type of asymptotics seems new to the high-dimensional change-point literature and may be more broadly adopted in change-point testing and other high-dimensional problems. Secondly, our asymptotic theory covers both scenarios, the weak coordinate-wise dependence via $\rho$ mixing, and strong coordindate-wise dependence under the framework of “randomly scaled $\rho$-mixing sequence” (RSRM) in Chakraborty and Chaudhuri, (2017). The process convergence associated with spatial signed U-process we develop in this paper further facilitates the application of our test under sequential asymptotics where $n$, in addition to $p$, also goes to infinity. In particular, we have developed novel theory to establish the process convergence result under the RSRM framework. In general, this requires to show the finite dimensional convergence and asymptotic equicontinuity (tightness). For the tightness, we derive a bound for the eighth moment of the increment of the sample path based on a conditional argument under the sequential asymptotics, which is new to the literature. Using this new technique, we provide the unconditional limiting null distribution of the test statistic for the fixed-$n$ and growing-$p$ case. This is stronger than the results in Chakraborty and Chaudhuri, (2017) which is a conditional limiting null distribution. Thirdly, We extend our test to estimate multiple changes by combining the p-value based on the fixed-$n$ asymptotics and the seeded binary segmentation (SBS) (Kovács et al.,, 2020). The use of fixed-$n$ asymptotics is especially recommended due to the fact that in these popular generic segmentation algorithms such as WBS (Fryzlewicz,, 2014) and SBS, test statistics over many intervals of small/moderate lengths are calculated and the sequential asymptotics is not accurate in approximating the finite sample distribution, as compared to its fixed-$n$ counterpart. The superiority and robustness of our estimation algorithm is corroborated in a small simulation study.

The rest of the paper is organized as follows. In Section 2, we define the spatial signed SN test. Section 3 studies the asymptotic behavior of the test under both the null and local alternatives. Extensions to estimating multiple change-points are elaborated in Section 4. Numerical studies for both testing and estimation are relegated to Section 5. Section 6 contains a real data example and Section 7 concludes. All proofs with auxiliary lemmas are given in the appendix. Throughout the paper, we denote $\to_{p}$ as the convergence in probability, $\overset{\mathcal{D}}{\rightarrow}$ as the convergence in distribution and $\rightsquigarrow$ as the weak convergence for stochastic processes. The notations $\bm{1}_{d}$ and $\bm{0}_{d}$ are used to represent vectors of dimension $d$ whose entries are all ones and zeros, respectively. For $a,b\in\mathbb{R}$, denote $a\wedge b=\min(a,b)$ and $a\vee b=\max(a,b)$. For a vector $a\in\mathbb{R}^{d}$, $\|a\|$ denotes its Euclidean norm. For a matrix $A$, $\|A\|_{F}$ denotes its Frobenius norm.

Let $\{X_{i}\}_{i=1}^{n}$ be a sequence of i.i.d $\mathbb{R}^{p}$-valued random vectors with mean $0$ and covariance $\Sigma$. We assume that the observed data $\{Y_{i}\}_{i=1}^{n}$ satisfies $Y_{i}=\mu_{i}+X_{i}$, where $\mathbb{E}X_{i}=\bm{0}_{p}$ and $\mu_{i}\in\mathbb{R}^{p}$ is the mean at time $i$. We are interested in the following testing problem:

In (1), under the null, the mean vectors are constant over time while under the alternative, there is one change-point at unknown time point $k^{*}$.

Let $S(X)=X/\|X\|\mathbf{1}(X\neq\mathbf{0})$ denote the spatial sign of a vector $X$. Consider the following spatial signed SN test statistic:

where for $1\leq l\leq k<m\leq n$,

Here, the superscript ^(s) is used to highlight the role of spatial sign plays in constructing the testing statistic. In contrast to Wang et al., (2022), where they introduced the test statistic

with $D(k;l,m)$ and $W_{n}(k;l,m)$ defined in the same way as (3) but without spatial sign.

As pointed out by Wang et al., (2022), the limiting distribution of (properly standardized) $D(k;1,n)$ relies heavily on the covariance (correlation) structure of $Y_{i}$, which is typically unknown in practice. One may replace it with a consistent estimator, and this is indeed adopted in high dimensional one-sample or two-sample testing problems, see, for example, Chen and Qin, (2010) and Chakraborty and Chaudhuri, (2017). Unfortunately, in the context of change-point testing, the unknown location $k^{*}$ makes this method practically unreliable. For our spatial sign based test, the limiting distribution of (properly standardized) $D^{(s)}(k;1,n)$ depends on unknown nuisance parameters that could be a complex functional of the underlying distribution. To this end, following Wang et al., (2022) and Zhang et al., (2021), we propose to adopt the SN technique in Shao and Zhang, (2010) to avoid the consistent estimation of unknown nuisance parameters. SN technique was initially developed in Shao, (2010) and Shao and Zhang, (2010) in the low dimensional time series setting and its main idea is to use an inconsistent variance estimator (i.e. self-normalizer) which is based on recursive subsample test statistic, so that the limiting distribution is pivotal under the null. See Shao, (2015) for a recent review.

We first introduce the concept of $\rho$-mixing, see e.g. Bradley, (1986). Typical $\rho$-mixing sequences include i.i.d sequences, $m$-dependent sequences, stationary strong ARMA processes and many Markov chain models.

In real applications, in addition to change-point testing, another important task is to estimate the number and locations of these change-points. In this section, we assume there are $m\geq 1$ change-points and are denoted by ${\bm{k}}=(k_{1},k_{2},\cdots,k_{m})\subset\{1,2,\cdots,n\}$. A commonly used algorithm for many practitioners would be binary segmentation (BS), where the data segments are recursively split at the maximal points of the test statistics until the null of no change-points is not rejected for each segment. However, as criticized by many researchers, BS tends to miss potential change-points when non-monotonic change patterns are exhibited. Hence, many algorithms have been proposed to overcome this drawback. Among them, wild binary segmentation (WBS) by Fryzlewicz, (2014) and its variants have become increasingly popular because of their easy-to-implement procedures. The main idea of WBS is to perform BS on randomly generated sub-intervals so that some sub-intervals can localize at most one change-point (with high probability). As pointed out by Kovács et al., (2020), WBS relies on randomly generated sub-intervals and different researchers may obtain different estimates. Hence, Kovács et al., (2020) propose seeded binary segmentation (SBS) algorithm based on deterministic construction of these sub-intervals with relatively cheaper computational costs so that results are replicable. To this end, we combine the spatial signed SN test with SBS to achieve the task of multiple change-point estimation, and we call it SBS-SN^(s). We first introduce the concept of seeded sub-intervals.

Let $\alpha\in[1/2,1)$ be a decay parameter, denote $\mathcal{I}_{\alpha}(n)$ as the set of seeded intervals based on Definition 4.1. For each sub-interval $(a,b)\in\mathcal{I}_{\alpha}(n)$, we calculate the spatial signed SN test

where $D^{(s)}(k;a,b))$ and $W_{b-a+1}^{(s)}(k;a,b)$ are defined in (3) and (4). We obtain the p-value of the sub-interval test statistic $T^{(s)}(a,b)$ based on the fixed-$n$ asymptotic distribution $\mathcal{T}_{b-a+1}$. SBS-SN^(s) then finds the smallest p-value evaluated at all sub-intervals and compare it with a predetermined threshold level $\zeta_{p}$. If the smallest p-value is also smaller than $\zeta_{p}$, denote the corresponding sub-interval where the smallest p-value is achived as $(a^{*},b^{*})$ and estimate the change-point by $\hat{k}=\arg\max_{k\in\{a^{*}+3,\cdots,b^{*}-4\}}\frac{(D^{(s)}(k;a^{*},b^{*}))^{2}}{W_{b^{*}-a^{*}+1}^{(s)}(k;a^{*},b^{*})}$. Once a change-point is identified, SBS-SN^(s) then divides the data sample into two subsamples accordingly and apply the same procedure to each of them. The process is implemented recursively until no change-point is detected. Details are provided in Algorithm 1.

Our SBS-SN^(s) algorithm differs from WBS-SN algorithm in Wang et al., (2022) and Zhang et al., (2021) in two aspects. First, WBS-SN is built on WBS, which relies on randomly generated intervals while SBS relies on deterministic intervals. As documented in Kovács et al., (2020), WBS is computationally more demanding than SBS. Second, the threshold used in WBS-SN is universal for each sub-interval, depends on the sample size $n$ and dimension $p$ and needs to be simulated via extensive Monte Carlo simulations. Generally speaking, WBS-SN requires simulating a new threshold each time for a new dataset. By contrast, our estimation procedure is based on p-values under the fixed-$n$ asymptotics, which takes into account the interval length $b-a+1$ for each sub-interval $(a,b)$. When implementing either WBS or SBS, inevitably, there will be intervals of small lengths. Hence, the universal threshold may not be suitable as it does not take into account the effect of different interval lengths. In order to alleviate the problem of multiple testing, we may set a small threshold number for $\zeta_{p}$, such as 0.0001 or 0.0005. Furthermore, the WBS-SN requires to specify a minimal interval lentgh which can affect the finite sample performance. In this work, when generating seed sub-intervals as in Definition 4.1, the lengths of these intervals are set as integer values times 10 to reduce the computational cost for simulating fixed-$n$ asymptotic distribution $\mathcal{T}_{n}$. Therefore, we only require the knowledge of $\{\mathcal{T}_{n}\}_{n=10,20,\cdots}$ for SBS-SN^(s) to work, which can be simulated once for good and do not change with a new dataset.

This section examines the finite sample behavior of the proposed tests and multiple change-point estimation algorithm SBS-SN^(s) via simulation studies.