Adaptive Inference for Change Points in High-Dimensional Data

In this subsection, we shall develop our test statistics for one change point alternative, i.e.,

$$\mathcal{H}_{1}:\mu_{1}=\cdots=\mu_{k_{1}}\neq\mu_{k_{1}+1}=\cdots=\mu_{n},$$

where $k_{1}$ is unknown. In Wang et al. (2019), a U-statistic based approach was developed and their test targets at the dense alternative since the power is a monotone function of $\sqrt{n}\|\Delta\|_{2}/\|\Sigma\|_{F}^{1/2}$, where $\Delta$ is the difference between pre-break and post-break means, i.e., $\Delta=\mu_{n}-\mu_{1}$, and $\Sigma$ is the covariance matrix of $X_{i}$. Thus their test may not be powerful if the change in mean is sparse and $\|\Delta\|_{2}$ is small. Note that several tests have been developed to capture sparse alternatives as mentioned in Section 1. In practice, when there is no prior knowledge of the alternative for a given data set at hand, it would be helpful to have a test that can be adaptive to different types and magnitudes of the change. To this end, we shall adopt the $L_{q}$ norm-based approach, as initiated by Xu et al. (2016) and He et al. (2020), and develop a class of test statistics indexed by $q\in 2\mathbb{N}$, and then combine these tests to achieve the adaptivity.

Denote $X_{i}=(X_{i,1},\ldots,X_{i,p})^{T}$. For any positive even number $q\in 2\mathbb{N}$, consider the following two-sample U-statistic of order $(q,q)$,

$$T_{n,q}(k)=\frac{1}{P_{q}^{k}P_{q}^{n-k}}\sum_{l=1}^{p}\sum^{*}_{1\leq i_{1},\ldots,i_{q}\leq k}\sum^{*}_{k+1\leq j_{1},\ldots,j_{q}\leq n}\left(X_{i_{1},l}-X_{j_{1},l}\right)\cdots\left(X_{i_{q},l}-X_{j_{q},l}\right),$$

for any $k=q,\cdots,n-q$. Simple calculation shows that $\mathbb{E}\left[T_{n,q}(k)\right]=0$ for any $k=q,\cdots,n-q$ under the null hypothesis, and $\mathbb{E}\left[T_{n,q}(k_{1})\right]=\|\Delta\|_{q}^{q}$ under the alternative. When $q\in 2\mathbb{N}+1$ (i.e., $q$ is odd) and under the alternative, $\mathbb{E}\left[T_{n,q}(k_{1})\right]=\sum_{j=1}^{p}\delta_{j}^{q}\not=\|\Delta\|_{q}^{q}$ where $\Delta=(\delta_{1},\cdots,\delta_{p})^{T}$. This is the main reason we focus on the statistics corresponding to even $q$s since for an odd $q$, $\sum_{j=1}^{p}\delta_{j}^{q}=0$ does not imply $\Delta=0$.

If the change point location $k_{1}=\lfloor\tau_{1}n\rfloor$, $\tau_{1}\in(0,1)$ is known, then we would use $T_{n,q}(k_{1})$ as our test statistic. As implied by the asymptotic results shown later, we have that under the null,

$$\left(\frac{\tau_{1}(1-\tau_{1})}{n\|\Sigma\|_{q}}\right)^{q/2}\frac{T_{n,q}(k_{1})}{\sqrt{q!}}\stackrel{{\scriptstyle D}}{{\rightarrow}}N(0,1),$$

under suitable moment and weak dependence assumptions on the components of $X_{t}$. In practice, a typical approach is to replace $\|\Sigma\|_{q}$ by a ratio-consistent estimator, which is available for $q=2$ [see Chen and Qin (2010)], but not for general $q\in 2\mathbb{N}$. In practice, the location $k_{1}$ is unknown, which adds additional complexity to the variance estimation and motivates Wang et al. (2019) to use the idea of self-normalization [Shao (2010), Shao and Zhang (2010)] in the case $q=2$. Self-normalization is a nascent inferential method [Lobato (2001), Shao (2010)] that has been developed for low and fixed-dimensional parameter in a low dimensional time series. It uses an inconsistent variance estimator to yield an asymptotically pivotal statistic, and does not involve any tuning parameter or involves less number of tuning parameters compared to traditional procedures. See Shao (2015) for a comprehensive review of recent developments for low dimensional time series. There have been two recent extensions to the high-dimensional setting: Wang and Shao (2019) adopted a one sample U-statistic with trimming and extended self-normalization to inference for the mean of high-dimensional time series; Wang et al. (2019) used a two sample U-statistic and extended the self-normalization (SN)-based change point test in Shao and Zhang (2010) to high-dimensional independent data. Both papers are $L_{2}$ norm based, and this seems to be the first time that a $L_{q}$-norm based approach is extended to high-dimensional setting via self-normalization.

Following Wang et al. (2019), we consider the following self-normalization procedure. Define

$\displaystyle U_{n,q}(k;s,m)$

$\displaystyle=\sum_{l=1}^{p}\sum^{*}_{s\leq i_{1},\ldots,i_{q}\leq k}\sum^{*}_{k+1\leq j_{1},\ldots,j_{q}\leq m}\left(X_{i_{1},l}-X_{j_{1},l}\right)\cdots\left(X_{i_{q},l}-X_{j_{q},l}\right),$

which is an un-normalized version of $T_{n,q}$ applied to the subsample $(X_{s},\cdots,X_{m})$. Let

$$W_{n,q}(k;s,m):=\frac{1}{m-s+1}\sum_{t=s+q-1}^{k-q}U_{n,q}(t;s,k)^{2}+\frac{1}{m-s+1}\sum_{t=k+q}^{m-q}U_{n,q}(t;k+1,m)^{2},$$

The self-normalized statistic is given by

$$\widetilde{T}_{n,q}:=\max_{k=2q,\ldots,n-2q}\frac{U_{n,q}(k;1,n)^{2}}{W_{n,q}(k;1,n)}.$$

Before presenting our main theorem, we need to make the following assumptions.

To derive the limiting null distribution for $\widetilde{T}_{n,q}$, we need to define some useful intermediate processes. Define

	$\displaystyle D_{n,q}(r;[a,b])$	$\displaystyle=U_{n,q}(\lfloor nr\rfloor;\lfloor na\rfloor+1,\lfloor nb\rfloor)$
		$\displaystyle=\sum_{l=1}^{p}\sum^{*}_{\lfloor na\rfloor+1\leq i_{1},\ldots,i_{q}\leq\lfloor nr\rfloor}\sum^{*}_{\lfloor nr\rfloor+1\leq j_{1},\ldots,j_{q}\leq\lfloor nb\rfloor}\left(X_{i_{1},l}-X_{j_{1},l}\right)\cdots\left(X_{i_{q},l}-X_{j_{q},l}\right),$

for any $0\leq a<r<b\leq 1$. Note that under the null, $X_{i}$’s have the same mean. Therefore, we can rewrite $D_{n,q}$ as

	$\displaystyle D_{n,q}(r;[a,b])=$	$\displaystyle\sum_{l=1}^{p}\sum^{*}_{\lfloor na\rfloor+1\leq i_{1},\ldots,i_{q}\leq\lfloor nr\rfloor}\sum^{*}_{\lfloor nr\rfloor+1\leq j_{1},\ldots,j_{q}\leq\lfloor bn\rfloor}\left(X_{i_{1},l}-X_{j_{1},l}\right)\cdots\left(X_{i_{q},l}-X_{j_{q},l}\right)$
	$\displaystyle=$	$\displaystyle\sum_{l=1}^{p}\sum^{*}_{\lfloor na\rfloor+1\leq i_{1},\ldots,i_{q}\leq\lfloor nr\rfloor}\sum^{*}_{\lfloor nr\rfloor+1\leq j_{1},\ldots,j_{q}\leq\lfloor bn\rfloor}\left(Z_{i_{1},l}-Z_{j_{1},l}\right)\cdots\left(Z_{i_{q},l}-Z_{j_{q},l}\right)$
	$\displaystyle=$	$\displaystyle\sum_{c=0}^{q}(-1)^{q-c}\binom{q}{c}P^{\lfloor nr\rfloor-\lfloor na\rfloor-c}_{q-c}P^{\lfloor nb\rfloor-\lfloor nr\rfloor-q+c}_{c}S_{n,q,c}(r;[a,b]).$

In the above expression, considering the summand for each $c=0,1,\ldots,q,$ we can define, for any $0\leq a<r<b\leq 1$,

$$S_{n,q,c}(r;[a,b])=\sum_{l=1}^{p}\sum^{*}_{\lfloor na\rfloor+1\leq i_{1},\cdots,i_{c}\leq\lfloor nr\rfloor}\sum^{*}_{\lfloor nr\rfloor+1\leq j_{1},\cdots,j_{q-c}\leq\lfloor nb\rfloor}\left(\prod_{t=1}^{c}Z_{i_{t},l}\prod_{s=1}^{q-c}Z_{j_{s},l}\right),$$

if $\lfloor nr\rfloor\geq\lfloor na\rfloor+1$ and $\lfloor nb\rfloor\geq\lfloor nr\rfloor+1,$ and 0 otherwise.

For illustration, consider the case when $a_{1}<a_{2}<r_{1}<r_{2}<b_{1}<b_{2}$ and $c_{1}\leq c_{2}$. We have

$$\operatorname{cov}\left(Q_{q,c_{1}}(r_{1};[a_{1},b_{1}]),Q_{q,c_{2}}(r_{2};[a_{2},b_{2}])\right)=\binom{c_{2}}{c_{1}}c_{1}!(q-c_{1})!(r_{1}-a_{2})^{c_{1}}(r_{2}-r_{1})^{c_{2}-c_{1}}(b_{1}-r_{2})^{q-c_{2}},$$

which implies, for example,

$$\operatorname{var}\left[Q_{q,c}(r;[a,b])\right]=c!(q-c)!(r-a)^{c}(b-r)^{q-c}.$$

The proof of Theorem 2.1 is long and is deferred to the supplement.

It can be derived that the $G_{q}(\cdot;[\cdot,\cdot])$ is a Gaussian process with the following covariance structure:

$$\operatorname{var}[G_{q}(r;[a,b])]=\sum^{q}_{c=0}\binom{q}{c}^{2}c!(q-c)!(r-a)^{2q-c}(b-r)^{q+c}=q!(r-a)^{q}(b-r)^{q}(b-a)^{q}.$$

When $r_{1}=r_{2}=r$,

$\displaystyle\operatorname{cov}(G_{q}(r;[a_{1},b_{1}]),G_{q}(r;[a_{2},b_{2}]))=q!(r-A)^{q}(b-r)^{q}(B-a)^{q}.$

When $r_{1}\not=r_{2}$,

$$\operatorname{cov}(G_{q}(r_{1};[a_{1},b_{1}]),G_{q}(r_{2};[a_{2},b_{2}]))=q![(r-A)(b-R)(B-a)-(A-a)(R-r)(B-b)]^{q},$$

where $(r,R,a,A,b,B,c,C)$ is defined in Theorem 2.1. The limiting null distribution $\widetilde{T}_{q}$ is pivotal and its critical values can be simulated as done in Wang et al. (2019) for the case $q=2$. The simulated critical values and their corresponding realizations for $q=2,4,6$ are available upon request. For a practical reason, we did not pursue the larger $q$, such as $q=8,10$, since larger $q$ corresponds to more trimming on the two ends and the finite sample performance when $q=6$ is already very promising for detecting sparse alternatives, see Section 4. An additional difficulty with larger $q$ is the associated computation cost and complexity in its implementation.

Let $I$ be a set of $q\in 2\mathbb{N}$ (e.g. {2,6}). Since $\widetilde{T}_{n,q}$s are asymptotically independent for different $q\in I$ under the null, we can combine their corresponding $p$-values and form an adaptive test. For example, we may use $p_{ada}=\min_{q\in I}p_{q}$, where $p_{q}$ is the $p$-value corresponding to $\widetilde{T}_{n,q}$, as a new statistic. Its $p$-value is equal to $1-(1-p_{ada})^{|I|}$. Suppose we want to perform a level-$\alpha$ test, it is equivalent to conduct tests based on $\widetilde{T}_{n,q},\forall q\in I$ at level $1-(1-\alpha)^{1/|I|}$, and reject the null if one of the statistics exceeds its critical value. Therefore, we only need to compare each $\widetilde{T}_{n,q}$ with its $(1-\alpha)^{1/|I|}$-quantile of the corresponding limiting null distribution.

As we explained before, a smaller $q$ (say $q=2$) tends to have higher power under the dense alternative, which is also the main motivation for the proposed method in Wang et al. (2019). On the contrary, a larger $q$ has a higher power under the sparse alternative, as $\lim_{q\rightarrow\infty}\|\Delta_{n}\|_{q}=\|\Delta_{n}\|_{\infty}$. Therefore, with the adaptive test, we can achieve high power under both dense and sparse alternatives with asymptotic size still equal to $\alpha$. This adaptivity will be confirmed by our asymptotic power analysis presented in Section 2.4 and simulation results presented in Section 4.

Liu et al. (2020) recently studied the detection of a sparse change in the high-dimensional mean vector under the Gaussian assumption as a minimax testing problem. Let $\rho^{2}=\min(k_{1},n-k_{1})\|\Delta_{n}\|_{2}^{2}$. In the fully dense case, i.e., when $\|\Delta_{n}\|_{0}=p$, where $\|\Delta_{n}\|_{0}$ denotes the $L_{0}$ norm, Theorem 8 in Liu et al. (2020) stated that the minimax rate is given by $\rho^{2}\asymp\|\Sigma\|_{F}\sqrt{\log\log(8n)}\vee\|\Sigma\|_{s}\log\log(8n)$. Thus under the assumption that $k_{1}/n=\tau_{1}\in(0,1)$, the $L_{2}$-norm based test in Wang et al. (2019) achieves the rate optimality up to a logarithm factor. Consequently, any adaptive test based on $I$ is rate optimal (up to a logarithm factor) as long as $2\in I$.

In the special case $\Sigma=I_{p}$, the minimax rate is given by

$$\rho^{2}\asymp\left\{\begin{array}[]{lr}\sqrt{p\log\log(8n)}&\text{ if }d\geq\sqrt{p\log\log(8n)}\\ d\log\left(\frac{ep\log\log(8n)}{d^{2}}\right)\vee\log\log(8n)&\text{ if }d<\sqrt{p\log\log(8n)}\end{array}\right..$$

Recall that $\Delta_{n}=\delta\cdot(\bm{1}_{d},\bm{0}_{p-d})^{T}$. In the sparse setting $d=o(\sqrt{p})$ and under the assumptions that $d\asymp p^{-v}$,$v\in(0,1/2)$ and $d>\log\log(8n)$, the minimax rate is $d$ (up to a logarithm factor), which corresponds to $\delta\asymp n^{-1/2}$. Our $L_{q}$-norm based test is not minimax rate optimal since the detection boundary is $(\sqrt{p}/d)^{1/q}n^{-1/2}$, which gets closer to $n^{-1/2}$ as $q\in 2\mathbb{N}$ gets larger. In the dense setting $\sqrt{p}=o(d)$ and under the assumptions that $d\asymp p^{-v}$,$v\in(1/2,1)$ and $d>\sqrt{p\log\log(8n)}$, the minimax rate is $\sqrt{p}$ (up to a logarithm factor), which corresponds to $\delta\asymp p^{1/4}/\sqrt{nd}$. Therefore the $L_{2}$-norm based test in Wang et al. (2019) is again rate optimal (up to a logarithm factor).

In this subsection, we propose to estimate the location of a change point assuming that the data is generated from the following single change-point model,

$$X_{t}=\mu_{1}+\Delta_{n}{\bf 1}(t>k^{*})+Z_{t},~{}t=1,\cdots,n,$$

where $k^{*}=k_{1}=\lfloor\tau^{*}n\rfloor$ is the location of change point. In the literature, it is common to focus on the convergence rate of the estimators of the relative location $\tau^{*}\in(0,1)$, that is, we shall focus on the convergence rate of $\hat{\tau}=\hat{k}/n$, where $\hat{k}$ is an estimator for $k^{*}$.

Given the discussions about size and power properties of the SN-based test statistic in Section 2, it is natural to use the argmax of the test statistic as the estimator for $k^{*}$. That is, we define

$$\hat{k}=\operatorname{argmax}_{k=2q,\ldots,n-2q}\frac{U_{n,q}(k;1,n)^{2}}{W_{n,q}(k;1,n)}.$$

To present the convergence rate for $\hat{\tau}$, we shall introduce the following assumptions.

Let $\gamma_{n,q}=n^{q/2}\|\Delta_{n}\|_{q}^{q}/\|\Sigma\|_{q}^{q/2}$ so $\gamma_{n,2}=n\|\Delta_{n}\|^{2}/\|\Sigma\|_{F}$. We have the following convergence rate of $\hat{\tau}$ for the case $q=2$.