# Sparse change detection in high-dimensional linear regression

Fengnan Gao111Fudan University and Shanghai Center for Mathematical Sciences.
Email<EMAIL_ADDRESS>and Tengyao Wang222London School of Economics and
Political Science. Email<EMAIL_ADDRESS>
(())

###### Abstract

We introduce a new methodology ‘charcoal’ for estimating the location of
sparse changes in high-dimensional linear regression coefficients, without
assuming that those coefficients are individually sparse. The procedure works
by constructing different sketches (projections) of the design matrix at each
time point, where consecutive projection matrices differ in sign in exactly
one column. The sequence of sketched design matrices is then compared against
a single sketched response vector to form a sequence of test statistics whose
behaviour shows a surprising link to the well-known CUSUM statistics of
univariate changepoint analysis. The procedure is computationally attractive,
and strong theoretical guarantees are derived for its estimation accuracy.
Simulations confirm that our methods perform well in extensive settings, and a
real-world application to a large single-cell RNA sequencing dataset showcases
the practical relevance.

##  1 Introduction

The past twenty years have witnessed rapid development of statistical
methodologies for high-dimensional data sets, where the number of variables of
interest is often of comparable or even larger order of magnitude than the
number of observations available. The most prominent example, perhaps, is the
line of work on sparse linear regression, which started from the seminal work
of Tibshirani (1996), and was developed and generalized subsequently by many
others (see, e.g. Fan and Lv (2010); Bühlmann and van de Geer (2011) for a
general overview of this area). In many of these works, the primary focus was
on how to exploit the sparsity of the regression coefficients for their
successful estimation, and to achieve this, a homogeneous data generating
mechanism was often assumed for simplicity of analysis.

However, it is usually unrealistic in large, high-dimensional data sets to
assume that the data generating mechanism holds true throughout. In fact,
heterogeneity is the norm rather than exception in Big Data applications.
Several attempts have been made to handle data heterogeneity in high-
dimensional linear models. For instance, Städler et al. (2010); Yin et al.
(2018); Krishnamurthy et al. (2019) considered the problem of learning mixture
of sparse linear regression, where the regression coefficient vector is
sampled from a small set of sparse regression coefficients. When observations
have a temporal structure, one common way to handle heterogeneity is to break
the sequence of observations into shorter time segments on which the data are
more homogeneous. This line of thinking is the driving force behind the recent
revival of interest in changepoint analysis, which dates back to the early
work of Page (1955), but has seen growing interest recently in high-
dimensional settings, see e.g., Lévy-Leduc and Roueff (2009); Bai (2010);
Zhang et al. (2010); Horváth and Hušková (2012); Cho and Fryzlewicz (2014);
Jirak (2015); Cho (2016); Wang and Samworth (2018); Enikeeva and Klopp (2021);
Liu et al. (2021); Follain et al. (2022); Chen et al. (2022).

In our linear regression setting, such a changepoint setup means that the
sequence of regression coefficient vectors has a piecewise constant structure.
More precisely, for an unknown sequence of changepoints
$0<z_{1}<\cdots<z_{\nu}<n$ (for notational convenience, we also define
$z_{0}:=0$ and $z_{\nu+1}:=n$) and an unknown sequence of a regression
coefficients $(\beta^{(r)}:1\leq r\leq\nu+1)$, we assume that the data
$(x_{t},y_{t})\in\mathbb{R}^{p}\times\mathbb{R}$, $1\leq t\leq n$ are
generated according to the following model

$y_{t}=x_{t}^{\top}\beta_{t}+\epsilon_{t},\quad\text{where
$\beta_{t}=\beta^{(r)}$ for $z_{r-1}<{t}\leq z_{r}$, $1\leq r\leq\nu+1$},$ (1)

and $(\epsilon_{t})_{1\leq t\leq n}$ are the observational errors distributed
as $N(0,\sigma^{2}I_{n})$ conditionally independent of $(x_{t})_{1\leq t\leq
n}$. The goal is to locate the changepoints $z_{1},\ldots,z_{\nu}$ upon
observing the response vector $Y=(y_{1},\dots,y_{n})^{\top}$ and the design
matrix $X=(x_{1},\dots,x_{n})^{\top}$.

Classically, when the dimension $p$ is far smaller than $n$, Bai (1997); Bai
and Perron (1998); Julious (2001) showed that a least-square-based approach
works well in the above changepoint problem, which is equivalent to maximum-
likelihood estimations under Gaussianity assumptions. Specifically, for a
given $\nu$, the maximum likelihood estimator finds the optimal partition of
$\\{1,\ldots,n\\}$ into $\nu+1$ segments such that residual sum of squares
from the least-square fit within each segment is minimized. The least-square
(maximum-likelihood) fit from different choices of $\nu$ can then be compared
using for instance the Bayesian Information Criterion (BIC) to choose the best
$\nu$, which is often solved algorithmically via dynamic programming.

In the high-dimensional setting, the above maximum-likelihood/least-square
approach no longer works. Several works have appeared to analyse such
regression changepoint problems in the high-dimensional context, see for
instance Rinaldo et al. (2021); Wang et al. (2021); Lee et al. (2016); Kaul et
al. (2019) and references therein. However, in addition to the modelling
assumption in (1), these works also impose the additional assumption that all
regression coefficients $(\beta^{(r)}:1\leq r\leq\nu+1)$ are individually
sparse. Given a hypothesized set of changepoints, this additional assumption
allows them to form estimators of $\beta^{(r)}$, $1\leq r\leq\nu+1$, which are
in turn used to form goodness-of-fit statistics for the set of hypothesized
changepoints.

A major difference between this work and the aforementioned existing line of
works is that we do not assume that the regression coefficients within
different stationary segments are individually sparse. Instead, we make the
less stringent assumption that the difference in the regression coefficient
vectors before and after each change, i.e.,
$\theta^{(r)}:=(\beta^{(r+1)}-\beta^{(r)})/2$, are sparse in the sense that
$\|\theta^{(r)}\|_{0}\leq k$, for $r=1,\ldots,\nu$. We would argue that this
is a more natural assumption, since it is the change in the regression
coefficients, rather than the pre- and post-change coefficients themselves,
that is the quantity of interest in this statistical problem. Practically, the
assumption that all regression coefficients are sparse can be violated in
applications. For instance, Kraft and Hunter (2009) argued that in genetic
studies, “many, rather than few, variant risk alleles are responsible for the
majority of the inherited risk of each common disease”, leading to non-sparse
regression coefficients. However, in such examples, the task of detecting
sparse changes in these regression coefficients over time can still be of
interest in, e.g., identifying different development stages in gene regulatory
networks in species (Hatleberg and Hinman, 2021). Furthermore, our ‘sparsity-
in-change’ assumption is also more in line with the assumptions made in the
high-dimensional change-in-mean problem (see, e.g. Cho and Fryzlewicz, 2014;
Jirak, 2015; Wang and Samworth, 2018), where the pre- and post-change mean
vectors are regarded as nuisance parameters and sparsity assumptions only need
to be placed on vectors of changes for successful detection and localization
of the changepoints.

Allowing for dense pre- and post-change regression coefficients makes the
changepoint estimation problem considerably more challenging. In particular,
the general strategy employed by existing works on high-dimensional regression
changepoints that relies on forming good estimators of
$(\beta^{(r)}:r\in\\{1,\ldots,\nu\\})$ will unlikely be successful here. Our
first contribution in this paper is to propose a novel methodology, which we
call charcoal (changepoint in regression via a complementary-sketching
algorithm), and works by forming a projected response vector and a sequence of
projected design matrices to eliminate the dense nuisance parameter. For
simplicity of exposition, we consider the single changepoint scenario, where
$\beta_{i}=\beta^{(1)}\mathbbb{1}_{\\{i\leq
z_{1}\\}}+\beta^{(2)}\mathbbb{1}_{\\{i>z_{1}\\}}$. Under the hypothesis that
the true change takes place at time $t$, we have

$\begin{cases}Y_{(0,t]}=X_{(0,t]}\beta^{(1)}+\epsilon_{(0,t]},\\\
Y_{(t,n]}=X_{(t,n]}\beta^{(2)}+\epsilon_{(t,n]},\end{cases}$ (2)

where the subscript $(0,t]$ indicates the concatenation of relating quantities
on $(1,\dots,t)$ and $(t,n]$ that of $(t+1,\dots,n)$. We may think of (2) as a
two-sample problem with different regression coefficients before and after
$t$. We assume throughout the paper that $n>p$ for otherwise it is impossible
to estimate the change when both pre- and post-change parameters are dense
(see further discussion at the beginning of Section 3). By invoking the
complementary sketching method of Gao and Wang (2022), we can find matrices
$A_{(0,t]}\in\mathbb{R}^{t\times(n-p)}$ and
$A_{(t,n]}\in\mathbb{R}^{(n-t)\times(n-p)}$ such that
$(A_{(0,t]}^{\top},A_{(t,n]}^{\top})^{\top}$ has orthogonal columns spanning
the orthogonal complement of the range of $X$. By forming the projected design
matrix $W_{t}:=A_{(0,t]}^{\top}X_{(0,t]}-A_{(t,n]}^{\top}X_{(t,n]}$ and the
projected response $Z:=A_{(0,t]}^{\top}Y_{(0,t]}+A_{(t,n]}^{\top}Y_{(t,n]}$,
we can eliminate the possibly dense nuisance parameter
$\zeta:=(\beta^{(2)}+\beta^{(1)})/2$ and conduct tests on
$\theta^{(1)}=(\beta^{(2)}-\beta^{(1)})/2$ is zero against that it is non-zero
and sparse.

In light of the true changepoint at time $z_{1}$, the hypothesized model (2)
is only correctly specified when $t=z_{1}$. The further $t$ is away from
$z_{1}$, the less different the two samples $(X_{(0,t]},Y_{(0,t]})$ and
$(X_{(t,n]},Y_{(t,n]})$ are, since one of the samples will be further
contaminated by the data points assigned to the wrong segment by the
hypothesized changepoint. Intuitively, we would expect the aforementioned two-
sample test statistics to peak around $t=z_{1}$, which can thus be used to
estimate the location of the single changepoint. Unfortunately, while good for
testing, these two-sample test statistics have variances too large for
accurate changepoint localization. Nevertheless, the general idea of using
complementary sketching to eliminate nuisance parameters is valid. We
introduce in Section 2 several alternative statistics based on the sketched
design $W_{t}$ and response $Z$ that do lead to good changepoint estimation
performance. In particular, we will show in Section 3 that a variant of the
charcoal procedure achieves a rate of convergence of order
$\sqrt{k/(n\|\theta^{(1)}\|_{2}^{2})}$, up to logarithmic factors. In the
course of investigating the theoretical properties, we have developed new
results in understanding the asymptotic behaviour of the sketched design
matrices by generalizing existing matrix-variate Beta distribution to rank-
deficient cases (Lemma 9 and Corollary 10) and extended sub-Gaussian bounds of
Beta random variables to the matrix variate case (Lemma 12), both of which may
be of independent interest.

###  1.1 Outline of the paper

We present the methodology in detail in Section 2, including several
algorithms that all use the complementary sketching idea. Section 3 provides
theoretical performance guarantees to the slight variants of those proposed in
Section 2. In the first part of Section 4, we conduct numerical experiments on
the charcoal methodology over a comprehensive range of settings for both
single and multiple changepoint estimation tasks and compare our methods with
other changepoint localization methods in the high-dimensional linear
regression context. In the second part of Section 4, we study a real data
example to identify changes for each gene in terms of its interaction with
other genes in the gene regulatory network across various development stages
of T cells. Section 5 collects the proofs of the main results while we gather
the proofs of the ancillary results in Section 6.

###  1.2 Notation

For a positive integer $p$, $[p]=\\{1,\dots,p\\}$ consists of all positive
integers not exceeding $p$. For vector $v=(v_{1},\dots,v_{p})^{\top}$,
$\operatorname{diag}(v)$ is a $p\times p$ matrix such that
$(\operatorname{diag}(v))_{i,j}=\mathbbb{1}_{\\{i=j\\}}v_{i}$ for $i,j\in[p]$.
We follow the usual definitions of
$\|v\|_{0}=\sum_{i\in[p]}\mathbbb{1}_{\\{v_{i}\neq 0\\}}$,
$\|v\|_{2}=(\sum_{i\in[p]}v_{i}^{2})^{1/2}$, $\|v\|_{1}=\sum_{i\in[p]}|v_{i}|$
and $\|v\|_{\infty}=\max_{i\in[p]}|v_{i}|$.

Given a matrix $A\in\mathbb{R}^{n\times m}$, we make it a convention that
$A=(A_{i,j})_{i\in[n],j\in[m]}=[A_{1}\mid\dots\mid
A_{m}]=(a_{1},\dots,a_{n})^{\top}$, where $a_{i}$ is the transpose of the
$i$th row of $A$ and $A_{j}$ is the $j$th column of $A$. Given any set
$S\subseteq\mathbb{R}$, we write $A_{S}$ to be the submatrix of $A$ with row
indices in $S$. For instance, given positive integers $s,t$ such that $1\leq
s<t\leq n$, $A_{(s,t]}:=(a_{s+1},\dots,a_{t})^{\top}$. We define the usual
norms for $A$ as follows
$\|A\|_{\mathrm{op}}:=\sup_{v\in\mathbb{R}^{m}:\|v\|_{2}=1}\|Av\|_{2}$ and
$\|A\|_{\max}:=\max_{i\in[n],j\in[m]}|A_{i,j}|$. Assuming $n=m$ in $A$,
$\operatorname{diag}(A)$ is an $n\times n$ matrix such that
$(\operatorname{diag}(A))_{i,j}:=\mathbbb{1}_{i=j}A_{i,j}$ for $i\in[n]$ and
$\mathrm{tr}(A):=\sum_{i\in[n]}A_{i,i}$.

For $n\geq m$, $\mathbb{O}^{n\times m}:=\\{O\in\mathbb{R}^{n\times
m}:O^{\top}O=I_{m}\\}$. We define
$\mathcal{S}^{p-1}:=\\{v\in\mathbb{R}^{p}:\|v\|_{2}=1\\}$ and the $k$-sparse
unit ball as $B_{0}(k):=\\{v\in\mathbb{R}^{p}:\|v\|_{2}\leq 1,\|v\|_{0}\leq
k\\}$.

##  2 Methodology

In this section, we describe in detail our charcoal algorithm for identifying
the changepoints in the problem setup of (1).

###  2.1 Single changepoint estimation

We start by focusing on the setting of a single changepoint estimation, i.e.,
$\nu=1$, which captures the essence of the difficulty of this problem. For
simplicity, we denote $z:=z_{1}$ for the location of the only changepoint and
write $m:=n-p$. The main idea is to use data-driven projections to sketch the
design matrix and the response vector to eliminate the effect of the nuisance
parameters.

Recall the data generating model (1). At each time point $t\in[n-1]$, we
perform a two-sample test for the equality of regression coefficients before
and after $t$ using data points $(x_{i},y_{i})_{i=1}^{t}$ and
$(x_{i},y_{i})_{i=t+1}^{n}$ respectively. Motivated by Gao and Wang (2022),
this can be achieved by constructing a matrix $A\in\mathbb{O}^{n\times m}$
whose columns span the orthogonal complement of the column space of $X$. We
then define for any $t\in[n-1]$

$\displaystyle W_{t}$
$\displaystyle:=A_{(0,t]}^{\top}X_{(0,t]}-A_{(t,n]}^{\top}X_{(t,n]}=2A_{(0,t]}^{\top}X_{(0,t]}\in\mathbb{R}^{m\times
p}$ $\displaystyle Z$
$\displaystyle:=A_{(0,t]}^{\top}Y_{(0,t]}+A_{(t,n]}^{\top}Y_{(t,n]}=A^{\top}Y\in\mathbb{R}^{m}.$

We define $\theta=(\beta^{(1)}-\beta^{(2)})/2$,
$\zeta=(\beta^{(1)}+\beta^{(2)})/2$ and $\xi=A^{\top}\epsilon\sim
N_{m}(0,\sigma^{2}I_{m})$. By the model construction, we have

$\displaystyle Z$
$\displaystyle=A^{\top}_{(0,z]}Y_{(0,z]}+A^{\top}_{(z,n]}Y_{(z,n]}=A^{\top}_{(0,z]}(X_{(0,z]}\beta^{(1)}+\epsilon_{(0,z]})+A^{\top}_{(z,n]}(X_{(z,n]}\beta^{(2)}+\epsilon_{(z,n]})$
(3)
$\displaystyle=A_{(0,z]}^{\top}X_{(0,z]}(\theta+\zeta)-A_{(z,n]}^{\top}X_{(z,n]}(\theta-\zeta)+\xi=W_{z}\theta+\xi,$

whence we have eliminated the nuisance parameter $\zeta$, and obtain the
sketched data in the form of $(Z,(W_{t})_{t\in[n-1]})$. By (3) and the
sparsity assumption on $\theta$, $Z$ can be approximated by a sparse linear
combination of the columns of $W_{z}$. Therefore, the changepoint localization
problem is reduced to finding $t$ such that $W_{t}$ forms a ‘best’ sparse
linear approximation to $Z$.

As mentioned in the introduction, a naive way to achieve this would be based
on the two-sample test statistics introduced in Gao and Wang (2022).
Specifically, let $Q=(Q_{1},\ldots,Q_{n-1})^{\top}$ be defined such that

$Q_{t}:=\\{\operatorname{diag}(W_{t}^{\top}W_{t})\\}^{-1/2}W_{t}^{\top}Z.$

We view $Q_{t}$ as the vector of the correlations between columns of $W_{t}$
and $Z$, where we naturally seek to find the time point $t$ such that such
correlations are as large as possible. To take into account of possible
observational errors, we first remove small entries of $Q_{t}$ via an
entrywise hard-thresholding operation $\operatorname*{hard}(Q_{t},\lambda)$
for
$\operatorname*{hard}(v,\lambda):(v_{i})_{i=1}^{p}\mapsto(v_{i}\mathbbb{1}_{\\{|v_{i}|\geq\lambda\\}})_{i=1}^{p}$,
where the threshold level $\lambda$ is a tuning parameter. This allows us to
estimate the location of the changepoint via
$\hat{z}^{\operatorname*{hard}}:=\operatorname*{arg\,max}_{t\in[n-1]}\|\operatorname*{hard}(Q_{t},\lambda)\|_{2}$.
Note that $\|\operatorname*{hard}(Q_{t},\lambda)\|_{2}$ is the statistic from
Gao and Wang (2022) to test whether the two samples $(X_{(0,t]},Y_{(0,t]})$
and $(X_{(t,n]},Y_{(t,n]})$ have the same regression coefficient against the
alternative that there is a sparse difference. As is argued before, we expect
that two-sample testing statistics gives the strongest signal against the null
of no change at $t=z$ — the only point where the two-sample problem is
correctly specified.

However, the changepoint estimator $\hat{z}^{\operatorname*{hard}}$ is less
than ideal in practice, as the discontinuity of the hard-thresholding function
creates large variabilities in the test statistics. Moreover, the theoretical
guarantees given in Gao and Wang (2022) becomes increasingly inapplicable for
test statistics away from the true changepoint as one of the two samples
contains a mixture of data both before and after the change. Coupled with the
fact that $W_{t}$ has large variance when $t$ is close to the boundary,
$Q_{t}$ may have a large number of entries above the hard-thresholding level
$\lambda=2\sqrt{\log p}$ as recommended in Gao and Wang (2022). Empirically,
this is evidenced by high variance of the test statistics near the two
endpoints of the interval for changepoint detection, as shown in Figure 1.
Quite often, this boundary effect may overwhelm the main signal near the true
changepoint, leading to a spurious changepoint being estimated near the
boundary. One way to alleviate the instability problem of
$\hat{z}^{\operatorname*{hard}}$ is to replace the hard-thresholding in
$\hat{z}^{\operatorname*{hard}}$ by a soft-thresholding operation on each
entry of $Q$. The changepoint is then estimated by
$\hat{z}^{\mathrm{soft}}=\operatorname*{arg\,max}_{t\in[n-1]}\|\operatorname*{soft}(Q_{t},\lambda)\|_{2}$,
where
$\mathrm{soft}(v,\lambda):(v_{i})_{i=1}^{p}\mapsto(\mathrm{sign}{(v_{i})}\max(|v_{i}|-\lambda,0))_{i=1}^{p}$
with a tuning parameter $\lambda$. The continuity of the soft-thresholding
function reduces the variance in the test statistics, and in the ensuing
changepoint estimator. However, as also shown in Figure 1, the sequence of
test statistics $(\|\operatorname*{soft}(Q_{t},\lambda)\|_{2})_{t\in[n-1]}$
could still exhibit undesirably large, although less so than
$(\|\operatorname*{hard}(Q_{t},\lambda)\|_{2})_{t\in[n-1]}$, variations when
$t$ is close to the boundary.

Figure 1:  Visualization of different test statistics for changepoint
localization. We compare the hard-thresholded test statistics
$\|\operatorname*{hard}(Q_{t},\lambda)\|_{2}$ used to estimate
$\hat{z}^{\operatorname*{hard}}$ (denoted by hard), the soft-thresholded test
statistics $\|\mathrm{soft}(Q_{t},\lambda)\|_{2}$ used to estimate
$\hat{z}^{\mathrm{soft}}$ (denoted by soft) and the projected statistics
$|\hat{v}^{\top}Q_{t}|$ in Algorithm 1 (denoted by proj) over two random
realizations. Here, $n=600$, $p=200$, $\|\beta_{1}-\beta_{2}\|_{0}=10$,
$\|\beta_{1}-\beta_{2}\|_{2}=8$, and the true change takes place at $z=180$,
as indicated by the dashed lines. In both panels, we observe that both
$\|\operatorname*{hard}(Q_{t},\lambda)\|_{2}$ and
$\|\mathrm{soft}(Q_{t},\lambda)\|_{2}$ exhibit relatively strong boundary
effect.

To avoid such boundary effect, we propose instead to aggregate the test
statistics $(Q_{t})_{t}$ via a projection-based approach. The key insight here
is that, away from the boundary, the matrix $Q=(Q_{1},\ldots,Q_{n-1})$ can be
well-approximated by a rank-one matrix whose leading left singular vector is
proportional to $\theta$. Hence, by first estimating $\theta/\|\theta\|_{2}$
via the leading left singular vector of $\mathrm{soft}(Q,\lambda)$, we can
aggregate each vector of correlation $Q_{t}$ along the direction of $\hat{v}$
and estimate the changepoint by
$\hat{z}:=\operatorname*{arg\,max}_{t\in[\alpha
n,(1-\alpha)n]}|\hat{v}^{\top}Q_{t}|$. This approach is summarized in
Algorithm 1. We allow Algorithm 1 to output both the changepoint estimator
$\hat{z}$ and a test statistic $H_{\max}$, which can be used in our multiple
changepoint algorithm to determine if an estimated changepoint is spurious.

To compute the sequence $(Q_{t})_{t\in[n-1]}$ in Algorithm 1, observe that the
same $A$ and $Z$ can be used for all $t\in[n-1]$ and hence only need to be
computed once. It is worth noting that we exploit the structures of the
sketched designs $(W_{t})_{t}$ to greatly simplify their computations. Recall
that $a_{t}\in\mathbb{R}^{n-p}$ is the $t$th row vector of $A$, i.e.,
$A=(a_{1},\dots,a_{n})^{\top}$. $W_{t}$ are computed via the simple iterative
scheme $W_{0}=0$ and $W_{t}=W_{t-1}+2a_{t}x_{t}^{\top}$ for $t\in[n-1]$. As a
common measure, we introduce the burn-in parameter ${\alpha}$ so that we forgo
the possibilities of having changepoints in
$(0,{\alpha}n)\cup((1-{\alpha})n,n)$.

While the main focus of our current work is the changepoint localization
problem, we remark that Algorithm 1 can be easily adapted to test the
existence of a single changepoint in the sequence of regression coefficients.
Specifically, we can construct the test

$\psi_{\alpha,\lambda,T}:=\mathbbb{1}\Bigl{\\{}\max_{t\in[\alpha
n,(1-\alpha)n]}\|\operatorname*{soft}(Q_{t},\lambda)\|\geq T\Bigr{\\}},$ (4)

where $T$ is some appropriate threshold.

Input: $X\in\mathbb{R}^{n\times p},Y\in\mathbb{R}^{n}$ satisfying $n-p>0$, a
soft threshold level $\lambda\geq 0$, burn-in parameter ${\alpha}\geq 0$

1 Set $m\leftarrow n-p$

2 Form $A\in\mathbb{O}^{n\times m}$ with columns orthogonal to the column
space of $X$

3 Compute $Z\leftarrow A^{\top}Y$

4 Set $W_{0}=\mathbf{0}_{m\times p}$

5 for  _$1\leq t\leq n-1$_ do

6 Compute $W_{t}\leftarrow W_{t-1}+2a_{t}x_{t}^{\top}$

7 Compute
$Q_{t}=\\{\operatorname{diag}(W_{t}^{\top}W_{t})\\}^{-1/2}W_{t}^{\top}Z$

8

9 end for

10Form $Q:=(Q_{\lfloor\alpha
n\rfloor},\dots,Q_{\lceil(1-\alpha)n\rceil})^{\top}$

11 Compute $\hat{v}\leftarrow$ the leading left singular vector of
$\operatorname*{soft}(Q,\lambda)$

Output: $\hat{z}:=\operatorname*{arg\,max}_{\alpha n\leq
t\leq(1-\alpha)n}|\hat{v}^{\top}Q_{t}|$ and $H_{\max}:=\max_{\alpha n\leq
t\leq(1-\alpha)n}\|\operatorname*{soft}(Q_{t},\lambda)\|$.

Algorithm 1 Pseudocode for changepoint estimation

Finally, we mention that another natural approach to find the $W_{t}$ whose
columns form the best sparse linear approximation of $Z=W_{z}\theta+\xi$ is to
fit a sparse linear model by regressing $Z$ against $W_{t}$ and compare the
goodness-of-fit across $t$ via the Bayesian Information Criterion (BIC). We
choose the BIC for the model selection purpose, though it is conceivably
straightforward to apply any other model selection criteria. The pseudocode
for this procedure is given in Algorithm 2. Specifically, for appropriately
chosen $(\lambda_{t})_{t\in[n-1]}$, we compute first the Lasso solutions in
Step 2 and then the corresponding BICs in Step 2. In practice, the sequence of
regularizing parameters $(\lambda_{t})_{t}$ may be chosen via cross-validation
for each $t$.

Input: $X\in\mathbb{R}^{n\times p},Y\in\mathbb{R}^{n}$ satisfying $n>p$,
${\alpha}>0$ and a sequence $(\lambda_{t})_{t\in[n-1]}$

1 Follow Algorithm 1 until line 1

2 for  _$1\leq t\leq n-1$_ do

3 Compute $W_{t}\leftarrow W_{t-1}+2a_{t}x_{t}^{\top}$

4 Compute the Lasso estimator
$\hat{\theta}_{t}\leftarrow\operatorname*{arg\,min}_{v\in\mathbb{R}^{p}}\bigl{\\{}\frac{1}{2m}\|Z-W_{t}v\|_{2}^{2}+\lambda_{t}\|v\|_{1}\bigr{\\}}$

5 Compute
$H_{t}\leftarrow-(\|Z-W_{t}\hat{\theta}_{t}\|_{2}^{2}+\|\hat{\theta}_{t}\|_{0}\log
m)$

6

7 end for

Output: $\hat{z}:=\operatorname*{arg\,max}_{{\alpha}n\leq
t\leq(1-{\alpha})n}H_{t}$

Algorithm 2 Pseudocode for changepoint estimation with Lasso with BIC

Algorithm 1 has a computational complexity of $O(n^{2}p)$, with the most
computationally intensive step being its Step 1 to form the sketching matrix
$A$ (e.g. via a QR decomposition). For Algorithm 2, each Lasso step has a
computational cost of $O(k^{2}n)$ (Efron et al., 2004), leading to an overall
computational complexity of $O(n^{2}(p+k^{2}))$. It is remarkable that for
sparse signals ($k=O(\sqrt{p})$), the changepoint algorithms we proposed here
has essentially the same computational complexity as the complementary-
sketching-based two-sample test (Gao and Wang, 2022) for any hypothesized
changepoint location $t$.

###  2.2 Multiple changepoint estimation

The single changepoint estimation procedure described above can be combined
with a generic top-down multiple changepoint localization method, such as
binary segmentation (Vostrikova, 1981), wild binary segmentation (Fryzlewicz,
2014) and its variants (e.g., Baranowski et al., 2019; Kovács et al., 2020;
Fryzlewicz, 2020), to iteratively identify multiple changepoints. For
concreteness, we describe an approach combining Algorithm 1 with the
narrowest-over-threshold method of Baranowski et al. (2019). Algorithm 3 is a
slight generalisation of Baranowski et al. (2019, Algorithm 1). It takes as
input a single changepoint estimation procedure $\hat{z}$ and a testing
procedure $\psi$. When the data $D_{1},\ldots,D_{n}$ are the covariate-
response pair $(X_{i},Y_{i})_{i\in[n]}$, we may apply Algorithm 1 or 2 to
obtain $\hat{z}$ and define $\psi(X,Y):=\mathbbb{1}_{\\{H_{\max}>T\\}}$ for
some $T$ using the output $H_{\max}$ of Algorithm 1. However, note that both
Algorithms 1 and 2 require the number of observations to be larger than the
dimension for the complementary sketching to work. If this is not satisfied,
we simply define $\hat{z}(X,Y):=0$ and $\psi(X,Y):=0$.

Essentially, in Algorithm 3, we generate multiple intervals and run the single
changepoint algorithm on each interval to obtain candidate changepoint
estimates and test results. We choose the candidate changepoint associated
with the narrowest interval for which the test rejects the null, and add that
to the set of estimated changepoints. We then segment the data at this
estimated changepoint, and repeat the above process recursively on the data to
the left and right segments, using only intervals lying completely within each
segment. The process terminates when none of the tests reject the null.
Furthermore, for practical reasons, we recommend combining Algorithm 3 with
some second-stage refinements, for which we discuss in more details in Section
4.4.

Input: Data $D_{1},\ldots,D_{n}$, number of intervals $M$, burn-in parameter
$\varpi>0$, single changepoint estimation procedure $\hat{z}$ and a single
changepoint testing procedure $\psi$

1 Set $\hat{Z}\leftarrow\emptyset$ and generate a set of $M$ intervals
$\mathcal{M}:=\\{(s_{1},e_{1}],\ldots,(s_{M},e_{M}]\\}$ independently and
uniformly from $\\{(a,b]:0\leq a<b\leq n\\}$.

2 Run $\textnormal{{NOT}}(0,n)$ where NOT is defined below.

3 Let $\hat{\nu}\leftarrow|\hat{Z}|$ and sort elements of $\hat{Z}$ in
increasing order to yield $\hat{z}_{1}<\cdots<\hat{z}_{\hat{\nu}}$.

Output: $\hat{z}_{1},\ldots,\hat{z}_{\hat{\nu}}$

4 Function  _NOT(_$s$ , $e$_)_

5 Set
$\mathcal{R}^{(s,e]}\leftarrow\\{m:(s_{m},e_{m}]\subseteq(s,e],\,\psi(D_{(s_{m}+n\varpi,e_{m}-n\varpi]})=1\\}$

6 if  _$\mathcal{R}^{(s,e]}\neq\emptyset$_ then

7
$m_{0}\leftarrow\operatorname*{arg\,min}_{m\in\mathcal{R}^{(s,e]}}(e_{m}-s_{m})$

8 $b\leftarrow s_{m_{0}}+\hat{z}(D_{(s_{m},e_{m}]})$

9 $\hat{Z}\leftarrow\hat{Z}\cup\\{b\\}$

10 NOT(_$s$ , $b$_)

11 NOT(_$b$ , $e$_)

12

13 end if

14

15 end

Algorithm 3 Pseudocode for multiple changepoint estimation

##  3 Theoretical guarantees

In this section, we establish theoretical guarantees for the changepoint
procedures proposed in Section 2. We start by focusing on the single
changepoint estimation problem. For simplicity of analysis, we will assume
that the noise variance $\sigma^{2}$ is known in this section, which by scale
invariance can be further assumed to be equal to 1. We discuss practical
aspects of estimating $\sigma^{2}$ in Section 4.1. We first present two
conditions, which we will need to establish the results in this section.

######  Condition 1.

All entries of $X$ are independent standard normals.

######  Condition 2.

$n,z,p$ satisfy that $n>p$, $z/n\to\tau\in(0,1)$ and $(n-p)/n\to\eta\in(0,1)$
as $\min(z,n,p)\to\infty$.

The design Condition 1 requires that the rows of the design matrix
$X=(x_{1},\dots,x_{n})^{\top}$ follow the isotropic Gaussian distribution.
Condition 2 specifies the asymptotic regime we work in. Note that the
assumption $n>p$ is necessary, since otherwise, even if $z$ is known, it is
impossible to test if $\|\theta\|_{2}=0$ against a sparse alternative (see the
discussion of condition (C2) in Gao and Wang (2022)). The key ingredient of
our theoretical analysis is the following proposition, which shows that
$W_{t}^{\top}W_{z}$ is close to a multiple of identity in terms of their
actions on sparse vectors. We impose both conditions 1 and 2 only to enable
the application of the existing random matrix theory on the limiting spectral
measure of matrix-variate Beta distributions in the proof of Proposition 1. In
principle, even if the above conditions are violated, the theoretical results
in the rest of the section hold for any data $(X,Y)$ such that (5) is
satisfied. In particular, we remark that the empirical study in Section 4.3
has demonstrated that our methodology exhibits good finite-sample performance
even when the above two conditions do not hold.

######  Proposition 1.

Suppose that Conditions 1 and 2 are satisfied and define

$g(t;z):=\begin{cases}4t(n-z)(n-p)/n^{2}&\text{if }1\leq t\leq z,\\\
4z(n-t)(n-p)/n^{2}&\text{if }z<t\leq n-1.\end{cases}$

There exists a constant $C_{\tau,\eta}>0$, depending only on $\tau$ and $\eta$
such that with probability 1, for any fixed $v\in\mathcal{S}^{p-1}$ and
$\ell\in[p]$, we have for all but finitely many $p$’s that

$\sup_{t\in[n-1]}\sup_{u\in
B_{0}(\ell)}u^{\top}\bigl{\\{}W_{t}^{\top}W_{z}-g(t;z)I_{p}\bigr{\\}}v\leq
C_{\tau,\eta}\sqrt{\ell n\log p}.$ (5)

Note that we suppress the dependence on $n$ and $t$ in the notation of
$g(t;z)$. Taking $\ell=1$ in the above proposition, we would expect
$\operatorname{diag}(W_{t}^{\top}W_{t})$ to concentrate around
$g(t;t)I_{p}=4t(n-t)(n-p)n^{-2}I_{p}$ for each $t$. This would allow us to
approximate the test statistics
$Q_{t}=\\{\mathrm{diag}(W_{t}^{\top}W_{t})\\}^{-1/2}(W_{t}^{\top}W_{z}\theta+W_{t}^{\top}\xi)$
for a fixed $t$. However, due to a lack of non-asymptotic probabilistic bounds
in random matrix theory on the convergence of the spectral measures of matrix-
variate Beta random matrices, we are unable to establish the said convergence
of $\operatorname{diag}(W_{t}^{\top}W_{t})$ uniformly over $t\in[n-1]$. As
such, we instead show the theoretical results for a slightly modified variant
of Algorithm 1, where we replace the definition of $Q_{t}$ by

$Q_{t}:=\sqrt{\frac{n}{t(n-t)}}W_{t}^{\top}Z.$

We will henceforth refer to the above variant of Algorithm 1 as Algorithm 1′.
It is worth noting that the latter is merely a proof device, and in practice
we always recommend applying Algorithm 1. Empirically, the primed variant has
a slightly worse but comparable estimation accuracies than Algorithm 1, which
can be seen in Table 1 in Section 4.1.

With this alternative choice of $Q_{t}$, after removing the perturbation
introduced by observational errors $\xi$’s with an appropriate soft-
thresholding tuning parameter $\lambda$, we expect

$H_{t}:=\|\mathrm{soft}(Q_{t},\lambda)\|_{2}\approx
h_{t}:=\|\theta\|_{2}\gamma_{t},$

where $\gamma_{t}$ is defined by

$\gamma_{t}:=g(t;t)\sqrt{\frac{n}{t(n-t)}}=\begin{cases}\frac{4(n-p)}{n}\sqrt{\frac{t}{n(n-t)}}(n-z)&\text{if
}1\leq t\leq z,\\\\[6.0pt] \frac{4(n-p)}{n}\sqrt{\frac{n-t}{nt}}z&\text{if
}z<t\leq n-1.\end{cases}$ (6)

Interestingly, $h_{t}$ (or $\gamma_{t}$) is proportional to the CUSUM
statistic in the univariate change-in-mean problem, whence $h_{t}$ attains its
maximum at $t=z$ (cf. Wang and Samworth, 2018, Equation (10)). By exploiting
the above observation, we establish in Theorem 2 that the testing procedure
(4) is capable of determining whether a (single) changepoint is present in the
regression data, as mentioned in Section 2.

Recall $\zeta:=(\beta^{(1)}+\beta^{(2)})/2$ and
$\theta:=(\beta^{(1)}-\beta^{(2)})/2$, where we regard $\zeta$ as a possibly
dense nuisance parameter and wish to localize the changepoint only assuming
the difference parameter $\theta$ is sparse, i.e., $\|\theta\|_{0}\leq k$ for
some unknown but fixed $k$ typically much smaller than $p$.

######  Theorem 2.

Assume Conditions 1 and 2 and that data $(X,Y)$ are generated according to (1)
with $\nu=1$. Suppose that $\|\theta\|_{0}\leq k$ satisfies $(k\log p)/n\to 0$
and that $\min(\tau,1-\tau)\geq{\alpha}$ for some known ${\alpha}$. There
exists
$c_{\tau,\eta,\alpha},C_{\tau,\eta,\alpha},c^{\prime}_{\tau,\eta,\alpha}>0$,
depending only on $\tau,\eta,\alpha$, such that for
$\lambda=c_{\tau,\eta,\alpha}\log p$, $T=C_{\tau,\eta,\alpha}\sqrt{k}\log p$,
the following holds.

  1. 1.

If $\theta=0$, then
$\psi_{\alpha,\lambda,T}(X,Y)\xrightarrow{\mathrm{a.s.}}0$.

  2. 2.

If $\|\theta\|_{2}>c^{\prime}_{\tau,\eta,\alpha}\frac{\sqrt{k}\log
p}{\sqrt{n}}$, then
$\psi_{\alpha,\lambda,T}(X,Y)\xrightarrow{\mathrm{a.s.}}1$.

We now turn our attention to the estimation in Algorithm 1. The key of
understanding the performance of Algorithm 1′ lies in an analysis of the
estimated projection vector $\hat{v}$ in Step 1 of the algorithm. By
Proposition 1, we expect $Q=(Q_{1},\ldots,Q_{n-1})$ to be well-approximated by
the rank-one matrix $\theta\gamma^{\top}$, where
$\gamma=(\gamma_{t})_{t\in[n-1]}$ is defined in (6). Thus, the oracle
projection direction to aggregate $Q$ is along $\theta/\|\theta\|_{2}$. The
following proposition shows that the estimated projection direction $\hat{v}$
is well-aligned with this oracle direction.

######  Proposition 3.

Assume Conditions 1 and 2 and that data $(X,Y)$ are generated according to (1)
with $\nu=1$. Suppose that $k\leq p/2$ and that
$\min(\tau,1-\tau)\geq{\alpha}$ for some known ${\alpha}$. There exists
$c_{\tau,\eta,{\alpha}},C_{\tau,\eta,{\alpha}}>0$, depending only on
$\tau,\eta,{\alpha}$, such that if
$\lambda>c_{\tau,\eta,{\alpha}}\max(1,\|\theta\|_{2})\log p$, then the
projection direction estimator $\hat{v}$ in Algorithm 1′ satisfies with
probability 1 for all but finitely many $p$’s that

$\sin\angle(\hat{v},\theta)\leq
C_{\tau,\eta,{\alpha}}\frac{\lambda\sqrt{k}}{\sqrt{n}\|\theta\|_{2}}.$

Equipped with Proposition 3, we are now in a position to state the convergence
rate of the changepoint estimator from Algorithm 1′.

######  Theorem 4.

Assume Conditions 1 and 2 and that data $(X,Y)$ are generated according to (1)
with $\nu=1$. Suppose that $\|\theta\|_{2}\leq 1$, $k\leq p/2$ and that
$\min(\tau,1-\tau)\geq{\alpha}>0$ for some known ${\alpha}$. There exists
$c_{\tau,\eta,{\alpha}}>0$, depending only on $\tau,\eta,{\alpha}$, such that
if $\lambda>c_{\tau,\eta,{\alpha}}\log p$, then the output $\hat{z}$ of
Algorithm 1′ with input $(X,Y)$, $\lambda$ and ${\alpha}$ satisfies with
probability 1 for all but finitely many $p$’s that

$\frac{|\hat{z}-z|}{n}\lesssim_{\tau,\eta,{\alpha}}\frac{\lambda^{2}\sqrt{k}}{\sqrt{n}\|\theta\|_{2}^{2}}.$

Theorem 4 shows that with a tuning parameter choice of order $\log p$, and
when $\|\theta\|_{2}$ is bounded (which is the more difficult regime for
estimation), Algorithm 1′ produces a consistent changepoint estimator with a
rate of convergence of order $k^{1/2}n^{-1/2}\|\theta\|_{2}^{-2}$ up to
logarithmic factors. However, in light of the testing viewpoint of Theorem 2,
in which it is possible to test apart the null of no change against a sparse
alternative if $\sqrt{n}\|\theta\|_{2}/\sqrt{k}$, up to logarithmic factors,
is sufficiently large, the rate in the above theorem appears to have an extra
factor of $\|\theta\|_{2}^{-1}$. This additional factor is likely to arise
from the technical difficulty of controlling the weak, though complex,
dependence between the estimated projection direction $\hat{v}$ and the
sketched Gaussian noises $(W_{t}^{\top}\xi:t\in[n-1])$. Indeed, the following
theorem shows that if $\hat{v}$ is estimated from an independent sample, then
the estimator from Algorithm 1 has a rate of convergence that agrees with what
is prescribed in Theorem 2, up to logarithmic factors.

######  Theorem 5.

Assume the same conditions as in Theorem 4. Let $(\tilde{X},\tilde{Y})$ be an
independent copy of $(X,Y)$. Let $Q$ be the matrix constructed in Step 8 of
Algorithm 1′ with input $(X,Y)$, $\lambda$ and ${\alpha}$. Suppose $\hat{v}$
is computed in Step 9 of Algorithm 1′ with input $(\tilde{X},\tilde{Y})$,
$\lambda$ and ${\alpha}$. Then
$\hat{z}:=\operatorname*{arg\,max}_{t\in[n-1]}|\hat{v}^{\top}Q_{t}|$ satisfies
with probability 1 for all but finitely many $p$’s that

$\frac{|\hat{z}-z|}{n}\lesssim_{\tau,\eta,{\alpha}}\frac{\lambda\sqrt{k}\log
p}{\sqrt{n}\|\theta\|_{2}}.$

This additional independent sample $(\tilde{X},\tilde{Y})$ may be obtained in
reality via a sample-splitting scheme. For example, we may take all odd time
points to construct the $Q$ matrix, and then use the even time points to
estimate the projection direction $\hat{v}$. However, such sample splitting is
necessary only from a technical viewpoint, and the algorithm typically
performs better without sample splitting in practice.

We remark that the rate in the above Theorem 5 is slower compared to the usual
results from change-in-mean problems, where rates of order
$n^{-1}\|\theta\|_{2}^{-2}$ are achievable under appropriate conditions (see
e.g. Wang and Samworth, 2018; Verzelen et al., 2020). Our slower rate arises
essentially from the approximation step in Proposition 1, which would not be
needed in a change-in-mean problem. It remains to be seen if the estimation
rate can be improved via alternative and possibly more refined analysis
routes.

We now turn our attention to theoretical guarantees in the multiple
changepoint setting. The following theorem shows that provided that we have a
good single changepoint estimation and testing procedure in any changepoint
problem, combining the narrowest-over-threshold with the single change
procedures yields a multiple changepoint estimation procedure of similar
accuracy with theoretical guarantees.

######  Theorem 6.

Let $D_{1},\ldots,D_{n}$ be a data sequence with changepoints
$0=z_{0}<z_{1}<\cdots<z_{\nu}<z_{\nu+1}=n$. We assume that $z_{i}-z_{i-1}\geq
n\Delta_{\tau}$ for all $i\in[\nu+1]$. Let $\mathcal{M}$ be defined as in
Algorithm 3. Write
$\mathcal{I}_{0}:=\\{(s,e]\in\mathcal{M}:(s+n\varpi,e-n\varpi]\cap\\{z_{1},\ldots,z_{\nu}\\}=\emptyset\\}$
and for $i\in[\nu]$, define
$\mathcal{I}_{i}:=\\{(s,e]\in\mathcal{M}:s\in[z_{i}-n\Delta_{\tau}/2,z_{i}-n\Delta_{\tau}/3],e\in[z_{i}+n\Delta_{\tau}/3,z_{i}+n\Delta_{\tau}/2]\\}$
and
$\tilde{\mathcal{I}}_{i}:=\\{(s,e]\in\mathcal{M}:\min\\{z_{i}-s,e-z_{i}\\}\geq
n\Delta_{\tau}/6\text{ and }e-s\leq n\Delta_{\tau}\\}$. Let $\hat{z}$ and
$\psi$ be the single changepoint estimation and testing procedure used in
Algorithm 3. Define the events

$\displaystyle\Omega_{0}$
$\displaystyle:=\\{\forall\,i\in[\nu],\exists\,m\in[M],\text{ s.t.
}(s_{m},e_{m}]\in\mathcal{I}_{i}\\},$ $\displaystyle\Omega_{1}$
$\displaystyle:=\\{\text{$\psi(D_{(s+n\varpi,e-n\varpi]})=0$ for all
$(s,e]\in\mathcal{I}_{0}$}\\},$ $\displaystyle\Omega_{2}$
$\displaystyle:=\bigcap_{i\in[\nu]}\bigl{\\{}\text{$\psi(D_{(s+n\varpi,e-n\varpi]})=1$
for all $(s,e]\in\mathcal{I}_{i}$}\bigr{\\}}$

and for some $\phi_{1},\ldots,\phi_{\nu}>0$,

$\Omega_{3}:=\bigcap_{i\in[\nu]}\bigl{\\{}\text{$|\hat{z}(D_{(s,e]})-(z_{i}-s)|\leq
n\phi_{i}$ for all $(s,e]\in\tilde{\mathcal{I}}_{i}$}\bigr{\\}}.$

Let $\hat{z}_{1},\ldots,\hat{z}_{\hat{\nu}}$ be the output of Algorithm 3 with
inputs $D_{1},\ldots,D_{n}$, $M>0$, $\varpi=\Delta_{\tau}/6$, $\hat{z}$ and
$\psi$. Assume further $\phi:=\max_{i\in[\nu]}\phi_{i}<\varpi$. We have on
$\Omega_{0}\cap\Omega_{1}\cap\Omega_{2}\cap\Omega_{3}$ that

$\hat{\nu}=\nu\text{ and }n^{-1}|\hat{z}_{i}-z_{i}|\leq\phi_{i}\text{ for all
$i\in[\nu]$}.$

In particular, we have

$\mathbb{P}\bigl{(}\hat{\nu}=\nu\text{ and
}n^{-1}|\hat{z}_{i}-z_{i}|\leq\phi_{i}\text{ for all $i\in[\nu]$}\bigr{)}\geq
1-\mathbb{P}(\Omega_{1}^{c})-\mathbb{P}(\Omega_{2}^{c})-\mathbb{P}(\Omega_{3}^{c})-\nu
e^{-\Delta_{\tau}^{2}M/36}.$

Note that the theorem is valid for any generic multiple changepoint estimation
that combines valid single changepoint estimation and testing procedures and
the top-down narrowest-over-threshold multiple changepoint estimation
paradigm. Hence, it can be applied in contexts other than the linear
regression setting here. The statement of Theorem 6 is slightly stronger than
the usual results on narrowest-over-threshold procedures, where $\phi_{i}$ are
taken identical.

Applying the above theorem to our specific problem, we extend the single
changepoint estimation result in Theorem 5 to the multiple changepoint setting
and establish the estimation accuracy of Algorithm 3. We first give the
following condition, which is the equivalent of Condition 2 in the multiple
changepoint setting.

######  Condition 3.

$n,p$ satisfy $n>p$ and that $(n-p)/n\rightarrow\eta$ as
$\min(n,p)\rightarrow\infty$. Assume further the changepoints satisfy
$0=z_{0}<z_{1}<\dots<z_{\nu}<z_{\nu+1}=n$ and $z_{i}-z_{i-1}>n\Delta_{\tau}$
for all $i\in[\nu+1]$ and $z_{i}/n\rightarrow\tau_{i}$ for $0\leq i\leq\nu+1$.

Due to the asymptotic nature of our theoretical results in the above, we
consider a sequence of Algorithm 3. To facilitate proof, we study a specific
coupling of the random intervals $\\{(s_{m},e_{m}]:m\in[M]\\}$ generated
across this sequence as follows:

$\displaystyle(\tilde{s}_{m},\tilde{e}_{m}]\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\mathrm{Unif}\bigl{(}\\{(\tilde{a},\tilde{b}]:0\leq
a<b\leq 1\\}\bigr{)},\text{ }\forall\;m\in[M],$ (7)
$\displaystyle\text{$s_{m}=\lfloor n\tilde{s}_{m}\rfloor$, $e_{m}=\lceil
n\tilde{e}_{m}\rceil$ for $m\in[M]$ and $n\in\mathbb{N}$}.$

Note that the intervals generated by (7) have the same law as those generated
in Algorithm 3.

######  Corollary 7.

Let $X$ and $Y$ be generated by (1) and write $D_{i}:=(X_{i},Y_{i})$ for
$i\in[n]$. Assume Conditions 1 and 3 hold. There exist
$c,C,c^{\prime},C^{\prime}>0$, depending only on $\alpha,\Delta_{\tau},\eta$,
such that the following holds. For $\alpha<1/6$, $\lambda=c\log p$,
$T=C\sqrt{k}\log p$, let $\hat{z}=\hat{z}_{\alpha,\lambda}$ be the sample-
splitted version of the single changepoint estimator defined in Algorithm 1′
and $\psi=\psi_{\alpha,\lambda,T}$ be the testing procedure defined in (4). If
$\Delta_{\tau}>3(1-\eta)$, $\frac{c^{\prime}\sqrt{k}\log
p}{\sqrt{n}}\leq\|\theta^{(i)}\|_{2}\leq 1$ and $\|\theta^{(i)}\|_{0}\leq k$
such that $\frac{k\log p}{n}\rightarrow 0$, then the output
$\hat{z}_{1},\ldots,\hat{z}_{\hat{\nu}}$ of Algorithm 3 with intervals
$\\{(s_{m},e_{m}]:m\in[M]\\}$ generated according to (7), inputs
$(D_{i})_{i\in[n]}$, $M>0$, $\varpi=\tau/6$, $\hat{z}_{\alpha,\lambda}$ and
$\psi_{\alpha,\lambda,T}$ satisfies with probability $1-\nu
e^{-\Delta_{\tau}^{2}M/36}$ that for all but finitely many $p$’s,

$\hat{\nu}=\nu\text{ and
}\frac{|\hat{z}_{i}-z_{i}|}{n}\leq\frac{C^{\prime}\lambda\sqrt{k}\log
p}{\sqrt{n}\|\theta^{(i)}\|_{2}}\text{ for all $i\in[\nu]$}.$

##  4 Numerical study

The implementation of our single- and multiple-changepoint algorithms are both
available in our GitHub repository.333https://github.com/gaofengnan/charcoal

###  4.1 Tuning parameter choice and comparison of variants

Theoretical analysis in Section 3 have assumed that the noise variance
$\sigma^{2}$ is known. In practice, we may obtain an upward-biased estimator
$\tilde{\sigma}$ as the median absolute deviation of entries of the $Q$
matrix. We note

$\displaystyle\mathrm{Var}(Q_{t})$
$\displaystyle=\mathbb{E}[\mathrm{Var}\\{\tilde{W}_{t}^{\top}(W_{z}\theta+\xi)\mid
X\\}]+\mathrm{Var}[\mathbb{E}\\{\tilde{W}_{t}^{\top}(W_{z}\theta+\xi)\mid
X\\}]$
$\displaystyle=\sigma^{2}\mathbb{E}(\tilde{W}_{t}^{\top}\tilde{W}_{t})+\mathrm{Var}(\tilde{W}_{t}^{\top}W_{z}\theta).$

Since $\mathbb{E}(\tilde{W}_{t}^{\top}\tilde{W}_{t})$ has all diagonal entries
equal to 1, every entry of $Q$ has a marginal variance of at least
$\sigma^{2}$.

(a) Varying $\rho$

(b) Varying $k$

(c) Varying dimensions

(d) Varying change locations

Figure 2:  The effect of choosing different $\lambda$ in a series of studies.
Our recommendation of $\lambda:=0.5\log p$ is marked in dashed vertical line
in each panel. Unless specified otherwise in each one, the panels share the
parameters $n=1200,p=1000,\tau=0.3,k=3,\rho=2$.

Algorithm 1 requires a soft-thresholding tuning parameter $\lambda>0$ as an
input. The theoretical results in Section 3 suggests using
$\lambda=c\sigma\log p$ for some $c>0$. We investigate here the performance of
our algorithm at different soft-thresholding levels $\lambda$. Specifically,
we computed the logarithmic average loss $|\hat{z}-z|$ of Algorithm 1 over 100
Monte Carlo repetitions for parameter settings of $n\in\\{600,1200,2400\\}$,
$p\in\\{n/3,5n/6\\}$, $\tau\in\\{0.2,0.3,0.4,0.5\\}$,
$k\in\\{1,3,10,\lfloor\sqrt{p}\rfloor,\lfloor 0.1p\rfloor,\lfloor
0.2p\rfloor,p\\}$, $\rho:=\|\theta\|_{2}\in\\{1,2,4,8\\}$, $\sigma=1$ and
various choices of $c\in[0.1,5]$. In all our simulations here and below, we
sample the vector of change in the regression coefficients $\theta$ uniformly
from the set $\\{v:\|v\|_{0}=k,\|v\|_{2}=\rho\\}$, and generate dense pre-
change vector from $N_{p}(0,\max\\{1,\rho^{2}\\}I_{p})$. Figure 2 illustrates
part of the simulation results where we vary one aspect of the parameters at a
time. From the figure, we see that a choice of $c=0.5$ provides good
statistical performance across the parameter settings considered, and we will
henceforth adopt this choice of $\lambda=0.5\tilde{\sigma}\log(p)$ in our
subsequent numerical studies.

We now compare the statistical performance of various versions of
complementary-sketching-based approaches proposed in the paper, including
Algorithms 1 and 2 from Section 2 and the slight variant Algorithm 1′
mentioned in Section 3 to facilitate theoretical analysis. For a demonstrative
purpose, we have also included the naive hard- and soft-thresholded
changepoint estimators $\hat{z}^{\operatorname*{hard}}$ and
$\hat{z}^{\operatorname*{soft}}$ mentioned just above Algorithm 1. We use the
$\lambda$ choice suggested in the previous paragraph for
$\hat{z}^{\operatorname*{hard}}$ and $\hat{z}^{\operatorname*{soft}}$,
Algorithm 1 and its variant, and choose $\lambda_{t}$ in Algorithm 2 via a
five-fold cross-validation for each $t\in[n-1]$. Empirical observations
suggest that Algorithms 1, 1′ and 2 work well without any burn-in (i.e.
$\alpha=0$). However, both $\hat{z}^{\operatorname*{soft}}$ and
$\hat{z}^{\operatorname*{hard}}$ do suffer from more serious boundary effects,
as seen in the large root mean squared errors in Table 1. In addition,
Algorithm 1 has roughly the same but slightly better estimation accuracy
compared to its primed variant. This justifies our recommendation of Algorithm
1 over its primed variant, and the similarity in performance further
consolidates the relevance of our theoretical analysis on the primed variant
as a proof device.

$n$ | $p$ | $k$ | $\rho$ | $\hat{z}^{\operatorname*{soft}}$ | $\hat{z}^{\operatorname*{hard}}$ | Alg$1$ | Alg$1^{\prime}$ | Alg$2$
---|---|---|---|---|---|---|---|---
$600$ | $200$ | $3$ | $1$ | $50.1$ | $178.35$ | $12.6$ | $40.79$ | $22.34$
|  |  | $2$ | $3.98$ | $88.76$ | $3.1$ | $5.01$ | $4.58$
|  |  | $4$ | $25.37$ | $92.3$ | $2.14$ | $4.33$ | $2.29$
|  | $14$ | $1$ | $31.7$ | $122.11$ | $37.36$ | $83.59$ | $129.71$
|  |  | $2$ | $6.4$ | $88.04$ | $6.14$ | $7.32$ | $9.83$
|  |  | $4$ | $3.88$ | $58.92$ | $3.82$ | $4.67$ | $2.95$
$1200$ | $400$ | $3$ | $1$ | $11.24$ | $196.81$ | $10.95$ | $13.34$ | $19.54$
|  |  | $2$ | $3.15$ | $136.02$ | $3.69$ | $4.58$ | $4.96$
|  |  | $4$ | $35.65$ | $101.17$ | $1.67$ | $4.36$ | $2.09$
|  | $20$ | $1$ | $39.17$ | $146.98$ | $18.14$ | $37.63$ | $168.56$
|  |  | $2$ | $5.06$ | $144.5$ | $7.29$ | $7.59$ | $15.51$
|  |  | $4$ | $4.72$ | $106.84$ | $2.72$ | $2.93$ | $4.47$
Table 1: Comparisons of performances of $\hat{z}^{\operatorname*{soft}}$,
$\hat{z}^{\operatorname*{hard}}$, Algorithms 1, 1′ and 2 in various settings
in terms of root mean squared error. No burn-in is applied anywhere
($\alpha=0$). The true change takes place at $0.3n$.

###  4.2 Comparisons with other methods

From the discussion above, we recommend using Algorithms 1 and 2 for their
robustness against the choice of the burn-in parameter $\alpha$, and for the
more accurate estimation of Algorithm 1 over its primed variant. We will
henceforth focus on Algorithms 1 and 2, which we call
$\mathsf{charcoal}_{\mathsf{proj}}$ and $\mathsf{charcoal}_{\mathsf{lasso}}$,
respectively. In this section, we compare the performance of
$\mathsf{charcoal}_{\mathsf{proj}}$ and $\mathsf{charcoal}_{\mathsf{lasso}}$
with existing approaches in the literature. Specifically, we will compare
against the VPBS algorithm of Rinaldo et al. (2021), two-sided Lasso-based
approaches of Lee et al. (2016) (LSS) and Leonardi and Bühlmann (2016) (LB),
and a two-stage refinement approach of Kaul et al. (2019) (KJF). We have used
the authors’ own implementation for VPBS and KJF, and Kaul et al. (2019)’s
implementation for LSS. We have implemented LB ourselves using the recommended
tuning parameter choices as in Leonardi and Bühlmann (2016). It is worth
noting that none of the four existing methods in the literature were designed
to estimate changes in the regression coefficients when both the pre- and
post-change coefficients are dense.

We compare the performance in terms of mean absolute loss of various methods
in a single changepoint estimation task for $n\in\\{600,1200\\}$,
$p\in\\{n/3,5n/6\\}$, $\tau\in\\{0.1,0.3,0.5\\}$,
$k\in\\{3,\lfloor\sqrt{p}\rfloor,p\\}$, $\rho\in\\{1,2,4,8\\}$. Table 2 shows
a representative subset of these simulation results. We see that none of VPBS,
LB, KJF and LSS show any sign of consistent estimation as their average loss
do not decrease as the signal strength increases. On the other hand, both
$\mathsf{charcoal}_{\mathsf{proj}}$ and $\mathsf{charcoal}_{\mathsf{lasso}}$
have shown highly promising performance in various settings. It is surprising
that $\mathsf{charcoal}_{\mathsf{proj}}$ and
$\mathsf{charcoal}_{\mathsf{lasso}}$ also seem to work even when the vector of
change is dense. We notice that $\mathsf{charcoal}_{\mathsf{lasso}}$ shows
better estimation accuracy when either the signal strength $\rho$ is high or
the vector of change $\theta$ is dense.

$p$ | $k$ | $\rho$ | $\mathsf{charcoal}_{\mathsf{proj}}$ | $\mathsf{charcoal}_{\mathsf{lasso}}$ | VPBS | LB | KJF | LSS
---|---|---|---|---|---|---|---|---
$400$ | $3$ | $1$ | 7.2 | $13.2$ | $452.4$ | $556.1$ | $238.8$ | $472.2$
|  | $2$ | 2.2 | $3.5$ | $476.3$ | $569.2$ | $239.3$ | $364.1$
|  | $4$ | 1.1 | $1.5$ | $434.2$ | $532.8$ | $239.1$ | $272.1$
|  | $8$ | 0.7 | $0.8$ | $326.3$ | $496.8$ | $239.1$ | $310.8$
| $20$ | $1$ | 12.4 | $85.4$ | $422.7$ | $528.8$ | $238.9$ | $479.5$
|  | $2$ | 3.0 | $9.2$ | $494.9$ | $546.8$ | $238.9$ | $284.5$
|  | $4$ | 2.0 | $2.6$ | $431.9$ | $553.1$ | $239.1$ | $268.5$
|  | $8$ | $1.9$ | 0.8 | $356.2$ | $513.3$ | $239.3$ | $261.5$
| $400$ | $1$ | 162.2 | $344.2$ | $477.8$ | $569.8$ | $238.8$ | $429.9$
|  | $2$ | 46.3 | $338.4$ | $504.0$ | $583.2$ | $238.8$ | $252.4$
|  | $4$ | $25.3$ | 13.3 | $446.3$ | $554.1$ | $238.9$ | $285.6$
|  | $8$ | $20.7$ | 3.0 | $355.6$ | $487.6$ | $239.1$ | $250.1$
$1000$ | $3$ | $1$ | 60.7 | $113.3$ | $241.6$ | $429.5$ | $237.2$ | $227.3$
|  | $2$ | 8.3 | $11.8$ | $243.4$ | $441.4$ | $239.0$ | $228.2$
|  | $4$ | 2.9 | $4.0$ | $239.5$ | $366.9$ | $243.9$ | $230.6$
|  | $8$ | $2.4$ | 1.4 | $235.1$ | $245.1$ | $262.2$ | $230.7$
| $31$ | $1$ | $300.3$ | $364.9$ | $233.4$ | $440.1$ | $238.8$ | 227.4
|  | $2$ | 71.7 | $140.9$ | $242.5$ | $469.5$ | $238.9$ | $228.3$
|  | $4$ | $16.0$ | 12.5 | $251.3$ | $358.4$ | $238.9$ | $224.5$
|  | $8$ | $13.7$ | 4.6 | $244.5$ | $249.0$ | $238.2$ | $230.1$
| $1000$ | $1$ | $275.5$ | $359.8$ | $232.6$ | $483.0$ | $239.3$ | 231.8
|  | $2$ | $256.9$ | $320.8$ | $238.4$ | $447.4$ | $238.9$ | 229.2
|  | $4$ | $224.1$ | 91.0 | $242.7$ | $378.2$ | $239.1$ | $228.0$
|  | $8$ | $194.5$ | 39.6 | $246.4$ | $253.5$ | $242.4$ | $226.7$
Table 2: Average loss of various changepoint methods under different settings.
Other parameters: $n=1200$, $z=360$. The method with the least average loss in
each line is marked in bold.

###  4.3 Model misspecification

While we have focused on the Gaussian Orthogonal Ensemble (GOE) design (i.e.
$X$ has independent $N(0,1)$ entries) and Gaussian noise in the theoretical
analysis, our methodology can be applied in more general settings. In this
subsection, we investigate the robustness of the estimation accuracy of our
method to deviations from this Gaussian distributional assumptions.
Specifically for $n=1,200$, $p=400$, $\tau=0.3$, $k=20$ and
$\rho\in\\{1.5^{0},\dots,1.5^{8}\\}$, we varied the design matrix $X$ to have
either $N_{p}(0,\Sigma)$ rows, where $\Sigma=(0.7^{|i-j|})_{i,j\in[p]}$ has an
autoregressive Toeplitz structure, or independent Rademacher entries. We also
vary the noise distribution to take $t_{4}$, $t_{6}$, centred
$\mathrm{Exp}(1)$ or Rademacher distributions. Overall, we see from Figure 3
that the performance of $\mathsf{charcoal}_{\mathsf{lasso}}$ is robust to both
non-GOE design matrices and discrete, heavy-tailed or skewed noise
distributions. Similar results hold for the
$\mathsf{charcoal}_{\mathsf{proj}}$ method.

(a) Varying design matrices

(b) Varying noise distributions

Figure 3:  Robustness to varying design matrices and noise distributions.

###  4.4 Multiple changepoints

As mentioned in Section 2, our charcoal algorithms can be easily combined with
generic multiple changepoint methods to deal with multiple changepoints, and
we proposed a specific version in Algorithm 3 of such a multiple changepoint
localization procedure. We run Algorithm 3 with $\varpi=0$ and $M=200$. For
the single changepoint estimator input $\hat{z}$ in Algorithm 3, we employ
Algorithm 1 with the recommended value of $\lambda$ in Section 4.1 and the
burn-in parameter $\alpha=0.05$. For the testing procedure input $\psi$ of
Algorithm 3, we run Algorithm 1 to obtain output $H_{\max}$ and define
$\psi(X,Y)=\mathbbb{1}_{\\{H_{\max}>T\\}}$, where the testing threshold $T$ is
chosen by a Monte Carlo simulation as follows. We generate $B=1000$ pairs of
$(X,Y)$ under model (1) with $\nu=0$ (i.e. no changepoint), and run Algorithm
1 with these synthetic $(X,Y)$ pairs and the same $\lambda$ and $\alpha$
choices as above. This would return $B$ test statistics
$(H_{\max}^{b})_{b\in[B]}$, which are used to estimate an upper $0.01/M$
quantile by fitting a generalized extreme value distribution (Smith, 1985).

While Algorithm 3 already produces good estimators of multiple changepoints,
its performance can be further improved by the following post-processing
procedures. Such post-processing has previously been described in e.g.,
Fryzlewicz (2018, Section 2 of online supplement).

Specifically, after we obtain an initial candidate set of changepoints using
Algorithm 3, we iteratively run the test prescribed by (4) on the largest time
interval containing each candidate changepoint as the only estimated
changepoint, and remove that candidate changepoint if the test is non-
significant. For the remaining candidate changepoints
$\hat{z}_{1},\ldots,\hat{z}_{\hat{\nu}}$, we refine their estimated locations
in two steps. We first perform a ‘midpoint’ refinement, where we use (for
instance) Algorithm 2 to output a refined estimator $\tilde{z}_{i}$ based on
data
$\\{(x_{t},y_{t}):t\in((\hat{z}_{i-1}+\hat{z}_{i})/2,(\hat{z}_{i}+\hat{z}_{i+1})/2]\\}$
for each $i\in[\hat{\nu}]$. Here, we use the convention that $\hat{z}_{0}=0$
and $\hat{z}_{\hat{\nu}+1}=n$. Using midpoints between successive estimated
changepoints ensures that each
$((\hat{z}_{i-1}+\hat{z}_{i})/2,(\hat{z}_{i}+\hat{z}_{i+1})/2]$ contains with
high probability at most one true changepoint. However, it does not use the
full data available around each true changepoint. As such, we perform a second
refinement step after this, where we use Algorithm 2 to output a further
refined estimator $\hat{z}_{i}^{\mathrm{refined}}$ based on data
$\\{(x_{t},y_{t}):t\in(\tilde{z}_{i-1}+\alpha n,\tilde{z}_{i+1}-\alpha n]\\}$
for each $i\in[\hat{\nu}]$ with $\alpha$ being the burn-in parameter as in
Algorithm 3. Again, we use the convention that $\tilde{z}_{0}=0$ and
$\tilde{z}_{\hat{\nu}+1}=n$. For both refinement steps, we may also use
Algorithm 1 in place of Algorithm 2, and they have very similar performances
in our numerical experiment. For definiteness and simpler presentation, we
employ Algorithm 2 for both refinement steps in the following numerical
experiments.

$n$ | $p$ | $k$ | $\rho_{\min}$ |  $\hat{\nu}-\nu$ value counts | Haus | ARI
---|---|---|---|---|---|---
|  |  |  | $-3$ | $-2$ | $-1$ | $0$ | $1$ |  |
$1200$ | $200$ | $3$ | 0.8 | $0$ | $0$ | $96$ | $4$ | $0$ | $292.8$ | $0.742$
|  |  | 1.2 | $0$ | $0$ | $22$ | $78$ | $0$ | $75.4$ | $0.918$
|  |  | 1.6 | $0$ | $0$ | $0$ | $98$ | $2$ | $8.8$ | $0.978$
|  | $10$ | 0.8 | $0$ | $2$ | $97$ | $1$ | $0$ | $304.9$ | $0.710$
|  |  | 1.2 | $0$ | $0$ | $42$ | $55$ | $3$ | $141.1$ | $0.856$
|  |  | 1.6 | $0$ | $0$ | $1$ | $96$ | $3$ | $18.0$ | $0.960$
|  | $100$ | 0.8 | $3$ | $67$ | $30$ | $0$ | $0$ | $591.7$ | $0.303$
|  |  | 1.2 | $0$ | $4$ | $88$ | $8$ | $0$ | $319.3$ | $0.611$
|  |  | 1.6 | $0$ | $0$ | $52$ | $46$ | $2$ | $217.1$ | $0.759$
$2400$ | $400$ | $3$ | 0.8 | $0$ | $0$ | $25$ | $75$ | $0$ | $155.3$ | $0.881$
|  |  | 1.2 | $0$ | $0$ | $0$ | $100$ | $0$ | $14.3$ | $0.975$
|  |  | 1.6 | $0$ | $0$ | $0$ | $100$ | $0$ | $10.1$ | $0.983$
|  | $10$ | 0.8 | $0$ | $15$ | $53$ | $32$ | $0$ | $376.9$ | $0.720$
|  |  | 1.2 | $0$ | $0$ | $2$ | $98$ | $0$ | $37.3$ | $0.945$
|  |  | 1.6 | $0$ | $0$ | $1$ | $99$ | $0$ | $21.0$ | $0.970$
|  | $100$ | 0.8 | $42$ | $57$ | $1$ | $0$ | $0$ | $1154.9$ | $0.184$
|  |  | 1.2 | $0$ | $32$ | $54$ | $14$ | $0$ | $647.0$ | $0.457$
|  |  | 1.6 | $0$ | $0$ | $14$ | $84$ | $2$ | $376.9$ | $0.658$
Table 3: Summary of results of multiple changepoint estimations under (M1) and
(M2) described in Section 4.4 with $\rho\in\\{0.8,1.2,1.6\\}$ and
$k\in\\{3,10,100\\}$. The first nine rows of the table corresponds to setting
(M1) and the last nine rows corresponds to (M2). The columns of
$\hat{\nu}-\nu$ tabulates the difference in number of estimated and true
changepoints over 100 Monte Carlo repetitions. The ‘Haus’ and ‘ARI’ columns
measure the average Hausdorff distance and the average adjusted rand index
between the discovered partition and the true partition over 100 repetitions.

We assume that the regression noise level $\sigma$ is known and consider the
following two multiple changepoint specifications in our simulations: (M1)
$n=1200$, $p=200$, $\nu=3$ and three changepoints are located at
$z=(z_{1},z_{2},z_{3})=(240,540,900)$ with signal sizes
$(\|\theta^{(1)}\|_{2},\|\theta^{(2)}\|_{2},\|\theta^{(3)}\|_{2})=\rho_{\min}\times(1,1.5,2)$
and sparsity
$\|\theta^{(1)}\|_{0}=\|\theta^{(2)}\|_{0}=\|\theta^{(3)}\|_{0}=k$
respectively for various $\rho_{\min}$ and $k$; (M2) $n=2400$, $p=400$,
$\nu=4$ and four changepoints are located at
$z=(z_{1},z_{2},z_{3},z_{4})=(720,1320,1800,2160)$ with signal sizes
$(\|\theta^{(1)}\|_{2},\|\theta^{(2)}\|_{2},\|\theta^{(3)}\|_{2},\|\theta^{(4)}\|_{2})=\rho_{\min}\times(1,1.15,1.45,2.18)$
and sparsity
$\|\theta^{(1)}\|_{0}=\|\theta^{(2)}\|_{0}=\|\theta^{(3)}\|_{0}=\|\theta^{(4)}\|_{0}=k$
respectively for various $\rho_{\min}$ and $k$.

Note that for (M2), the signal sizes are chosen such that

$\|\theta^{(i)}\|_{2}^{2}\frac{(z_{i}-z_{i-1})(z_{i+1}-z_{i})(z_{i+1}-z_{i-1}-p)}{(z_{i+1}-z_{i-1})^{2}}$

is approximately constant for each $i\in[\nu]$, which according to Gao and
Wang (2022) means that the effective signal-to-noise ratio of testing for each
changepoint $z_{i}$ within the interval $(z_{i-1},z_{i+1}]$ is almost
constant. Table 3 reports the multiple changepoint estimation performances for
both (M1) and (M2) with $\rho_{\min}\in\\{0.8,1.2,1.6\\}$ and
$k\in\\{3,10,100\\}$. The multiple changepoint estimation accuracy is measured
in terms of the difference between the number of estimated and true
changepoints, the average Hausdorff distance between the sets
$\\{z_{i}:i\in[n]\\}$ and $\\{\hat{z}^{\mathrm{refined}}_{i}:i\in[\hat{n}]\\}$
and finally the average adjusted Rand index (ARI) (Rand, 1971) of the
estimated segments against the truth, over 100 Monte Carlo repetitions. We see
from Table 3 that the promising single changepoint estimation performance of
our methodology carries over to the multiple changepoint settings.

Figure 4 visualizes the simulation results by showing the histograms of
estimated changepoints in four of the parameter settings shown in Table 3. It
is worth noting that in the bottom two panels of the figure, where the
effective signal-to-noise ratios are chosen to be approximately constant for
all the four changepoints, we indeed see a similar number of times in
identifications of each changepoint.

(a) $n=1200$, $p=200$, $k=3$, $\rho_{\min}=0.8$

(b) $n=1200$, $p=200$, $k=3$, $\rho_{\min}=1.6$

(c) $n=2400$, $p=400$, $k=10$, $\rho_{\min}=0.8$

(d) $n=2400$, $p=400$, $k=10$, $\rho_{\min}=1.6$

Figure 4: Histogram of estimated changepoint locations in four settings. The
true changes take place at $z=(240,540,900)$ for the $(n,p)=(1200,200)$
specifications in the above two panels with the signal strengths at respective
changes being $\rho_{\min}\times(1,1.5,2)$. For the two panels below with
$(n,p)=(2400,400)$, $z=(720,1320,1800,2160)$ with
$\rho_{\min}\times(1,1.15,1.45,2.175)$. The locations of true changes are
marked in lightly-coloured dashed vertical lines in each plot.

###  4.5 Real data example

In this subsection, we showcase how the charcoal algorithm can be applied to a
single-cell gene expression dataset from Suo et al. (2023). The original data
consists of the logarithmic normalized gene expression levels of 3211 genes
measured in 11853 cells along the developmental trajectory from proliferating
double positive cells (DP(P) cells) to quiescent double positive cells (DP(Q)
cells), $\alpha\beta$T entry cells (ABT cells) and finally to CD4+ T cells.
These cells have been ordered in pseudotime according to their development
stage in Suo et al. (2023), which we use as our timeline (see Figure 5). We
are interested in understanding the change in the gene regulatory networks
along this time trajectory. We can estimate the changepoints by modelling the
logarithmic normalized expressions using Gaussian graphical models and seek
changes in the nodewise regression coefficients of each gene against the
remaining genes. To speed up the computation, we preprocess the data by
subsampling $1/3$ of the original cells and only using genes that have non-
zero expression in at least $5\%$ of the cells. Our preprocessed data is
available on the GitHub repository. The changepoints are estimated using
Algorithm 1 with tuning parameters chosen as suggested in Section 4.1. In
Table 4, we list the genes that reported most significant test statistics in
their nodewise regression coefficients along this pseudotime trajectory. From
Figure 5, we see that most of the changes are identified immediately before
the boundary between the DP(P) and DP(Q) boundary, and most of the associated
genes TK1, CKAP2L, TTK, ARHGEF39, DEPDC1, SPC25, GTSE1, HMMR, CENPA are well-
known regulators for cell proliferation in biology (Bitter et al., 2020; Mills
et al., 1992; Zhou et al., 2018; Zhang et al., 2019; Guo et al., 2016, see,
e.g.). The change in nodewise regression coefficient of the RAG2 gene occurred
immediately before the DP(Q) and ABT boundary, which agrees with the existing
literature that RAG2 is a regulator for T cell development (Kalman et al.,
2004).

Figure 5: Ordered pseudotime of cells in the real data example of Section 4.5. Each plot point represents a cell, coloured by its annotated cell type. Vertical lines corresponds to estimated changepoint locations of the most significant changes in the nodewise regression coefficients as described in Section 4.5. Gene | Changepoint | Top interacting partners
---|---|---
TK1 | 495 | PTP4A1, KIF20B, CENPF, SPRY1, ZWINT,
|  | CALCRL, RHPN1, LYPD3, COMMD3, LINC00672
CKAP2L | 430 | TOGARAM1, DEPDC1B, AP001816.1, TUBA1C,
|  | FANCI, IGF2BP2, REC8, TOB1, ZMAT3, STK11IP
RAG2 | 1178 | SMPD3, AL365440.2, LZTFL1, AEBP1, HIST1H2BJ,
|  | MTSS1, CD1C, CSNK2B, CASC15, SLC29A1
TTK | 407 | UBE2S, KNL1, CDC20, TRAV19, DDIT3, AC023157.3,
|  | AC012360.3, IL23A, DNTT, USP53
ARHGEF39 | 396 | HJURP, CD1A, SLC25A25, CCDC152, MBTD1, PON1,
|  | H1F0, RNF125, APH1B, DDX3Y
DEPDC1 | 444 | AL138899.1, USPL1, RIPK4, SERPINF1, EPHB6,
|  | MTSS1, SLC8A1-AS1, SLC5A3, HDAC4, SGK1
SPC25 | 442 | ATF3, ITGAE, CDC42EP3, AC136475.5, EPS8, NINJ2,
|  | NDC80, ZNF280D, L3MBTL3, FBLN5
GTSE1 | 437 | MID1IP1, HIST1H2AG, GADD45G, PSRC1, FBLN5,
|  | HIST1H2BN, CCDC171, ARMH1, SERTAD2, DDX3Y
HMMR | 516 | TAX1BP3, LAIR1, SERP2, LANCL2, MANEA-DT,
|  | TIMP1, CSRNP1, TSGA10, CKAP5, RGS16
CENPA | 503 | TRBV7-3, SOCS1, FRMD4B, CDKN1A, FXYD2,
|  | PTPN12, NLGN4X, NINL, KLRG1, SWT1
Table 4: List of genes with most significant changes in the nodewise
regression coefficients, together with their changepoint locations and top 10
interacting partners.

##  5 Proof of main results

###### Proof of Proposition 1.

Define $\kappa_{1}=\kappa_{1}(n,z,p):=z(n-z)(n-p)/n^{3}$, which under
Condition 2 is $O(1)$. We decompose

$\displaystyle W_{t}^{\top}W_{z}-\frac{4tn\kappa_{1}}{z}I_{p}$
$\displaystyle=W_{t}^{\top}W_{z}-\frac{t}{z}W_{z}^{\top}W_{z}+\frac{t}{z}W_{z}^{\top}W_{z}-\frac{4tn\kappa_{1}}{z}I_{p}$
$\displaystyle={4t}\biggl{\\{}\biggl{(}\frac{S_{0,t}}{t}-\frac{S_{0,z}}{z}\biggr{)}S_{0,n}^{-1}S_{z,n}+\frac{n}{z}\biggl{(}\frac{W_{z}^{\top}W_{z}}{4n}-\kappa_{1}I_{p}\biggr{)}\biggr{\\}}.$

Define
$\Omega_{\nu}:=\\{\|S_{0,n}^{-1}S_{z,n}\|_{\mathrm{op}}\|S_{0,z}/z\|_{\mathrm{op}}\leq\nu\\}$.
By Lemma 14, we have a $(1/2)$-net $\mathcal{N}$ of $B_{0}(\ell)$ with
cardinality at most $(5ep/\ell)^{\ell}$ such that we have

$\displaystyle\mathbb{P}($ $\displaystyle\sup_{u\in
B_{0}(\ell)}|u^{\top}(S_{0,t}/t-S_{0,z}/z)S_{0,n}^{-1}S_{z,n}v|\geq x)$
$\displaystyle\leq\mathbb{E}[\mathbb{P}(2\sup_{u\in\mathcal{N}}|u^{\top}(S_{0,t}/t-S_{0,z}/z)S_{0,n}^{-1}v|\geq
x\mid
S_{0,z},S_{z,n})\mathbbb{1}_{\Omega_{\nu}}]+\mathbb{P}(\Omega_{\nu}^{\mathrm{c}})$
$\displaystyle\leq
52\biggl{(}\frac{5ep}{\ell}\biggr{)}^{\ell}\exp\\{-t^{2}x^{2}/(32z\nu^{2}\\})+\mathbb{P}(\Omega_{\nu}^{\mathrm{c}}),$

where the first inequality holds by that $\Omega_{\nu}$ is measurable with
respect to the $\sigma$-algebra generated by $(S_{0,z},S_{z,n})$ and the
second by Lemma 12 and a union bound. Define $\lambda_{\min}(A)$ and
$\lambda_{\max}(A)$ for any generic symmetric matrix $A$ to be the smallest
and largest eigenvalues of $A$, respectively. By Wainwright (2019, Theorem
6.1), we have

$\displaystyle\mathbb{P}\biggl{(}\lambda_{\max}\biggl{(}\frac{S_{z,n}}{n-z}\biggr{)}\geq
2\bigl{(}1+\sqrt{p/(n-z)}\bigr{)}^{2}\biggr{)}\leq\exp\biggl{(}-\frac{(n-z)\bigl{(}1+\sqrt{p/(n-z)}\bigr{)}^{2}}{8}\biggr{)},$
$\displaystyle\mathbb{P}\biggl{(}\lambda_{\max}\biggl{(}\frac{S_{0,z}}{z}\biggr{)}\geq
2\bigl{(}1+\sqrt{p/z}\bigr{)}^{2}\biggr{)}\leq\exp\biggl{(}-\frac{z\bigl{(}1+\sqrt{p/z}\bigr{)}^{2}}{8}\biggr{)},$
$\displaystyle\mathbb{P}\biggl{(}\lambda_{\min}\biggl{(}\frac{S_{0,n}}{n}\biggl{)}\leq\frac{(1-\sqrt{{p}/{n}})^{2}}{4}\biggr{)}\leq\exp\biggl{(}-\frac{n(1-\sqrt{p/n})^{2}}{8}\biggr{)}.$

We define
$\nu:=16[(n-z)/n](1+\sqrt{{p}/{z}})^{2}(1+\sqrt{{p}/(n-z)})^{2}(1-\sqrt{{p}/{n}})^{-2}$.
By a union bound, we arrive at

$\mathbb{P}(\Omega_{\nu}^{\mathrm{c}})=\mathbb{P}(\|S_{0,n}^{-1}S_{z,n}\|_{\mathrm{op}}\|S_{0,z}/z\|_{\mathrm{op}}\leq\nu)\leq\exp\biggl{(}-\frac{n\bigl{(}1-\sqrt{p/n}\bigr{)}^{2}}{8}\biggr{)}+2\exp(-p/8).$

Combining the above displays and setting $x:=16\nu\sqrt{\ell
z\log(5ep/\ell)}/t$, we have by a union bound that

$\displaystyle\mathbb{P}\biggl{(}\sup_{u\in
B_{0}(\ell)}|u^{\top}(S_{0,t}/t-S_{0,z}/z)S_{0,n}^{-1}v|\geq x\biggr{)}\leq
p^{-7}+e^{-{n(1-\sqrt{p/n})^{2}}/{8}}+2e^{-p/8}.$

Taking another union bound over $t\in[z]$, and by the Borel–Cantelli lemma, we
have with probability 1, for all but finitely many $p$’s that

$\sup_{t\in[z]}\sup_{u\in
B_{0}(\ell)}u^{\top}\bigl{(}W_{t}^{\top}W_{z}-(t/z)W_{z}^{\top}W_{z}\bigr{)}v\leq
C^{\prime}_{\tau,\eta}\sqrt{\ell p\log(ep/\ell)},$ (8)

for some constant $C^{\prime}_{\tau,\eta}>0$ that depends only on $\tau$ and
$\eta$. By Gao and Wang (2022, (22)), we have for all but finitely many $p$’s
that

$\sup_{u\in
B_{0}(\ell)}u^{\top}\biggl{(}\frac{W_{z}^{\top}W_{z}}{4n}-\kappa_{1}I_{p}\biggr{)}u\leq(4+o(1))\sqrt{\frac{(\ell+4)\log(10ep/\ell)}{n}}\bigl{\\{}(\kappa_{1}+\kappa_{2})\sqrt{n/p}+\kappa_{1}\bigr{\\}},$

where $\kappa_{2}>0$ is again a constant depending only on $\tau$ and $\eta$.
This, together with the first claim of Lemma 8, implies that with probability
1, for all but finitely many $p$’s, we have

$\sup_{t\in[z]}\sup_{u\in
B_{0}(\ell)}\frac{4tn}{z}u^{\top}\biggl{(}\frac{W_{z}^{\top}W_{z}}{4n}-\kappa_{1}I_{p}\biggr{)}v\leq
C^{\prime\prime}_{\tau,\eta}\sqrt{\ell p\log(ep/\ell)}.$ (9)

The conclusion follows by combining (8) and (9), and the corresponding
inequality for $t\in[n-1]\setminus[z]$. ∎

###### Proof of Theorem 2.

Applying Proposition 1, we have with probability 1 that for all but finitely
many $p$’s that

$\sup_{t\in[z]}\sup_{u\in
B_{0}(\ell)}u^{\top}\biggl{(}W_{t}^{\top}W_{z}-\frac{4t(n-z)(n-p)}{n^{2}}I_{p}\biggr{)}\frac{\theta}{\|\theta\|_{2}}\lesssim_{\tau,\eta}\sqrt{p\ell\log
p}.$ (10)

Let $S:=\mathrm{supp}(\theta)$. Taking $\ell=k$ in (10), we have

$\displaystyle\sup_{{\alpha}n\leq t\leq z}$
$\displaystyle\,\biggl{\|}\sqrt{\frac{n}{t(n-t)}}\biggl{(}\biggl{(}W_{t}^{\top}W_{z}-\frac{4t(n-z)(n-p)}{n^{2}}I_{p}\biggr{)}\theta\biggr{)}_{S}\biggr{\|}_{2}$
$\displaystyle\leq\sup_{{\alpha}n\leq t\leq z}\sup_{u\in
B_{0}(k)}u^{\top}\biggl{(}\sqrt{\frac{n}{t(n-t)}}W_{t}^{\top}W_{z}-\frac{4t^{1/2}(n-z)(n-p)}{n^{3/2}(n-t)^{1/2}}I_{p}\biggr{)}\theta$
$\displaystyle\lesssim_{\tau,\eta,{\alpha}}\sqrt{k\log p}\|\theta\|_{2}.$ (11)

By (10) and the second claim of Lemma 8, with probability 1 that for all but
finitely many $p$’s,

$\displaystyle\sup_{t\in[z]}\|(W_{t}^{\top}W_{z}\theta)_{S^{c}}\|_{2}$
$\displaystyle=\sup_{t\in[z]}\,\biggl{\|}\biggl{(}\biggl{\\{}W_{t}^{\top}W_{z}-\frac{4t(n-z)(n-p)}{n^{2}}I_{p}\biggr{\\}}\theta\biggr{)}_{S^{c}}\biggr{\|}_{2}$
$\displaystyle\lesssim_{\tau,\eta}p\sqrt{\log p}\|\theta\|_{2}.$

For any $Q\in\mathbb{O}^{p\times p}$, we have
$X\stackrel{{\scriptstyle\mathrm{d}}}{{=}}XQ^{\top}=:\tilde{X}$, and the
latter has the corresponding sketching matrix $\tilde{A}=AQ^{\top}$ because
$\tilde{A}^{\top}\tilde{A}=I_{n-p}$ and
$\tilde{A}^{\top}\tilde{X}=QA^{\top}XQ^{\top}=\mathbf{0}_{(n-p)\times p}$. As
such, for any $Q\in\mathbb{O}^{p\times p}$ such that $Q\theta=\theta$,
$Q(W_{t}^{\top}W_{z}\theta)=(QW_{t}Q^{\top})^{\top}(QW_{z}Q^{\top})\theta\stackrel{{\scriptstyle\mathrm{d}}}{{=}}W_{t}^{\top}W_{z}\theta$.
In particular, $(W_{t}^{\top}W_{z}\theta)_{S^{c}}$ is spherically symmetric on
$\mathbb{R}^{p-k}$. Hence, by Lemma 13 (with a choice of $\delta=2p^{-4}$),
with probability 1 we have for all but finitely many $p$’s that

$\displaystyle\sup_{{\alpha}n\leq t\leq
z}\,\biggl{\|}\sqrt{\frac{n}{t(n-t)}}(W_{t}^{\top}W_{z}\theta)_{S^{c}}\biggr{\|}_{\infty}$
$\displaystyle\leq\sup_{{\alpha}n\leq t\leq
z}\|(W_{t}^{\top}W_{z}\theta)_{S^{c}}\|_{2}\sqrt{\frac{n}{t(n-t)}}\frac{4\sqrt{\log
p}}{\sqrt{p-k}}$ $\displaystyle\lesssim_{\tau,\eta,{\alpha}}\|\theta\|_{2}\log
p$ (12)

Let $X=QT$ be the QR decomposition of $X$ and define
$B_{t}:=Q_{(0,t]}^{\top}Q_{(0,t]}$. By Equation (16) in the proof of Gao and
Wang (2022, Proposition 8), there exists $C_{\eta}>0$, depending only on
$\eta$ that for any fixed $p$ and $t$ and $j\in[p]$, with probability
$1-p^{-4}$, we have

$\frac{1}{n}(W_{t}^{\top}W_{t})_{j,j}\leq\frac{4\mathrm{tr}(B_{t}(I_{p}-B_{t}))}{p}+C_{\eta}\sqrt{\frac{\log
p}{n}}\leq 1+C_{\eta}\sqrt{\frac{\log p}{n}},$

where the final inequality follows from the fact that
$\|B_{t}\|_{\mathrm{op}}\leq 1$. Taking union bounds over $j\in[p]$ and
$t\in[z]$, and applying the Borel–Cantelli lemma, we have with probability 1
that for all but finitely many $p$’s,

$\sup_{\alpha n\leq t\leq
z}\sup_{j\in[p]}\,\frac{n}{t(n-t)}(W_{t}^{\top}W_{t})_{j,j}\lesssim_{\tau,{\alpha}}1.$

Furthermore, applying the Gaussian tail bound followed by a union bound, we
have with probability 1 for all but finitely many $p$’s that

$\sup_{\alpha n\leq t\leq
z}\biggl{\|}\sqrt{\frac{n}{t(n-t)}}W_{t}^{\top}\xi\biggr{\|}_{\infty}\leq
4\sqrt{\log p}\sup_{\alpha n\leq t\leq
z}\sup_{j\in[p]}\,\frac{n}{t(n-t)}(W_{t}^{\top}W_{t})_{j,j}\lesssim_{\tau,{\alpha}}\sqrt{\log
p}.$ (13)

We now work on the probability 1 event $\Omega$, such that (11), (12), (13)
all hold for all but finitely many $p$’s.

For sufficiently large $c_{\tau,\eta,{\alpha}}$, we have the right-hand side
of (12) and (13) are both dominated by $\lambda/2$. Hence, on $\Omega$, we
have for all $t\in[\alpha n,z]$ that

$\displaystyle H_{t}$
$\displaystyle=\biggl{\|}\operatorname*{soft}\biggl{(}\sqrt{\frac{n}{t(n-t)}}W_{t}^{\top}(W_{z}\theta+\xi),\lambda\biggr{)}\biggr{\|}_{2}$
$\displaystyle=\biggl{\|}\operatorname*{soft}\biggl{(}\sqrt{\frac{n}{t(n-t)}}(W_{t}^{\top}W_{z})_{S,S}\theta_{S}+\sqrt{\frac{n}{t(n-t)}}(W_{t}^{\top}\xi)_{S},\lambda\biggr{)}\biggr{\|}_{2}$

Writing
$\tilde{H}_{t}:=\bigl{\|}\sqrt{\frac{n}{t(n-t)}}(W_{t}^{\top}W_{z})_{S,S}\theta_{S}\bigr{\|}_{2}$,
we have by the triangle inequality and (13) that on $\Omega$,

$\sup_{t\in[\alpha
n,z]}\bigl{|}H_{t}-\tilde{H}_{t}\bigr{|}\lesssim_{\tau,{\alpha}}\sqrt{k}\lambda+\sqrt{k\log
p}.$ (14)

Recall the definition of $\gamma_{t}$ in (6) and write

$h_{t}:=\gamma_{t}\|\theta\|_{2}=\frac{4t^{1/2}(n-z)(n-p)}{n^{3/2}(n-t)^{1/2}}\|\theta\|_{2}=\frac{4(n-p)}{n}\sqrt{\frac{t}{n(n-t)}}(n-z)\|\theta\|_{2}.$

Then by (11), we also have on $\Omega$ that

$\sup_{t\in[{\alpha}n,z]}|\tilde{H}_{t}-h_{t}|\lesssim_{\tau,\eta,{\alpha}}\sqrt{k\log
p}\|\theta\|_{2}.$ (15)

By a symmetric argument, both (14) and (15) hold for $t\in[z,(1-\alpha)n]$ and
consequently for $t\in[\alpha n,(1-\alpha)n]$ (with perhaps a slightly
different constant). Consequently, we have that with probability 1 for all but
finitely many $p$,

$\sup_{t\in[\alpha n,(1-\alpha)n]}|H_{t}-h_{t}|\leq
C_{1}\bigl{(}\sqrt{k}\lambda+\sqrt{k\log p}\|\theta\|_{2}\bigr{)},$

where $C_{1}$ depends only on $\tau$, $\eta$ and $\alpha$.

If $\theta=0$, then $h_{t}=0$ for all $t$, and hence for sufficiently large
$C_{\tau,\eta,\alpha}$, we have with probability $1$ for all but finitely many
$p$’s that $\sup_{t\in[\alpha n,(1-\alpha)n]}|H_{t}|\leq T$ and thus the first
conclusion holds.

For the second conclusion, we have for some $C_{2}$ depending only on
$\tau,\eta$ and $\alpha$ that

$\displaystyle|H_{z}|\geq|h_{z}|-\sup_{t\in[\alpha
n,(1-\alpha)n]}|H_{t}-h_{t}|$ $\displaystyle\geq
C_{2}\sqrt{n}\|\theta\|_{2}-C_{1}(\sqrt{k}\lambda+\sqrt{k\log
p}\|\theta\|_{2})$ $\displaystyle\geq
C_{2}(1-o(1))\sqrt{n}\|\theta\|_{2}-C_{1}c_{\tau,\eta,\alpha}\sqrt{k}\log p.$

The signal size condition on $\|\theta\|_{2}$ then ensures that
$|H_{z}|\geq(C_{2}c^{\prime}_{\tau,\eta,\alpha}/2-C_{1}c_{\tau,\eta,\alpha})\sqrt{k}\log
p$. Hence, for sufficiently large $c^{\prime}_{\tau,\eta,\alpha}$, we can
ensure that with probability 1 for all but finitely many $p$’s, we have
$\max_{t\in[\alpha n,(1-\alpha)n]}|H_{t}|\geq|H_{z}|\geq T$, completing the
proof. ∎

###### Proof of Proposition 3.

Define for each $t\in[\alpha n,(1-\alpha)n]$ a vector
$\tilde{Q}_{t}\in\mathbb{R}^{p}$ such that
$(\tilde{Q}_{t})_{S^{\mathrm{c}}}:=0$ and

$(\tilde{Q}_{t})_{S}:=\sqrt{\frac{n}{t(n-t)}}(W_{t}^{\top}W_{z})_{S,S}\theta_{S}.$

By (12) and (13) and their symmetric results for $t\in[z,(1-\alpha)n]$, we
have with probability 1 for all but finitely many $p$’s that

$\sup_{\alpha n\leq
t\leq(1-\alpha)n}\|Q_{t}-\tilde{Q}_{t}\|_{\infty}\lesssim_{\tau,\eta,{\alpha}}\max(1,\|\theta\|_{2})\log
p$ (16)

Recall the definition of $\gamma=(\gamma_{t})_{t\in[n-1]}$ in (6). Applying
(10) and its symmetric result for $t\geq z$ with $\ell=1$, we have with
probability $1$ for all but finitely many $p$’s that

$\sup_{t\in[n-1]}\|\tilde{Q}_{t}-\theta\gamma_{t}\|_{\infty}\lesssim_{\tau,\eta,{\alpha}}\|\theta\|_{2}\sqrt{\log
p}.$ (17)

By symmetry, both (16) and (17) are still valid when we replace the supremum
over $t\in[z,(1-{\alpha})n]$ instead. Define
$\tilde{Q}:=(\tilde{Q}_{\lceil\alpha
n\rceil},\dots,\tilde{Q}_{\lfloor(1-\alpha)n\rfloor})^{\top}$ and
$\gamma:=(\gamma_{t})_{t\in[\alpha n,(1-\alpha)n]}$. We then have

$\|Q-\theta\gamma^{\top}\|_{\infty}\leq\|Q-\tilde{Q}\|_{\infty}+\|\tilde{Q}-\theta\gamma^{\top}\|_{\infty}\lesssim_{\tau,\eta,{\alpha}}\max(1,\|\theta\|_{2})\log
p.$ (18)

By Wang and Samworth (2018, Propositions 2 and 4 in the online supplement),
for $c_{\tau,\eta,{\alpha}}$ large enough such that
$\lambda\geq\|Q-\theta\gamma^{\top}\|_{\infty}$, we have

$\sin\angle(\hat{v},\theta)\lesssim_{\tau,\eta,{\alpha}}\frac{\lambda\sqrt{kn}}{\|\theta\|_{2}\|\gamma\|_{2}},$

whence the desired form follows by noting
$\|\gamma\|_{2}\asymp_{\tau,\eta,{\alpha}}n$ by Condition 2. ∎

###### Proof of Theorem 4.

To simplify exposition, all statements should be interpreted as valid with
probability 1 for all but finitely many $p$’s. Write
$v:=\theta/\|\theta\|_{2}$ for simplicity, and note that
$\|v\|_{0}=\|\theta\|_{0}\leq k$. Since the estimator $\hat{z}$ is unchanged
if we replace $\hat{v}$ by $-\hat{v}$ in Algorithm 1, we may assume without
loss of generality that $\rho:=\hat{v}^{\top}v\geq 0$. Our strategy is to view
$(\hat{v}^{\top}Q_{t}:t\in[\alpha n,(1-\alpha)n])$ as a perturbation of a
multiple of $(\gamma_{t}:t\in[\alpha n,(1-\alpha)n])$, which is maximized at
$z$. By (18), we may choose $c_{\tau,\eta,\alpha}$ large enough such that
$\|Q-\theta\gamma^{\top}\|_{\infty}\leq\lambda$. By Proposition 3, we then
have

$\displaystyle\max_{t\in[\alpha
n,(1-\alpha)n]}\bigl{|}\hat{v}^{\top}(Q_{t}-\theta\gamma_{t})\bigr{|}$
$\displaystyle\leq\|\hat{v}\|_{1}\|Q-\theta\gamma^{\top}\|_{\max}\leq(\|v\|_{1}+\|\hat{v}-v\|_{1})\lambda$
$\displaystyle\lesssim_{\tau,\eta,{\alpha}}(\sqrt{k}+\sqrt{p}\|\hat{v}-v\|_{2})\lambda.$

From Proposition 3, there exists $C_{\tau,\eta,{\alpha}}^{\prime}>0$,
depending only on $\tau,\eta,{\alpha}$, such that

$\|\hat{v}-v\|_{2}\leq 2\sin\angle(\hat{v},\theta)\leq
C_{\tau,\eta,{\alpha}}^{\prime}\frac{\lambda\sqrt{k}}{\sqrt{n}\|\theta\|_{2}},$
(19)

which implies that

$\max_{t\in[\alpha
n,(1-\alpha)n]}\bigl{|}\hat{v}^{\top}(Q_{t}-\theta\gamma_{t})\bigr{|}\lesssim_{\tau,\eta,\alpha}\frac{\lambda^{2}\sqrt{k}}{\|\theta\|_{2}}.$
(20)

We may further assume that

$\frac{\lambda\sqrt{k}}{\sqrt{n}\|\theta\|_{2}}\leq\frac{\lambda^{2}\sqrt{k}}{\sqrt{n}\|\theta\|_{2}^{2}}\leq\frac{1}{C^{\prime}_{\tau,\eta,{\alpha}}}$
(21)

for all $p$’s, since for $p$ where this is not satisfied the result is
trivially true. Then, $\sin\angle(\hat{v},\theta)\leq 1/2$ and thus
$\rho=\\{1-\sin^{2}\angle(\hat{v},\theta)\\}^{1/2}\geq 1/2$. Consequently,
from (20) and (21), increasing $C^{\prime}_{\tau,\eta,\alpha}$ if necessary,
we have

$\hat{v}^{\top}\theta\gamma_{z}=\rho\|\theta\|_{2}\frac{4z^{1/2}(n-z)^{1/2}(n-p)}{n^{3/2}}\geq
2\max_{t\in[\alpha
n,(1-\alpha)n]}\bigl{|}\hat{v}^{\top}(Q_{t}-\theta\gamma_{t})\bigr{|},$

which implies in particular that $\hat{v}^{\top}Q_{\hat{z}}>0$. Now, since
$z=\operatorname*{arg\,max}_{t\in[n-1]}\gamma_{t}$ and
$\hat{z}=\operatorname*{arg\,max}_{t\in[\alpha
n,(1-\alpha)n]}\hat{v}^{\top}Q_{t}$, we have from (20) that

$\hat{v}^{\top}\theta\gamma_{z}-\hat{v}^{\top}\theta\gamma_{\hat{z}}\leq\hat{v}^{\top}Q_{z}-\hat{v}^{\top}Q_{\hat{z}}+2\max_{t\in[\alpha
n,(1-\alpha)n]}\bigl{|}\hat{v}^{\top}(Q_{t}-\theta\gamma_{t})\bigr{|}\lesssim_{\tau,\eta,{\alpha}}\frac{\lambda^{2}\sqrt{k}}{\|\theta\|_{2}}.$
(22)

On the other hand, by Wang and Samworth (2018, Lemma 7), we have

$\inf_{t\in[z-\min\\{z,n-z\\}/2,z+\min\\{z,n-z\\}/2]}\frac{\hat{v}^{\top}\theta\gamma_{z}-\hat{v}^{\top}\theta\gamma_{t}}{|z-t|}\gtrsim_{\tau,\eta}\|\theta\|_{2}\sqrt{n}.$
(23)

We arrive at the conclusion by combining (22) and (23). ∎

###### Proof of Theorem 5.

As in the proof of Theorem 4, all statements are valid with probability 1 for
all but finitely many $p$’s, and we may assume without loss of generality that
$\hat{v}^{\top}v\geq 0$. Let $\tilde{Q}_{t}$ be as in the proof of Proposition
3. The main difference to the proof of Theorem 4 will be an improvement of
(20) using the independence between $\hat{v}$ and $Q_{t}-\theta\gamma_{t}$.
Specifically, since

$\displaystyle\\{t(n-t)/n\\}^{1/2}\hat{v}^{\top}(Q_{t}-\tilde{Q}_{t})$
$\displaystyle=\hat{v}^{\top}W_{t}^{\top}(W_{z}\theta+\xi)-\hat{v}_{S}^{\top}(W_{t}^{\top}W_{z})_{S,S}\theta_{S}$
$\displaystyle=\hat{v}_{S^{c}}^{\top}(W_{t}^{\top}W_{z}\theta)_{S^{c}}+\hat{v}^{\top}W_{t}^{\top}\xi,$

we have that

$\displaystyle\max_{t\in[\alpha
n,(1-\alpha)n]}\bigl{|}\hat{v}^{\top}(Q_{t}-\theta\gamma_{t})\bigr{|}$
$\displaystyle\leq\max_{t\in[\alpha
n,(1-\alpha)n]}\biggl{|}\hat{v}^{\top}(\tilde{Q}_{t}-\theta\gamma_{t})+\sqrt{\frac{n}{t(n-t)}}\hat{v}^{\top}(W_{t}^{\top}W_{z}\theta)_{S^{c}}\biggr{|}$
$\displaystyle\qquad+\max_{t\in[\alpha
n,(1-\alpha)n]}\biggl{|}\sqrt{\frac{n}{t(n-t)}}\hat{v}^{\top}W_{t}^{\top}\xi\biggr{|}$
(24)

We control the two terms on the right-hand side of (24) separately. By (17),
(12), the Cauchy–Schwarz inequality and finally (19), the first term from the
above display (24) is bounded by

$\displaystyle\|\hat{v}\|_{1}$
$\displaystyle\biggl{(}\|\tilde{Q}-\theta\gamma^{\top}\|_{\max}+\max_{t\in[\alpha
n,(1-\alpha)n]}\biggl{\|}\sqrt{\frac{n}{t(n-t)}}\|\hat{v}\|_{2}\|(W_{t}^{\top}W_{z}\theta)_{S^{c}}\|_{2}\biggr{\|}_{\infty}\biggr{)}$
$\displaystyle\lesssim_{\tau,\eta,{\alpha}}(\sqrt{k}\|v\|_{2}+\sqrt{p}\|\hat{v}-v\|_{2})\|\theta\|_{2}\log
p\lesssim_{\tau,\eta,{\alpha}}\lambda\sqrt{k}\log p.$ (25)

On the other hand, since $\hat{v}$, $W_{t}$ and $\xi$ are mutually
independent, we have $\hat{v}^{\top}W_{t}^{\top}\xi\mid(\hat{v},W_{t})\sim
N(0,\|W_{t}\hat{v}^{\top}\|_{2}^{2})$. By Gao and Wang (2022, (22)), we have

$\sup_{t\in[\alpha
n,(1-\alpha)n]}\biggl{\|}\frac{1}{n}W_{t}^{\top}W_{t}\biggr{\|}_{\mathrm{op}}\lesssim_{\tau,\eta,{\alpha}}1.$

Hence, by Gaussian tail bounds followed by a union bound, we have

$\max_{t\in[\alpha
n,(1-\alpha)n]}\biggl{|}\sqrt{\frac{n}{t(n-t)}}\hat{v}^{\top}W_{t}^{\top}\xi\biggr{|}\lesssim_{\tau,\eta,{\alpha}}\sqrt{\log
p}.$ (26)

Substituting (25) and (26) into (24), we have

$\max_{t\in[\alpha
n,(1-\alpha)n]}\bigl{|}\hat{v}^{\top}(Q_{t}-\theta\gamma_{t})\bigr{|}\lesssim_{\tau,\eta,{\alpha}}\lambda\sqrt{k}\log
p.$ (27)

Following the same argument as in the proof of Theorem 4, with (27) replacing
(20), we arrive at the following counterpart to (22):

$\hat{v}^{\top}\theta\gamma_{z}-\hat{v}^{\top}\theta\gamma_{\hat{z}}\leq\hat{v}^{\top}Q_{z}-\hat{v}^{\top}Q_{\hat{z}}+2\max_{t\in[\alpha
n,(1-\alpha)n]}\bigl{|}\hat{v}^{\top}(Q_{t}-\theta\gamma_{t})\bigr{|}\lesssim_{\tau,\eta,{\alpha}}\lambda\sqrt{k}\log
p.$ (28)

Combining (28) with (23), the proof is complete. ∎

###### Proof of Theorem 6.

First observe that since

$\mathbb{P}(\Omega_{0}^{c})\leq\sum_{i=1}^{\nu}\prod_{m=1}^{M}\Bigl{(}1-\mathbb{P}\bigl{(}(s_{m},e_{m}]\in\mathcal{I}_{i}\bigr{)}\Bigr{)}\leq\nu\bigl{(}1-\Delta_{\tau}^{2}/36\bigr{)}^{M}\leq\nu\exp(-\Delta_{\tau}^{2}M/36),$

the second conclusion follows immediately from the first one. To establish the
first conclusion, we henceforth work on the event
$\Omega_{0}\cap\Omega_{1}\cap\Omega_{2}\cap\Omega_{3}$.

For $0\leq s<e\leq n$, we define the following sets

$\displaystyle\mathcal{M}^{(s,e]}$ $\displaystyle:=\\{m\in[M]:s\leq
s_{m}<e_{m}\leq e\\},$ $\displaystyle\mathcal{R}^{(s,e]}$
$\displaystyle:=\\{m\in\mathcal{M}^{(s,e]}:\psi(D_{(s_{m}+n\varpi,e_{m}-n\varpi]})=1\\},$
$\displaystyle\mathcal{Z}^{(s,e]}$
$\displaystyle:=\\{i\in[\nu]:z_{i}\in(s,e]\\},$
$\displaystyle\mathcal{Z}_{\mathrm{good}}^{(s,e]}$
$\displaystyle:=\\{i\in[\nu]:z_{i}\in(s,e],\min\\{z_{i}-s,e-z_{i}\\}\geq
n\Delta_{\tau}/2\\},$ $\displaystyle\mathcal{Z}_{\mathrm{bad}}^{(s,e]}$
$\displaystyle:=\\{i\in[\nu]:z_{i}\in(s,e],\min\\{z_{i}-s,e-z_{i}\\}<n\phi_{i}\\}.$

Note that on the event $\Omega_{0}$, we can associate each true changepoint
$z_{i}$ with an $m_{i}\in[M]$ such that
$(s_{m_{i}},e_{m_{i}}]\in\mathcal{I}_{i}$. On $\Omega_{2}$, we have

$\bigl{\\{}m_{i}:i\in\mathcal{Z}_{\mathrm{good}}^{(s,e]}\bigr{\\}}\subseteq\mathcal{R}^{(s,e]}.$
(29)

Recall the assumption $\phi<\varpi$. For any $(s_{0},e_{0}]\subset(s,e]$ such
that
$(s_{0},e_{0}]\cap\\{z_{i}:i\in[\nu],z_{i}\in(s,e]\\}\subseteq\mathcal{Z}_{\mathrm{bad}}^{(s,e]}$,
we have $(s_{0},e_{0}]\in\mathcal{I}_{0}$ and hence
$\psi(D_{(s_{0}+n\varpi,e_{0}-n\varpi]})=0$ on $\Omega_{1}$.

For any set of changepoints $\hat{Z}$, we can partition the original timeline
$(0,n]$ into $|\hat{Z}|+1$ segments, which we call segments induced by
$\hat{Z}$. We now prove by induction that as we update $\hat{Z}$ throughout
the recursion of Algorithm 3, for any $(s,e]$ induced by $\hat{Z}$, we have
$\mathcal{Z}^{(s,e]}=\mathcal{Z}^{(s,e]}_{\mathrm{good}}\cup\mathcal{Z}^{(s,e]}_{\mathrm{bad}}$.
The base case is trivially true as at the beginning of the algorithm,
$\hat{Z}=\emptyset$, so the only segment induced is $(0,n]$ so
$\mathcal{Z}^{(s,e]}=\mathcal{Z}^{(s,e]}_{\mathrm{good}}$ and
$\mathcal{Z}^{(s,e]}_{\mathrm{bad}}=\emptyset$ by our assumption that
$z_{i}-z_{i-1}\geq n\Delta_{\tau}$ for all $i\in[\nu+1]$.

Now assuming that the inductive hypothesis holds at some iteration of the
recursion in Algorithm 3. We show that the claimed statement still holds if a
new changepoint is estimated. Let $\hat{Z}$ be the set of changepoints
identified before this new changepoint, say $\hat{z}_{*}$, is added to it. We
necessarily have $\hat{z}_{*}=\hat{z}(D_{(s^{\prime},e^{\prime}]})+s^{\prime}$
for some $(s^{\prime},e^{\prime}]\in\mathcal{Z}^{(s,e]}$ where $(s,e]$ is
induced by $\hat{Z}$. From the inductive hypothesis, we have
$\mathcal{Z}^{(s,e]}=\mathcal{Z}^{(s,e]}_{\mathrm{good}}\cup\mathcal{Z}^{(s,e]}_{\mathrm{bad}}$
and $\mathcal{Z}^{(s,e]}_{\mathrm{good}}$ is necessarily non-empty, for
otherwise all changepoints in $(s,e]$ are within a distance of $n\phi$ to the
boundary of the interval, which together with the fact that $\phi<\varpi$,
implies that $\mathcal{M}^{(s,e]}\cap\mathcal{R}^{(s,e]}=\emptyset$,
contradicting the fact that a new changepoint is identified. Thus, there
exists some $i^{\prime}\in\mathcal{Z}^{(s,e]}_{\mathrm{good}}$, which by (29)
means that $m_{i^{\prime}}\in\mathcal{R}^{(s,e]}$. By the definition of
$m_{0}$ in Line 3 of Algorithm 3, we have $e_{m_{0}}-s_{m_{0}}\leq
e_{m_{i^{\prime}}}-s_{m_{i^{\prime}}}\leq n\Delta_{\tau}$. Thus, from the
condition of the theorem, we have that $(s_{m_{0}},e_{m_{0}}]$ contains at
most one changepoint. If
$(s_{m_{0}},e_{m_{0}}]\cap\\{z_{i}:i\in\mathcal{Z}^{(s,e]}\\}=\emptyset$, then
on $\Omega_{1}$, $\psi(D_{(s_{m_{0}}+n\varpi,e_{m_{0}}-n\varpi]})=0$,
contradicting $m_{0}\in\mathcal{R}^{(s,e]}$. If $(s_{m_{0}},e_{m_{0}}]$
contains a single changepoint $z_{i}$ for
$i\in\mathcal{Z}^{(s,e]}_{\mathrm{bad}}$, then since $\phi<\varpi$, we again
have on $\Omega_{1}$ that $\psi(D_{(s_{m_{0}}+n\varpi,e_{m_{0}}-n\varpi]})=0$,
a contradiction. By the inductive hypothesis, this implies that
$(s_{m_{0}},e_{m_{0}}]$ contains exactly one true change-point $z_{i_{0}}$ for
some $i_{0}\in\mathcal{Z}^{(s,e]}_{\mathrm{good}}$ and that
$\min\\{e_{m_{0}}-z_{i_{0}},z_{i_{0}}-s_{m_{0}}\\}\geq n\varpi$. Hence,
$(s_{m_{0}},e_{m_{0}}]\in\tilde{\mathcal{I}}_{i_{0}}$, and thus on
$\Omega_{3}$, we have $|\hat{z}_{*}-z_{i_{0}}|\leq n\phi_{i_{0}}$.

We finally check that the two new segments induced by
$\hat{Z}\cup\\{\hat{z}_{*}\\}$, say $(\hat{z}_{\mathrm{left}},\hat{z}_{*}]$
and $(\hat{z}_{*},\hat{z}_{\mathrm{right}}]$ for
$\hat{z}_{\mathrm{left}}<\hat{z}_{*}<\hat{z}_{\mathrm{right}}$, still satisfy
the inductive hypothesis. By symmetry, we may assume without loss of
generality that $\hat{z}_{*}\leq z_{i_{0}}$. Since
$|z_{i_{0}}-\hat{z}_{*}|\leq n\phi$, we have
$i_{0}\in\mathcal{Z}^{(\hat{z}_{\mathrm{left}},\hat{z}_{*}]}_{\mathrm{bad}}$.
For any $i\in\mathcal{Z}^{(\hat{z}_{\mathrm{left}},\hat{z}_{*}]}$ such that
$i<i_{0}$, we have $\hat{z}_{*}-z_{i}\geq z_{i_{0}}-z_{i}\geq n\Delta_{\tau}$,
and thus
$z_{i}\in\mathcal{Z}^{(\hat{z}_{\mathrm{left}},\hat{z}_{*}]}_{\mathrm{good}}\cup\mathcal{Z}^{(\hat{z}_{\mathrm{left}},\hat{z}_{*}]}_{\mathrm{bad}}$
by the inductive hypothesis, consequently,
$\mathcal{Z}^{(\hat{z}_{\mathrm{left}},\hat{z}_{*}]}=\mathcal{Z}^{(\hat{z}_{\mathrm{left}},\hat{z}_{*}]}_{\mathrm{good}}\cup\mathcal{Z}^{(\hat{z}_{\mathrm{left}},\hat{z}_{*}]}_{\mathrm{bad}}$.
Similarly, for $i\in\mathcal{Z}^{(\hat{z}_{*},\hat{z}_{\mathrm{right}}]}$, we
have $i>i_{0}$ and $z_{i}-\hat{z}_{*}\geq
z_{i}-z_{i_{0}}-(\hat{z}_{*}-z_{i_{0}})\geq n\Delta_{\tau}-\phi\geq
n\Delta_{\tau}/2$. Again by the inductive hypothesis,
$\mathcal{Z}^{(\hat{z}_{*},\hat{z}_{\mathrm{right}}]}=\mathcal{Z}^{(\hat{z}_{*},\hat{z}_{\mathrm{right}}]}_{\mathrm{good}}\cup\mathcal{Z}^{(\hat{z}_{*},\hat{z}_{\mathrm{right}}]}_{\mathrm{bad}}$.
This completes the induction.

As a consequence of the above inductive argument, we have shown that a new
changepoint will be identified in Algorithm 3 if and only if
$(s,e]\cap\\{z_{i}:i\in\mathcal{Z}^{(s,e]}_{\mathrm{good}}\\}\neq\emptyset$.
Thus, from the inductive claim, at the end of the recursion, each changepoint,
say $z_{i}$, must be less than $n\phi_{i}$ away from one of the end points of
the segments induced by $\hat{Z}$. This, as well as the assumption that
$z_{i}-z_{i-1}\geq n\Delta_{\tau}$ for all $i\in[\nu+1]$, means that
$|\hat{Z}|=\nu$ and that $|\hat{z}_{i}-z_{i}|\leq n\phi_{i}$ as desired. ∎

###### Proof of Corollary 7.

First we write $\phi_{i}=\frac{C^{\prime}\lambda\sqrt{k}\log
p}{\sqrt{n}\|\theta^{(i)}\|_{2}}$ and define $\Omega_{0}$, $\Omega_{1}$,
$\Omega_{2}$ and $\Omega_{3}$ as in Theorem 6. Observe that $\Omega_{1}$,
$\Omega_{2}$ and $\Omega_{3}$ has implicit dependence on $p$, whereas in the
specific coupling (7) considered in this theorem, $\Omega_{0}$ does not vary
with $p$. We have from the proof of Theorem 6 that $\mathbb{P}(\Omega_{0})\geq
1-\nu e^{-\Delta_{\tau}^{2}M/36}$. Hence, it suffices show that on
$\Omega_{0}$, we have for all but finitely many $p$’s that $\Omega_{1}$,
$\Omega_{2}$ and $\Omega_{3}$ hold simultaneously. We keep in mind that $M$ is
fixed and finite, and for the rest of the proof, we condition on a realization
of $(\tilde{s}_{m},\tilde{e}_{m})_{m=1}^{M}$ as in (7) such that $\Omega_{0}$
holds.

Let $\mathcal{I}_{0}$, $\mathcal{I}_{i}$ and $\tilde{\mathcal{I}}_{i}$ be
defined as in Theorem 6. We first establish $\Omega_{1}$ and $\Omega_{2}$. For
any interval $(s,e]$ with $e-s\leq p$, $\psi=0$ by definition. For every
$(s,e]\in\cup_{0\leq i\leq\nu}\,\mathcal{I}_{i}$ whose length is longer than
$p$, the fixed-ratio regime Condition 2 is true by the generating mechanism of
the intervals in (7), and it is straightforward to verify that Theorem 2
applies. As a result, there exist $c,c^{\prime},C$, which may depend on
$(s,e]$, such that the conclusion of Theorem 2 holds for each
$(s,e]\in\cup_{0\leq i\leq\nu}\,\mathcal{I}_{i}$ with $e-s>p$. Inspecting the
proof of Theorem 2 shows that we can take the maximum of all such
$c,c^{\prime},C$’s so that the conclusion of Theorem 2 holds for all intervals
in $\cup_{0\leq i\leq\nu}\,\mathcal{I}_{i}$ with length longer than $p$. As
such, we have, for all but finitely many $p$’s, $\Omega_{1}$ and $\Omega_{2}$
holds.

Now we turn to $\Omega_{3}$. Again, by reasoning similar to the above, we see
the conditions of Theorem 5 hold for each $(s,e]\in\tilde{\mathcal{I}}_{i}$,
and for the above-mentioned specific choices of $c$, which may depend on
$(s,e]$, the conclusion of Theorem 5 holds for each
$(s,e]\in\tilde{\mathcal{I}}_{i}$. We can again take the maximum of all such
$c$’s so that for all intervals in $\tilde{\mathcal{I}}_{i}$ for all
$i\in[\nu]$, the conclusion of Theorem 5 holds, i.e.,

$\frac{\hat{z}(D_{(s,e]})-(z_{i}-s)}{n}\leq
C^{\prime}_{(s,e]}\frac{\lambda\sqrt{k}\log p}{\sqrt{n}\|\theta\|_{2}}.$

Setting, e.g.,
$C^{\prime}=\max_{i\in[\nu]}\max_{(s,e]\in\tilde{\mathcal{I}}_{i}}C^{\prime}_{(s,e]}$,
we have, for all but finitely many $p$’s, $\Omega_{3}$ holds. Invoking Theorem
6 completes the proof. ∎

##  6 Ancillary results

We collect here the ancillary results and their proofs.

######  Lemma 8.

Fix $A\in\mathbb{R}^{p\times p}$ and $k\in[p]$. The following are true.

  1. 1.

If $A$ is symmetric, then $\sup_{u,v\in B_{0}(k)}u^{\top}Av\leq\sup_{v\in
B_{0}(2k)}v^{\top}Av$.

  2. 2.

$\sup_{v\in B_{0}(k)}\|Av\|_{2}\leq\sqrt{p/k}\sup_{u,w\in
B_{0}(k)}u^{\top}Aw$.

###### Proof.

For the first claim, fix $u,v\in B_{0}(k)$ and let $S$ and $T$ be their
respective support. Then by the symmetry of $A$, we have

$u^{\top}Av\leq\|A_{S,T}\|_{\mathrm{op}}\leq\|A_{S\cup T,S\cup
T}\|_{\mathrm{op}}=\sup_{w\in\mathcal{S}^{2k-1}}w^{\top}A_{S\cup T,S\cup
T}w\leq\sup_{v\in B_{0}(2k)}v^{\top}Av.$

The first claim then follows by taking supremum on the left-hand side.

For the second claim, define $\psi:=\sup_{u,w\in B_{0}(k)}u^{\top}Aw$. Write
$\binom{[p]}{k}:=\\{S\subseteq B:|S|=k\\}$. For any $v\in B_{0}(k)$, let
$\tilde{v}:=Av/\|Av\|_{2}$ and $T:=\mathrm{supp}(v)$. Then by the
Cauchy–Schwarz inequality, we have

$\displaystyle\|Av\|_{2}$
$\displaystyle=\tilde{v}^{\top}Av=\frac{1}{\binom{p-1}{k-1}}\sum_{S\in\binom{[p]}{k}}\tilde{v}_{S}^{\top}A_{S,T}v_{T}\leq\frac{1}{\binom{p-1}{k-1}}\sum_{S\in\binom{[p]}{k}}\|\tilde{v}_{S}\|_{2}\psi$
$\displaystyle\leq\biggl{\\{}\frac{1}{\binom{p-1}{k-1}}\sum_{S\in\binom{[p]}{k}}\|\tilde{v}_{S}\|_{2}^{2}\biggr{\\}}^{1/2}\biggl{\\{}\frac{1}{\binom{p-1}{k-1}}\sum_{S\in\binom{[p]}{k}}\psi^{2}\biggr{\\}}^{1/2}\leq\|\tilde{v}\|_{2}\cdot\sqrt{\frac{p}{k}}\psi=\sqrt{\frac{p}{k}}\psi.$

Taking supremum over $v$ on the left-hand side, we arrive at the conclusion. ∎

Suppose that $X=(x_{1},\dots,x_{n})^{\top}$ is generated by independent
$x_{i}\sim N_{p}(0,\Sigma)$ with $n\geq p$ and some positive definite
$\Sigma$. Write $S=X^{\top}X$ and $S_{1}=X_{(0,t]}^{\top}X_{(0,t]}$. By Mitra
(1970), for any well-defined function $\phi:S\mapsto\phi(S)$ such that
$\phi(S)S\phi(S)^{\top}=I_{p}$, $U=\phi(S)S_{1}\phi(S)^{\top}$ is said to have
a matrix-variate Beta distribution, i.e.
$U\sim\mathrm{Beta}_{p}(t/2,(n-t)/2)$. To the best of our knowledge, it is
unclear what happens when $n<p$ in the literature. The following Lemma 9 and
Corollary 10 effectively generalize the existing matrix-variate Beta
distribution to the rank-deficient case of $n<p$.

######  Lemma 9.

Suppose $X\in\mathbb{R}^{n\times p}$ has independent $N(0,1)$ entries and
write $r:=\min\\{n,p\\}$. There exists an almost surely unique way of writing
$X=QR$ as its QR decomposition, where $Q\in\mathbb{O}^{n\times r}$ and
$R\in\mathbb{R}^{r\times p}$ such that $R_{i,j}=0$ for all $i>j$ and
$R_{i,i}\geq 0$ for all $i\in[n]$. We have that $Q$ and $R$ are independent
and $Q\sim\mathrm{Unif}(\mathbb{O}^{n\times n})$. Furthermore,
$R_{i,i}^{2}\sim\chi^{2}(n-i+1)$ and $R_{i,j}\sim N(0,1)$ for $i\in[n]$ and
$j\in[p]$ with $i<j$. Then
$B:=Q_{(0,t]}^{\top}Q_{(0,t]}\sim\operatorname{Beta}_{r}(t/2,(n-t)/2)$ and is
independent of $X^{\top}X$.

###### Proof.

First we consider the case of $n\geq p$. Write the (almost surely) unique QR
decomposition of $X$ by $X=QR$ with $Q\in\mathbb{O}^{n\times p}$ and
$R\in\mathbb{R}^{p\times p}$ being an upper triangular matrix with
$R_{i,i}\geq 0$ for all $i\in[p]$. For any fixed $H\in\mathbb{O}^{n\times n}$,
$HX\stackrel{{\scriptstyle\mathrm{d}}}{{=}}X$, whence
$HQR\stackrel{{\scriptstyle\mathrm{d}}}{{=}}QR$. As such, the joint density of
$Q$ and $R$ is constant for every possible value of $Q\in\mathbb{O}^{p\times
n}$, whence $Q$ and $R$ are independent and
$Q\sim\mathrm{Unif}(\mathbb{O}^{p\times n})$. By (Muirhead, 2009, Theorem
3.2.14), we have $R_{i,i}^{2}\sim\chi^{2}(n-i+1)$ and $R_{i,j}\sim N(0,1)$ and
$R_{i,j}$ are independent for all $i\leq j$. We define $S:=X^{\top}X$ and
$S_{1}=X_{(0,t]}^{\top}X_{(0,t]}$ and $S_{2}=X_{(t,n]}^{\top}X_{(t,n]}$.
Define $S^{1/2}:=R^{\top}$, and by Mitra (1970),

$B:=Q_{(0,t]}^{\top}Q_{(0,t]}=S^{-1/2}S_{1}(S^{-1/2})^{\top}\sim\operatorname{Beta}_{p}(t/2,(n-t)/2).$

We note that $B$ as a function of $Q$ is independent of $X^{\top}X=R^{\top}R$,
by the independence of $Q$ and $R$.

Now we consider the case $n<p$. Write $X_{j}$ as the $j$th column of $X$.
Write $X=[X_{\mathrm{L}}\mid X_{\mathrm{R}}]$ where
$X_{\mathrm{L}}=[X_{1}\mid\dots\mid X_{n}]$ and
$X_{\mathrm{R}}:=[X_{n+1}\mid\dots\mid X_{p}]$. For $X_{\mathrm{L}}$ whose
rank is almost surely $n$, there exists a unique QR decomposition such that
$X_{\mathrm{L}}=QR_{\mathrm{L}}$. Take
$R_{\mathrm{R}}:=Q^{\top}X_{\mathrm{R}}$ and $R:=[R_{\mathrm{L}}\mid
R_{\mathrm{R}}]$, and we have $X=QR$, where both $Q$ and $R$ are almost surely
unique. By the same argument as the case of $n\geq p$, we have $Q$ and $R$ are
independent and $Q\sim\mathrm{Unif}(\mathbb{O}^{n\times n})$. Applying the
conclusion from the case of $n\geq p$ on $X_{\mathrm{L}}=QR_{\mathrm{L}}$, we
have $(R_{L})_{i,i}\sim\chi^{2}(n-i+1)$ and $(R_{L})_{i,j}\sim N(0,1)$.
Furthermore, since $R_{\mathrm{R}}=Q^{\top}X_{\mathrm{R}}$ where both $Q$ and
$X_{\mathrm{R}}$ is independent of $X_{\mathrm{L}}$ and $Q$ is independent of
$X_{\mathrm{R}}$, all entries of $R_{\mathrm{R}}$ are standard normals
independent of $R_{\mathrm{L}}$.

Applying the case of $n\geq p$ on $X_{\mathrm{L}}=QR_{\mathrm{L}}$, we have
$B:=Q_{(0,t]}^{\top}Q_{(0,t]}\sim\operatorname{Beta}_{n}(t/2,(n-t)/2).$ By the
same argument as before, $B$ is independent of $X^{\top}X$. ∎

######  Corollary 10.

Let $X=(x_{1},\dots,x_{n})^{\top}$ where $x_{i}\sim N_{p}(0,\Sigma)$ where
$\Sigma$ is a positive definite matrix. Write $S=X^{\top}X$ and
$S_{1}:=X_{(0,t]}^{\top}X_{(0,t]}$. Let $r:=\min\\{n,p\\}$ and
$\phi:\mathbb{R}^{p\times p}\to\mathbb{R}^{p\times r}$ be a function such that
$\phi(S)S\phi(S)^{\top}=I_{r}$ for any positive semi-definite matrix
$S\in\mathbb{R}^{p\times p}$ of rank $r$. Then
$\phi(S)S_{1}\phi(S)^{\top}\sim\mathrm{Beta}_{r}(t/2,(n-t)/2)$ and is
independent of $S=X^{\top}X$.

###### Proof.

By the positive definiteness of $\Sigma$, we find
$\Sigma^{1/2}\in\mathbb{R}^{p\times p}$ such that
$\Sigma=\Sigma^{1/2}(\Sigma^{1/2})^{\top}$. Find $\Sigma^{-1/2}$ such that
$\Sigma^{-1/2}\Sigma^{1/2}=I_{p}$. Define $y_{i}:=\Sigma^{-1/2}x_{i}\sim
N_{p}(0,I_{p})$ and $Y=X(\Sigma^{-1/2})^{\top}$ has independent $N(0,1)$
entries. Let $\tilde{S}=Y^{\top}Y=\Sigma^{-1/2}S(\Sigma^{-1/2})^{\top}$ and
$\tilde{S}_{1}=Y_{(0,t]}^{\top}Y_{(0,t]}=\Sigma^{-1/2}S_{1}(\Sigma^{-1/2})^{\top}$,
whence by the definition of $\phi$ we have
$\phi(S)S\phi(S)^{\top}=\phi(S)\Sigma^{1/2}\tilde{S}(\Sigma^{1/2})^{\top}\phi(S)^{\top}=I_{r}$.
As such, we define
$\tilde{\phi}(\tilde{S}):=\phi(\Sigma^{1/2}\tilde{S}(\Sigma^{-1/2})^{\top})\Sigma^{1/2}$,
which is a well-defined function, and have
$\tilde{\phi}(\tilde{S})\tilde{S}\tilde{\phi}({\tilde{S}})^{\top}=I_{r}$.
Since $\tilde{\phi}(\tilde{S})\tilde{S}_{1}\tilde{\phi}(\tilde{S})^{\top}$ can
be defined by $Y$ with no dependence on $\Sigma$, it suffices to work on the
case $\Sigma=I_{p}$, which we assume for the rest of the proof.

By Lemma 9, write the unique QR decomposition of $X$ by $X=QR$ with
$Q\in\mathbb{O}^{n\times r}$ and $R\in\mathbb{R}^{r\times p}$, whence
$S=R^{\top}R$. Write $B=Q_{(0,t]}^{\top}Q_{(0,t]}$. Since
$\phi(S)^{\top}R^{\top}R\phi(S)=I_{r}$, $H:=R\phi(S)\in\mathbb{O}^{r\times
r}$, whence
$\phi(S)^{\top}S_{1}\phi(S)=H^{\top}Q_{(0,t]}^{\top}Q_{(0,t]}H=H^{\top}BH$.
Since $B$ as a function of $Q$ is independent of $S$ and $R$, it is
independent of $H$. By noting $H^{\top}BH=\phi(S)S_{1}\phi(S)^{\top}$ and
$B\stackrel{{\scriptstyle\mathrm{d}}}{{=}}H^{\top}BH$ because
$QH\stackrel{{\scriptstyle\mathrm{d}}}{{=}}Q$ for any $H\in\mathbb{O}^{r\times
r}$. The independence between $\phi(S)S_{1}\phi(S)$ and $S$ follows from the
fact that the distribution of $\phi(S)S_{1}\phi(S)$ is invariant conditionally
on $S$. ∎

Recall that $J_{(a_{1},a_{2}]}:=(J_{a_{1}+1},\cdots,J_{a_{2}})^{\top}$ is the
submatrix of $J$ by taking only the $(a_{1}+1)$-th to $a_{2}$-th rows for any
matrix $J$. For the rest of the paper, we define shorthand

$S_{a_{1},a_{2}}:=\sum_{i=a_{1}+1}^{a_{2}}x_{i}x_{i}^{\top}=X^{\top}_{(a_{1},a_{2}]}X_{(a_{1},a_{2}]}$

Define the scalar quantity
$\eta(n,p):=(\mathbb{E}[x_{1}x_{1}^{\top}(n^{-1}S_{0,n})^{-1}x_{n}x_{n}^{\top}])_{1,1}$.

######  Lemma 11.

For all $t\leq z$, $W_{t}^{\top}W_{z}=4S_{0,t}S_{0,n}^{-1}S_{z,n}$, whence for
$z\in[n]$ and $t\in[z]$ we have

$\displaystyle\mathbb{E}[W_{t}^{\top}W_{z}]=4t(n-z)n\eta(n,p)I_{p}.$

Furthermore, under Condition 2, $\frac{\eta(n,p)}{(n-p)n^{-1}}\to 1$, i.e.,
$\eta(n,p)\to\eta$ as $n,p\to\infty$.

###### Proof.

By the construction of $A$, we have
$AA^{\top}=I_{n}-X(X^{\top}X)^{-1}X^{\top}$. We have

$\displaystyle W_{t}^{\top}W_{z}$
$\displaystyle=\begin{pmatrix}X_{(0,t]}^{\top}&-X_{(t,n]}^{\top}\end{pmatrix}\begin{pmatrix}{A_{(0,t]}}\\\
{A_{(t,n]}}\end{pmatrix}\begin{pmatrix}A_{(0,z]}^{\top}&A_{(z,n]}^{\top}\end{pmatrix}\begin{pmatrix}{X_{(0,z]}}\\\
{-X_{(z,n]}}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}X_{(0,t]}^{\top}&-X_{(t,n]}^{\top}\end{pmatrix}(I_{n}-X(X^{\top}X)^{-1}X^{\top})\begin{pmatrix}{X_{(0,z]}}\\\
{-X_{(z,n]}}\end{pmatrix}$
$\displaystyle=(S_{0,t}-S_{t,z}+S_{z,n})-(S_{0,t}-S_{t,z}-S_{z,n})S_{0,n}^{-1}(S_{0,t}+S_{t,z}-S_{z,n})$
$\displaystyle=2(S_{0,t}-S_{t,z}-S_{z,n})S_{0,n}^{-1}S_{z,n}+2S_{z,n}$
$\displaystyle=4S_{0,t}S_{0,n}^{-1}S_{z,n}.$

In particular, we have for all $z\in[n]$ and $t\in[z]$

$\displaystyle\mathbb{E}W_{t}^{\top}W_{z}$
$\displaystyle=4\sum_{i=1}^{t}\sum_{j=z+1}^{n}\mathbb{E}[x_{i}x_{i}^{\top}S_{0,n}^{-1}x_{j}x_{j}^{\top}]=4t(n-z)\mathbb{E}[x_{1}x_{1}^{\top}S_{0,n}^{-1}x_{n}x_{n}^{\top}],$

where we invoke the exchangeability of
$x_{i}x_{i}^{\top}S_{0,n}^{-1}x_{j}x_{j}^{\top}$ for all $1\leq i<j\leq n$ in
the second equality.

We first note
$\mathbb{E}[W_{z}^{\top}W_{z}]=4z(n-z)\mathbb{E}[x_{1}x_{1}^{\top}S_{0,n}^{-1}x_{n}x_{n}^{\top}]$,
and then show that $\mathbb{E}[W_{z}^{\top}W_{z}]$ must be a scale of $I_{p}$.
Since for any $U\in\mathbb{O}^{p\times p}$,
$x_{i}^{\top}U\stackrel{{\scriptstyle\mathrm{d}}}{{=}}x_{i}^{\top}$ for
$i\in[n]$, whence we have
$U^{\top}W_{z}^{\top}W_{z}U=4U^{\top}S_{0,t}U(U^{\top}S_{0,n}U)^{-1}U^{\top}S_{z,n}U\stackrel{{\scriptstyle\mathrm{d}}}{{=}}W_{z}^{\top}W_{z}$.
In particular, $W_{z}^{\top}W_{z}$ have identically distributed diagonals and
identically distributed off-diagonals. It suffices to verify that its off-
diagonals have zero mean.

Let $X=QT$ be the almost surely unique QR decomposition of $X$, where we only
take non-negative diagonal entries in $T$. By Equation (15) of Gao and Wang
(2022), we have $W_{z}^{\top}W_{z}=4T^{\top}V\Lambda(I_{p}-\Lambda)V^{\top}T$,
where $V\Lambda V^{\top}=Q_{(0,z]}^{\top}Q_{(0,z]}$ is the eigendecomposition
of $B:=Q_{(0,t]}^{\top}Q_{(0,t]}$. Note that
$V\sim\mathrm{Unif}(\mathbb{O}^{p\times p})$, $\Lambda$ and $T$ are mutually
independent and $T$ has independent entries with $T_{j,j}=t_{j}>0$ such that
$t_{j}^{2}\sim\chi^{2}_{n-j+1}$ and $T_{j,k}=z_{jk}\sim N(0,1)$ for $j\neq k$.
For off-diagonals, it suffices to have

$\displaystyle(\mathbb{E}[W_{z}^{\top}W_{z}])_{1,2}$
$\displaystyle=4\mathbb{E}\biggl{[}\sum_{j=1}^{p}t_{1}V_{j,1}\lambda_{j}(1-\lambda_{j})(t_{2}V_{j,2}+z_{12}V_{j,1})\biggr{]}$
$\displaystyle=4\sum_{j=1}^{p}\biggl{[}\mathbb{E}[t_{1}]\mathbb{E}[\lambda_{j}(1-\lambda_{j})](\mathbb{E}[t_{2}]\mathbb{E}[V_{j,1}V_{j,2}]+\mathbb{E}[z_{12}]\mathbb{E}[V_{j,1}^{2}])\biggr{]}=0,$

where $\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{p}$ are diagonal
elements of $\Lambda$ and
$\mathbb{E}[V_{j,1}V_{j,2}]=(1/p)\sum_{j=1}^{p}V_{j,1}V_{j,2}=0$ since
$V\in\mathbb{O}^{p\times p}$.

Lastly, under Condition 2, by Proposition 8 of Gao and Wang (2022),
$(W_{z}^{\top}W_{z})_{1,1}\xrightarrow{\mathrm{a.s.}}(4z(n-z)(n-p))/n^{3}$.
Noting $(\mathbb{E}[W_{z}^{\top}W_{z}])_{1,1}=4z(n-z)n\eta(n,p)$, we conclude
the convergence of $\eta(n,p)\to\eta$. ∎

######  Lemma 12.

Let
$x_{1},\ldots,x_{n}\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}N_{p}(0,I_{p})$.
Fix $z\in[n]$ and $t\in[z]$. For any nonrandom $u,w\in\mathcal{S}^{p-1}$, we
have

$\mathbb{P}\biggl{[}\frac{1}{t}\biggl{|}u^{\top}\biggl{\\{}S_{0,t}-\frac{t}{z}S_{0,z}\biggr{\\}}S_{0,n}^{-1}S_{z,n}w\biggr{|}\geq
x\biggm{|}S_{0,z},S_{z,n}\biggr{]}\leq
52\exp\biggl{\\{}-\frac{t^{2}x^{2}}{8z\|S_{0,n}^{-1}S_{z,n}\|_{\mathrm{op}}^{2}\|S_{0,z}/z\|_{\mathrm{op}}^{2}}\biggr{\\}}.$

###### Proof.

For notational simplicity, we use $\mathbb{P}^{z}$ and $\mathbb{E}^{z}$ to
denote the conditional probability and expectation with respect to the
$\sigma$-algebra generated by $(S_{0,z},S_{z,n})$. Note that
$S_{0,t}\mid(S_{0,z},S_{z,n})\stackrel{{\scriptstyle\mathrm{d}}}{{=}}S_{0,t}\mid
S_{0,z}$ and $\mathbb{E}^{z}[S_{0,t}]=(t/z)S_{0,z}$.

Let $r:=\min\\{z,p\\}$. We note $S_{0,z}$ has (almost surely) rank $r$, and
write $S_{0,z}=R^{\top}R$ for the (almost surely unique) Cholesky
decomposition of $S_{0,z}$ such that $R\in\mathbb{R}^{r\times p}$ is an upper-
triangular matrix with positive diagonal entries. Write
$R^{\dagger}\in\mathbb{R}^{p\times r}$ for the (almost surely unique)
Moore–Penrose pseudo-inverse of $R$ such that $RR^{\dagger}=I_{p}$. By
Corollary 10, the matrix
$B:=(R^{\dagger})^{\top}S_{0,t}R^{\dagger}\sim\mathrm{Beta}_{r}(t/2,(z-t)/2)$
has a matrix-variate Beta distribution, and is independent of $S_{0,z}$ with
$(t/z)I_{p}$ as its (conditional) mean. Observe that
$R^{\top}BR=(R^{\dagger}R)^{\top}S_{0,t}(R^{\dagger}R)=S_{0,t}$, since
$R^{\dagger}R$ is a (symmetric) orthogonal projection matrix onto the row
space of $X_{(0,z]}$, which contains the row space of $X_{(0,t]}$.

Define $v:=S_{0,n}^{-1}S_{z,n}w/\|S_{0,n}^{-1}S_{z,n}w\|_{2}$. Writing
$\tilde{u}:=Ru/\sqrt{z}$ and $\tilde{v}:=Rv/\sqrt{z}$, we have

$\displaystyle\mathbb{P}^{z}\biggl{[}\frac{1}{t}|u^{\top}\\{S_{0,t}-\mathbb{E}^{z}(S_{0,t})\\}$
$\displaystyle S_{0,n}^{-1}S_{z,n}w|\geq
x\biggr{]}\leq\mathbb{P}^{z}\biggl{[}\frac{1}{t}|u^{\top}\\{S_{0,t}-\mathbb{E}^{z}(S_{0,t})\\}v|\geq\frac{x}{\|S_{0,n}^{-1}S_{0,z}\|_{\mathrm{op}}}\biggr{]}$
$\displaystyle=\mathbb{P}^{z}\biggl{\\{}\biggl{|}\tilde{u}^{\top}\biggl{(}\frac{z}{t}B-I_{p}\biggr{)}\tilde{v}\biggr{|}\geq\frac{x}{\|S_{0,n}^{-1}S_{0,z}\|_{\mathrm{op}}}\biggr{\\}}$
$\displaystyle\leq\mathbb{P}^{z}\biggl{\\{}\biggl{|}\biggl{(}\frac{\tilde{u}}{\|\tilde{u}\|_{2}}\biggr{)}^{\top}\biggl{(}\frac{z}{t}B-I_{p}\biggr{)}\frac{\tilde{v}}{\|\tilde{v}\|_{2}}\biggr{|}\geq\frac{x}{\|S_{0,n}^{-1}S_{0,z}\|_{\mathrm{op}}\|S_{0,z}/z\|_{\mathrm{op}}}\biggr{\\}}.$

Write shorthand
$\psi:=\|S_{0,n}^{-1}S_{0,z}\|_{\mathrm{op}}\|S_{0,z}/z\|_{\mathrm{op}}$,
which is measurable with respect to the $\sigma$-algebra generated by
$(S_{0,z},S_{z,n})$. There exists an orthogonal matrix $U$ such that
$U\tilde{u}/\|\tilde{u}\|_{2}=e_{1}$ and $U\tilde{v}/\|\tilde{v}\|_{2}=\alpha
e_{1}+\beta e_{2}$ for real $\alpha$ and $\beta$ such that
$\alpha^{2}+\beta^{2}=1$, where $e_{j}$ denotes the $j$th standard basis
vector in $\mathbb{R}^{p}$. Using the fact that
$B\mid(S_{0,z},S_{z,n})\stackrel{{\scriptstyle\mathrm{d}}}{{=}}U^{\top}BU\mid(S_{0,z},S_{z,n})$,
we have for $J=\\{1,2\\}$ that

$\displaystyle\mathbb{P}^{z}\biggl{\\{}\biggl{|}\biggl{(}\frac{\tilde{u}}{\|\tilde{u}\|_{2}}\biggr{)}^{\top}\biggl{(}\frac{z}{t}B$
$\displaystyle-
I_{p}\biggr{)}\frac{\tilde{v}}{\|\tilde{v}\|_{2}}\biggr{|}\geq\frac{x}{\psi}\biggr{\\}}$
$\displaystyle\leq\mathbb{P}^{z}\biggl{\\{}\biggl{|}\frac{z}{t}B_{1,1}-1\biggr{|}\geq\frac{x}{\sqrt{2}\psi}\biggr{\\}}+\mathbb{P}^{z}\biggl{\\{}\frac{z}{t}|B_{1,2}|\geq\frac{x}{\sqrt{2}\psi}\biggr{\\}}$
$\displaystyle\leq
2\mathbb{P}^{z}\biggl{\\{}\biggl{\|}\biggl{(}B-\frac{t}{z}I_{p}\biggr{)}_{J,J}\biggr{\|}_{\mathrm{op}}\geq\frac{tx}{\sqrt{2}z\psi}\biggr{\\}},$

where the first inequality holds by noting
$|\alpha|+|\beta|\leq(2\alpha^{2}+2\beta^{2})^{1/2}=\sqrt{2}$. Note that
$\\{w\in\mathcal{S}^{p-1}:\mathrm{supp}(w)\subseteq J\\}$ is isomorphic to
$S^{1}$, which contains a $(1/4)$-net $\mathcal{N}$ of cardinality
$\lceil\frac{2\pi}{4\arcsin(1/8)}\rceil=13$. By Gupta and Nagar (1999, Theorem
5.3.12), for each $w\in\mathcal{N}$, we have
$w^{\top}Bw\sim\mathrm{Beta}(t/2,(z-t)/2)$. Hence, by Vershynin (2012, Lemma
5.4) and a union bound, we have

$\displaystyle\mathbb{P}^{z}\biggl{\\{}\biggl{\|}\biggl{(}B-\frac{t}{z}I_{p}\biggr{)}_{J,J}\biggr{\|}_{\mathrm{op}}\geq\frac{tx}{\sqrt{2}z\psi}\biggr{\\}}$
$\displaystyle\leq\mathbb{P}^{z}\biggl{\\{}\sup_{w\in\mathcal{N}}\bigl{|}w^{\top}Bw-t/z\bigr{|}\geq\frac{tx}{2\sqrt{2}z\psi}\biggr{\\}}$
$\displaystyle\leq
13\mathbb{P}^{z}\biggl{\\{}\bigl{|}B_{1,1}-t/z\bigr{|}\geq\frac{tx}{2\sqrt{2}z\psi}\biggr{\\}}$
$\displaystyle\leq
26\exp\biggl{\\{}-\frac{t^{2}x^{2}}{8z\psi^{2}}\biggr{\\}},$

where we have used Marchal and Arbel (2017, Theorem 2.1) in the final
inequality. ∎

######  Lemma 13.

Let $X=(X_{1},\ldots,X_{p})^{\top}$ be uniformly distributed on the sphere
$\mathcal{S}^{p-1}$. Then, for $\delta\geq e^{-p/16}$, we have

$\mathbb{P}\biggl{(}\|X\|_{\infty}>\sqrt{\frac{4\log(2/\delta)}{p}}\biggr{)}\leq
p\delta.$

###### Proof.

Let $Z_{1},\ldots,Z_{p}$ be independent $N(0,1)$ random variables, then
$X_{1}\stackrel{{\scriptstyle\text{d}}}{{=}}Z_{1}/(Z_{1}^{2}+\cdots+Z_{p}^{2})^{1/2}$.
By a standard Gaussian tail bound, we have

$\mathbb{P}\bigl{\\{}Z_{1}>\sqrt{2\log(1/\delta)}\bigr{\\}}\leq\delta.$

Moreover, $\sum_{j=1}^{p}Z_{j}^{2}\sim\chi^{2}_{p}$. Since we have $\delta\geq
e^{-p/16}$, by Laurent and Massart (2000, Lemma 1),

$\mathbb{P}\biggl{(}\sum_{j=1}^{p}Z_{j}^{2}<\frac{p}{2}\biggr{)}\leq\mathbb{P}\biggl{(}\sum_{j=1}^{p}Z_{j}^{2}<p-2\sqrt{p\log(1/\delta)}\biggr{)}\leq\delta.$

The result follows by combining the above two bounds and applying a union
bound. ∎

Recall $B_{0}(k)\subseteq\mathbb{R}^{p}$ is the set of $k$-sparse unit
vectors.

######  Lemma 14.

For any $A\in\mathbb{R}^{p\times p}$, and any $\epsilon\in(0,1)$, there exists
an $\epsilon$-net $\mathcal{N}_{\epsilon}$ of $B_{0}(k)$ of cardinality at
most $\\{(1+2/\epsilon)ep/k\\}^{k}$ such that

$\sup_{u\in
B_{0}(k)}u^{\top}Av\leq(1-\epsilon)^{-1}\max_{u\in\mathcal{N}_{\epsilon}}u^{\top}Av.$

###### Proof.

By Vershynin (2012, Lemma 5.2), for each subset $S\subseteq[p]$ of cardinality
$k$, there exists an $\epsilon$-net $\mathcal{N}_{S}$ of $\\{v\in
B_{0}(k):\mathrm{supp}(v)\subseteq S\\}$ of cardinality at most
$(1+2/\epsilon)^{k}$. Define
$\mathcal{N}_{\epsilon}:=\cup_{S\subseteq[p]:|S|=k}\mathcal{N}_{S}$, then
$|\mathcal{N}_{\epsilon}|\leq\binom{p}{k}(1+2\epsilon)^{k}\leq\\{(1+2\epsilon)ep/k\\}^{k}$.
For any fixed $x\in B_{0}(k)$, find $\tilde{x}\in\mathcal{N}_{\epsilon}$ such
that $\|x-\tilde{x}\|_{2}\leq\epsilon$ and $\|x-\tilde{x}\|_{0}\leq k$. Thus,

$\displaystyle x^{\top}Av$
$\displaystyle=(x-\tilde{x})^{\top}Av+\tilde{x}^{\top}Av$
$\displaystyle\leq\|x-\tilde{x}\|_{2}\sup_{u\in
B_{0}(k)}u^{\top}Av+\sup_{u\in\mathcal{N}_{\epsilon}}u^{\top}Av$
$\displaystyle\leq\epsilon\sup_{u\in
B_{0}(k)}u^{\top}Av+\sup_{u\in\mathcal{N}_{\epsilon}}u^{\top}Av.$

The desired result follows by taking supremum over $x\in B_{0}(k)$ above. ∎

## References

  * Bai (1997) Bai, J. (1997) Estimation of a change point in multiple regression models.  _Review of Economics and Statistics_ , 79, 551–563.
  * Bai (2010) Bai, J. (2010) Common breaks in means and variances for panel data.  _J. Econom._ , 157, 78–92.
  * Bai and Perron (1998) Bai, J. and Perron, P. (1998) Estimating and testing linear models with multiple structural changes.  _Econometrica_ , 47–78.
  * Baranowski et al. (2019) Baranowski, R., Chen, Y. and Fryzlewicz, P. (2019) Narrowest-over-threshold detection of multiple change points and change-point-like features.  _J. Roy. Statist. Soc., Ser. B_ , 81, 649–672.
  * Bitter et al. (2020) Bitter, E. E., Townsend, M. H., Erickson, R., Allen, C. and O’Neill, K. L. (2020) Thymidine kinase 1 through the ages: a comprehensive review.  _Cell & Bioscience_, 10, 1–16.
  * Bühlmann and van de Geer (2011) Bühlmann, P. and van de Geer, S. (2011) _Statistics for High-dimensional Data: Methods, Theory and Applications_.  Springer Science & Business Media.
  * Chen et al. (2022) Chen, Y., Wang, T. and Samworth, R. J. (2022) High-dimensional, multiscale online changepoint detection.  _J. Roy. Statist. Soc., Ser. B_ , 84, 234–266.
  * Cho (2016) Cho, H. (2016) Change-point detection in panel data via double CUSUM statistic.  _Electron. J. of Statist._ , 10, 2000–2038.
  * Cho and Fryzlewicz (2014) Cho, H. and Fryzlewicz, P. (2014) Multiple-change-point detection for high dimensional time series via sparsified binary segmentation.  _J. Roy. Statist. Soc., Ser. B_ , 77, 475–507.
  * Efron et al. (2004) Efron, B., Hastie, T., Johnstone, I. and Tibshirani, R. (2004) Least angle regression.  _Ann. Statist._ , 32, 407–499.
  * Enikeeva and Klopp (2021) Enikeeva, F. and Klopp, O. (2021) Change-point detection in dynamic networks with missing links.  _arXiv preprint_ , arXiv:2106.14470.
  * Fan and Lv (2010) Fan, J. and Lv, J. (2010) A selective overview of variable selection in high dimensional feature space.  _Statist. Sinica_ , 20, 101.
  * Follain et al. (2022) Follain, B., Wang, T. and Samworth, R. J. (2022) High-dimensional changepoint estimation with heterogeneous missingness.  _J. Roy. Statist. Soc., Ser. B_ , to appear.
  * Fryzlewicz (2014) Fryzlewicz, P. (2014) Wild binary segmentation for multiple change-point detection.  _Ann. Statist._ , 42, 2243–2281.
  * Fryzlewicz (2018) Fryzlewicz, P. (2018) Tail-greedy bottom-up data decompositions and fast multiple change-point detection.  _Ann. Statist._ , 46, 3390–3421.
  * Fryzlewicz (2020) Fryzlewicz, P. (2020) Detecting possibly frequent change-points: Wild Binary Segmentation 2 and steepest-drop model selection.  _J. Korean Statist. Soc._ , 49, 1027–1070.
  * Gao and Wang (2022) Gao, F. and Wang, T. (2022) Two-sample testing of high-dimensional linear regression coefficients via complementary sketching.  _Ann. Statist._ , 50, 2950 – 2972.  URL https://doi.org/10.1214/22-AOS2216.
  * Guo et al. (2016) Guo, L., Zhang, S., Zhang, B., Chen, W., Li, X., Zhang, W., Zhou, C., Zhang, J., Ren, N. and Ye, Q. (2016) Silencing GTSE-1 expression inhibits proliferation and invasion of hepatocellular carcinoma cells.  _Cell Biology and Toxicology_ , 32, 263–274.
  * Gupta and Nagar (1999) Gupta, A. K. and Nagar, D. K. (1999) _Matrix Variate Distributions_ , volume 104.  CRC Press.
  * Hatleberg and Hinman (2021) Hatleberg, W. L. and Hinman, V. F. (2021) Modularity and hierarchy in biological systems: Using gene regulatory networks to understand evolutionary change.  In S. F. Gilbert (ed.), _Evolutionary Developmental Biology_ , volume 141 of _Current Topics in Developmental Biology_ , 39–73. Academic Press.
  * Horváth and Hušková (2012) Horváth, L. and Hušková, M. (2012) Change-point detection in panel data.  _J. Time Series Anal._ , 33, 631–648.
  * Jirak (2015) Jirak, M. (2015) Uniform change point tests in high dimension.  _Ann. Statist._ , 43, 2451–2483.
  * Julious (2001) Julious, S. A. (2001) Inference and estimation in a changepoint regression problem.  _J. Roy. Statist. Soc., Ser. D_ , 50, 51–61.
  * Kalman et al. (2004) Kalman, L., Lindegren, M. L., Kobrynski, L., Vogt, R., Hannon, H., Howard, J. T. and Buckley, R. (2004) Mutations in genes required for T-cell development: IL7R, CD45, IL2RG, JAK3, RAG1, RAG2, ARTEMIS, and ADA and severe combined immunodeficiency: HuGE review.  _Genetics in Medicine_ , 6, 16–26.
  * Kaul et al. (2019) Kaul, A., Jandhyala, V. K. and Fotopoulos, S. B. (2019) An efficient two step algorithm for high dimensional change point regression models without grid search.  _J. Mach. Learn. Res._ , 20, (111), 1–40.
  * Kovács et al. (2020) Kovács, S., Li, H., Bühlmann, P. and Munk, A. (2020) Seeded binary segmentation: A general methodology for fast and optimal change point detection.  _arXiv preprint arXiv:2002.06633_.
  * Kraft and Hunter (2009) Kraft, P. and Hunter, D. J. (2009) Genetic risk prediction—are we there yet?  _N. Engl. J. Medcine_ , 360, 1701–1703.
  * Krishnamurthy et al. (2019) Krishnamurthy, A., Mazumdar, A., McGregor, A. and Pal, S. (2019) Sample complexity of learning mixture of sparse linear regressions.  _Adv. Neur. Inform. Proc. Sys._ , 32.
  * Laurent and Massart (2000) Laurent, B. and Massart, P. (2000) Adaptive estimation of a quadratic functional by model selection.  _Ann. Statist._ , 28, 1302–1338.
  * Lee et al. (2016) Lee, S., Seo, M. H. and Shin, Y. (2016) The lasso for high dimensional regression with a possible change point.  _J. Roy. Statist. Soc., Ser. B_ , 78, 193–210.
  * Leonardi and Bühlmann (2016) Leonardi, F. and Bühlmann, P. (2016) Computationally efficient change point detection for high-dimensional regression.  _arXiv preprint arXiv:1601.03704_.
  * Lévy-Leduc and Roueff (2009) Lévy-Leduc, C. and Roueff, F. (2009) Detection and localization of change-points in high-dimensional network traffic data.  _Ann. Appl. Statist._ , 637–662.
  * Liu et al. (2021) Liu, H., Gao, C. and Samworth, R. J. (2021) Minimax rates in sparse, high-dimensional change point detection.  _Ann. Statist._ , 49, 1081–1112.
  * Marchal and Arbel (2017) Marchal, O. and Arbel, J. (2017) On the sub-Gaussianity of the Beta and Dirichlet distributions.  _Electron. Commun. Probab._ , 22, 1–14.
  * Mills et al. (1992) Mills, G. B., Schmandt, R., McGill, M., Amendola, A., Hill, M., Jacobs, K., May, C., Rodricks, A.-M., Campbell, S. and Hogg, D. (1992) Expression of TTK, a novel human protein kinase, is associated with cell proliferation.  _Journal of Biological Chemistry_ , 267, 16000–16006.
  * Mitra (1970) Mitra, S. K. (1970) A density-free approach to the matrix variate beta distribution.  _Sankhyā: The Indian Journal of Statistics, Series A_ , 81–88.
  * Muirhead (2009) Muirhead, R. J. (2009) _Aspects of Multivariate Statistical Theory_.  John Wiley & Sons, Hoboken, New Jersey.
  * Page (1955) Page, E. S. (1955) A test for a change in a parameter occurring at an unknown point.  _Biometrika_ , 42, 523–527.
  * Rand (1971) Rand, W. M. (1971) Objective criteria for the evaluation of clustering methods.  _J. Amer. Statist. Assoc._ , 66, 846–850.
  * Rinaldo et al. (2021) Rinaldo, A., Wang, D., Wen, Q., Willett, R. and Yu, Y. (2021) Localizing Changes in High-Dimensional Regression Models.  In A. Banerjee and K. Fukumizu (eds.), _Proceedings of The 24th International Conference on Artificial Intelligence and Statistics_ , volume 130 of _Proceedings of Machine Learning Research_ , 2089–2097. PMLR.
  * Smith (1985) Smith, R. L. (1985) Maximum likelihood estimation in a class of nonregular cases.  _Biometrika_ , 72, 67–90.
  * Städler et al. (2010) Städler, N., Bühlmann, P. and van de Geer, S. (2010) $\ell_{1}$-penalization for mixture regression models.  _TEST_ , 19, 209–256.
  * Suo et al. (2023) Suo, C., Polanski, K., Dann, E., Lindeboom, R. G., Vilarrasa-Blasi, R., Vento-Tormo, R., Haniffa, M., Meyer, K. B., Dratva, L. M., Tuong, Z. K. et al. (2023) Dandelion uses the single-cell adaptive immune receptor repertoire to explore lymphocyte developmental origins.  _Nature Biotechnology_ , 1–12.
  * Tibshirani (1996) Tibshirani, R. (1996) Regression shrinkage and selection via the lasso.  _J. Roy. Statist. Soc., Ser. B_ , 58, 267–288.
  * Vershynin (2012) Vershynin, R. (2012) Introduction to the non-asymptotic analysis of random matrices. In Y. Eldar and G. Kutyniok (Eds.).  _Compressed Sensing, Theory and Applications_ , 210–268.
  * Verzelen et al. (2020) Verzelen, N., Fromont, M., Lerasle, M. and Reynaud-Bouret, P. (2020) Optimal change-point detection and localization.  _arXiv preprint arXiv:2010.11470_.
  * Vostrikova (1981) Vostrikova, L. Y. (1981) Detecting “disorder” in multidimensional random processes.  In _Doklady Akademii Nauk_ , volume 259(2), 270–274. Russian Academy of Sciences.
  * Wainwright (2019) Wainwright, M. J. (2019) _High-Dimensional Statistics: A Non-Asymptotic Viewpoint_.  Cambridge University Press, Cambridge.
  * Wang et al. (2021) Wang, D., Yu, Y. and Rinaldo, A. (2021) Optimal covariance change point localization in high dimensions.  _Bernoulli_ , 27, 554–575.
  * Wang and Samworth (2018) Wang, T. and Samworth, R. J. (2018) High dimensional change point estimation via sparse projection.  _J. Roy. Statist. Soc., Ser. B_ , 80, 57–83.
  * Yin et al. (2018) Yin, D., Pedarsani, R., Chen, Y. and Ramchandran, K. (2018) Learning mixtures of sparse linear regressions using sparse graph codes.  _IEEE Trans. Inform. Theory_ , 65, 1430–1451.
  * Zhang et al. (2019) Zhang, L., Du, Y., Xu, S., Jiang, Y., Yuan, C., Zhou, L., Ma, X., Bai, Y., Lu, J. and Ma, J. (2019) DEPDC1, negatively regulated by miR-26b, facilitates cell proliferation via the up-regulation of FOXM1 expression in TNBC.  _Cancer Letters_ , 442, 242–251.
  * Zhang et al. (2010) Zhang, N. R., Siegmund, D. O., Ji, H. and Li, J. Z. (2010) Detecting simultaneous changepoints in multiple sequences.  _Biometrika_ , 97, 631–645.
  * Zhou et al. (2018) Zhou, H., Cai, L., Zhang, X., Li, A., Miao, Y., Li, Q., Qiu, X. and Wang, E. (2018) ARHGEF39 promotes tumor progression via activation of Rac1/P38 MAPK/ATF2 signaling and predicts poor prognosis in non-small cell lung cancer patients.  _Laboratory Investigation_ , 98, 670–681.