# Estimation of high-dimensional change-points under a group sparsity
structure

Hanqing Cai and Tengyao Wang111Research supported by EPSRC grant EP/T02772X/1.
University College London

###### Abstract

Change-points are a routine feature of ‘big data’ observed in the form of
high-dimensional data streams. In many such data streams, the component series
possess group structures and it is natural to assume that changes only occur
in a small number of all groups. We propose a new change point procedure,
called groupInspect, that exploits the group sparsity structure to estimate a
projection direction so as to aggregate information across the component
series to successfully estimate the change-point in the mean structure of the
series. We prove that the estimated projection direction is minimax optimal,
up to logarithmic factors, when all group sizes are of comparable order.
Moreover, our theory provide strong guarantees on the rate of convergence of
the change-point location estimator. Numerical studies demonstrates the
competitive performance of groupInspect in a wide range of settings and a real
data example confirms the practical usefulness of our procedure.

##  1 Introduction

Modern applications routinely generate time-ordered high-dimensional datasets,
where many covariates are simultaneously measured over time. Examples include
wearable technologies recording the health state of individuals from multi-
sensor feedbacks (Hanlon and Anderson, 2009), internet traffic data collected
by tens of thousands of routers (Peng, Leckie and Ramamohanarao, 2004) and
functional Magnetic Resonance Imaging (fMRI) scans that record the time
evolution of blood oxygen level dependent (BOLD) chemical contrast in
different areas of the brain (Aston and Kirch, 2012). The explosion in number
of such high-dimensional data streams calls for methodological advances for
their analysis.

Change-point analysis is an essential statistical technique used in
identifying abrupt changes in a time series. Time points at which such abrupt
change occurs are called ‘change-points’. Through estimating the location of
change-points, we can divide the time series into shorter segments that can be
analysed using methods designed for stationary time series. Moreover, in many
applications, the estimated change-points indicate specific events that are
themselves of great interest. In the examples mentioned in the previous
paragraph, they can be used to raise alarms about abnormal health events,
detect distributed denial of service attacks on the network and pinpoint the
onset of certain brain activities.

Classical change-point analysis focuses on univariate time series. The current
state-of-art methods including Killick, Fearnhead and Eckley (2012); Frick,
Munk and Sieling (2014); Fryzlewicz (2014). However, classical univariate
change-point methods are often inadequate for high-dimensional datasets that
are routinely encountered in modern applications. When applied componentwise,
they are often sub-optimal as signals can spread over many components. As a
result, several new methodologies have been proposed to test and estimate
change-points in the high-dimensional settings. These include methods that
apply a simple $\ell_{2}$ or $\ell_{\infty}$ aggregation of test statistics
across different components (Horváth and Hušková, 2012; Jirak, 2015), and more
complex methods such as a scan-statistics based approach by Enikeeva and
Harchaoui (2019), the Sparsified Binary Segmentation algorithm by Cho and
Fryzlewicz (2015), the double CUSUM algorithm of Cho (2016) and a projection-
based approach by Wang and Samworth (2018).

To get around the issue of the curse of dimensionality, existing high-
dimensional change-point methods often assume that the signal of change
possesses some form of sparsity. For example, in the high-dimensional mean
change setting studied in Jirak (2015); Cho and Fryzlewicz (2015); Wang and
Samworth (2018); Enikeeva and Harchaoui (2019), it is assumed that the
difference in mean before and after a change-point is nonzero only in a small
subset of coordinates. While the sparsity assumption greatly reduces the
complexity of the original high-dimensional problem, it often does not capture
the the full extent of the structure in the vector of change available in real
data applications. For instance, in many applications, the coordinates of the
high-dimensional vectors are naturally clustered into groups and coordinates
within the same group tend to change together. At each change-point, only a
small number of groups will undergo a change. Such a group sparsity change-
point structure is useful in modelling many practical applications. Examples
include financial data stream where changes are often grouped by industry
sectors and a small number of sectors may experience virtually simultaneous
market shocks. Also, in functional magnetic resonance imaging data, voxels
belonging to the same brain functional regions tend to change simultaneously
over time. Similar group sparsity assumptions have been made in other
statistical problems including Yuan and Lin (2006); Wang and Leng (2008);
Simon et al. (2020).

In this work, we provide a new high-dimensional change-point methodology that
exploits the group sparsity structure of the changes. More precisely, given
pre-specified grouping information of all the coordinates, our algorithm,
named groupInspect (standing for group-based informative sparse projection
estimator of change-points), will first estimate a vector of projection that
is closely aligned with the true vector of change at each change-point. It
will then project the high-dimensional data series along this estimated
direction and apply a univariate change-point method on the projected series
to identify the location of the change. The above procedure can be combined
with a wild binary-segmentation algorithm (Fryzlewicz, 2014) to recursively
identify multiple change-points. We show that, in a single change-point
setting, the projection direction estimator employed in groupInspect has a
minimax optimal dependence, up to logarithmic factors, on both the $\ell_{0}$
sparsity parameter and the group-sparsity parameter, representing respectively
the number of nonzero elements and the number of nonzero groups in the vector
of change. Furthermore, groupInspect achieves a
$\sqrt{\log(n)}/(n\vartheta^{2})$ rate of convergence for the estimated
location of a single change-point, where $\vartheta$ denotes the $\ell_{2}$
norm of the vector of change, which up to logarithmic factors is minimax
optimal.

The outline of the paper is as follows. In Section 2, we describe the formal
setup of our problem. The groupInspect methodology is then introduced in
Section 3, with its theoretical performance guarantees provided in Section 4.
We illustrate the empirical performance of groupInspect via simulatinos and a
real-data example in Section 5. Proofs of all theoretical results are deferred
to Section 6, and ancillary results and their proofs are given in Section 7.

###  1.1 Notation

For any positive integer $n$, we write $[n]=\\{1,\dots,n\\}$. For a vector
$v=(v_{1},\ldots,v_{n})^{\top}\in\mathbb{R}^{n}$, we define
$\|v\|_{0}=\sum_{i=1}^{n}\mathbbm{1}_{\\{v_{i}\neq 0\\}}$,
$\|v\|_{\infty}=\max_{i\in[n]}|v_{i}|$ and
$\|v\|_{q}=\bigl{\\{}\sum_{i=1}^{n}(v_{i})^{q}\bigr{\\}}^{1/q}$ for any
positive integer $q$, and let
$\mathbb{S}^{n-1}=\\{v\in\mathbb{R}^{n}:\|v\|_{2}=1\\}$. For a matrix
$A\in\mathbb{R}^{p\times n}$, we write $\|A\|_{*}$ for its nuclear norm and
write $\|A\|_{\mathrm{F}}$ for its Frobenius norm.

For any $S\subseteq[n]$, we write $v_{S}$ for the $|S|$-dimensional vector
obtained by extracting coordinates of $v$ in $S$. For a matrix
$A\in\mathbb{R}^{p\times n}$, $J\in[p]$ and $S\in[n]$, we write $A_{J,S}$ for
the submatrix obtained by extracting rows and columns of $A$ indexed by $J$
and $S$ respectively. When $S=[n]$, we abbreviate $A_{J,[n]}$ by $A_{J}$. When
$S=\\{t\\}$ is a single element set, we slightly abuse notation and write
$A_{J,t}$ instead of $A_{J,\\{t\\}}$.

Given two sequences $(a_{n})_{n\in\mathbb{N}}$ and $(b_{n})_{n\in\mathbb{N}}$
such that $a_{n},b_{n}>0$ for all $n$, we write $a_{n}\lesssim b_{n}$ (or
equivalently $b_{n}\gtrsim a_{n}$) if $a_{n}\leqslant Cb_{n}$ for some
universal constant $C$.

##  2 Problem description

Let $X_{1},\ldots,X_{n}$ be independent random vectors with distribution:

$X_{t}\sim N_{p}(\mu_{t},\sigma^{2}I_{p}),\quad 1\leqslant t\leqslant n,$ (1)

which we can combine into a single data matrix $X\in\mathbb{R}^{p\times n}$.
We assume that the sequence of mean vectors $(\mu_{t})_{t=1}^{n}$ undergoes
changes at times $z_{i}\in\\{1,\ldots,n-1\\}$ for $i\in\\{1,\ldots,\nu\\}$, in
the sense that

$\mu_{z_{i}+1}=\cdots=\mu_{z_{i+1}}=:\mu^{(i)},\quad\forall\,i\in\\{0,\ldots,\nu\\},$
(2)

where we use the convention that $z_{0}=0$ and $z_{\nu+1}=n$. We assume that
consecutive change-points are sufficiently separated in the sense that

$\min\\{z_{i+1}-z_{i}:0\leqslant i\leqslant\nu\\}\geqslant n\tau.$

Suppose further that each of the $p$ coordinates belong to (at least) one of
$G$ groups. Specifically, let $\mathcal{J}_{g}$ denotes the set of indices
associated with the $g$th group for $g\in\\{1,\ldots,G\\}$, we have that

$\bigcup_{g=1}^{G}\mathcal{J}_{g}=[p].$ (3)

We assume that coordinates in the same group will tend to change together. We
will consider both the case of overlapping and non-overlapping groups. In the
latter scenario, each coordinate belongs to a unique group and
$(\mathcal{J}_{g})_{g\in[G]}$ forms a partition of $[p]$.

Our goal is to estimate the locations of change $z_{1},\ldots,z_{\nu}$ from
the data matrix $X$ and the pre-specified grouping information
$(\mathcal{J}_{g})_{g\in[G]}$. Motivated by Wang and Samworth (2018), the best
way to aggregate the component series so as to maximise the signal-to-noise
ratio around the $i$th change-point is to project the data along a direction
close to the vector of change $\theta^{(i)}=\mu^{(i)}-\mu^{(i-1)}$. Let
$v^{(i)}$ be the unit vector parallel to $\theta^{(i)}$:

$v^{(i)}=\theta^{(i)}/\|\theta^{(i)}\|_{2},$

which we will call the oracle direction for the $i$th change-point. We measure
the quality of any estimated projection direction $\hat{v}$ with the
Davis–Kahan sin $\theta$ loss (Davis and Kahan, 1970)

$L(\hat{v},v^{(i)})=\sqrt{1-(\hat{v}^{\top}v^{(i)})^{2}}$

and measure the quality of the subsequent location estimator $\hat{z}_{i}$ by
$\mathbb{E}|\hat{z}_{i}-z_{i}|$.

The difficulty of the estimation task depends on both the noise level $\sigma$
and the vector of change $\theta^{(i)}=\mu^{(i)}-\mu^{(i-1)}$. More precisely,
we assume that the change is localised in a small number of the $G$ groups as
defined in (3). Define $\phi:\mathbb{R}^{p}\to\mathbb{R}^{G}$ such that
$\phi(x)=(\|x_{\mathcal{J}_{1}}\|_{2},\|x_{\mathcal{J}_{2}}\|_{2},\ldots,\|x_{\mathcal{J}_{G}}\|_{2})^{\top}$,
we assume that

$\|\theta^{(i)}\|_{0}\leqslant k,\quad\|\phi(\theta^{(i)})\|_{0}\leqslant
s\quad\text{and}\quad\|\theta^{(i)}\|_{2}\geqslant\vartheta.$ (4)

##  3 Methodology

###  3.1 Single change-point estimation

Initially, we will consider estimation of a single change-point, where
$\nu=1$. This can be extended to estimate multiple change-points in
conjunction with top-down approaches such as wild binary segmentation, which
we will discuss in Section 3.2.

We define the CUSUM transformation $\mathcal{T}:\mathbb{R}^{p\times
n}\to\mathbb{R}^{p\times(n-1)}$ by

$\mathcal{T}(M)_{j,t}=\sqrt{\frac{t(n-t)}{n}}\bigg{(}\frac{1}{n-t}\sum_{r=t+1}^{n}M_{j,r}-\sum_{r=1}^{t}\frac{1}{t}M_{j,r}\bigg{)},$
(5)

and compute the CUSUM matrix $T=\mathcal{T}(X)$. As discussed in Section 2,
our general strategy is to use the matrix $T$ to estimate a projection
direction that is well-aligned with the direction of change, and then project
the data along this direction to estimate the change-point location from the
univariated projected series. More precisely, we would like to solve for

$\hat{v}\in\operatorname*{arg\,max}_{u\in\mathbb{S}^{p-1},\|\phi(u)\|_{0}\leqslant
s}\|u^{\top}T\|_{2},$ (6)

where, $\mathbb{S}^{p-1}=\\{x\in\mathbbm{R}^{p}:\|x\|_{2}=1\\}$. However, the
above optimisation problem is non-convex due to the group-sparsity constraint.
Consequently, we perform the following convex relaxation of the above problem.
We first note that the set of optimisers of (6) is equal to the set of leading
left singular vectors of

$\operatorname*{arg\,max}_{\begin{subarray}{c}M\in\mathbb{R}^{p\times(n-1)}:\|M\|_{*}=1,\mathrm{rank}(M)=1\\\
\sum_{g\in[G]}\mathbbm{1}_{\\{\|M_{\mathcal{J}_{g}}\|_{\mathrm{F}}\neq
0\\}\leqslant s}\end{subarray}}\langle M,T\rangle,$

We relax the above matrix-variate optimisation problem by dropping the
combinatorial rank constraint, and replacing the nuclear norm constraint set
by the larger Frobenius norm set of
$\mathcal{S}=\\{M\in\mathbb{R}^{p\times(n-1)}:\|M\|_{F}\leqslant 1\\}$. The
constraint that $M$ has at most $s$ groups of non-zero rows can be written as
an $\ell_{0}$ constraint on the vector of Frobenius norms of such submatrices,
i.e. $\|(\|M_{\mathcal{J}_{g}}\|_{F}:g\in\\{1,\ldots,G\\})\|_{0}\leqslant s$.
Motivated by the group lasso penalty (Yuan and Lin, 2006), we replace this
group sparsity constraint with a _group norm_ penalty, where the group norm
for a matrix $M\in\mathbb{R}^{p\times(n-1)}$ is defined as

$\|M\|_{\mathrm{grp}}=\sum_{g=1}^{G}p_{g}^{1/2}\|M_{\mathcal{J}_{g}}\|_{2,1},$

where $\|M_{\mathcal{J}_{g}}\|_{2,1}$ is the sum of column $\ell_{2}$ norms of
the submatrix $M_{\mathcal{J}_{g}}$ and $p_{g}=|\mathcal{J}_{g}|$. Overall, we
obtain the following optimisation problem:

$\hat{M}\in\operatorname*{arg\,max}_{M\in\mathcal{S}}\big{\\{}\langle
T,M\rangle-\lambda\|M\|_{\mathrm{grp}}\bigr{\\}},$ (7)

where $\lambda\in[0,\infty)$ is a regularization parameter.

If the groups are non-overlapping, in the sense that
$\mathcal{J}_{g}\cap\mathcal{J}_{g^{\prime}}=\emptyset$ for all $g\neq
g^{\prime}$, then we see from Proposition 5 that (7) has a closed form
solution

$\hat{M}=\frac{T-R^{*}}{\left\lVert T-R^{*}\right\rVert_{\mathrm{F}}},$ (8)

where
$R_{\mathcal{J}_{g},t}^{*}=T_{\mathcal{J}_{g},t}\min\big{\\{}\frac{\lambda
p_{g}^{1/2}}{\|T_{\mathcal{J}_{g},t}\|_{2}},1\big{\\}}$.

For overlapping groups, (7) can be optimised using Frank–Wolfe algorithm
(Frank and Wolfe, 1956), as described in Algorithm 1. We first compute the
gradient of the objective function which is the step 4 in Algorithm 1. We then
project the $\hat{M}$ back onto $\mathcal{S}$.

After solving the optimization problem, we can obtain the estimated projection
direction $\hat{v}$ by computing the leading left singular vector of
$\hat{M}$. Then, we project the data along $\hat{v}$ to obtain a univariate
series for which existing one-dimensional change-point estimation methods
apply. Specifically, we perform the CUSUM transformation over the projected
data series, and locate the change-point by the maximum absolute value of the
CUSUM vector. The full procedure is described in Algorithm 2.

Input: $T\in\mathbb{R}^{p\times(n-1)}$, grouping
$(\mathcal{J}_{g})_{g\in[G]}$, $\lambda>0$ and $\epsilon>0$.

1 Initialise $\hat{M}^{[0]}=T/\|T\|_{\mathrm{F}}$ and $t=0$.

2 repeat

3 $t\leftarrow t+1$

4 Compute
$G^{[t]}=(G_{1}^{[t]},\ldots,G_{p}^{[t]})^{\top}\in\mathbb{R}^{p\times(n-1)}$
such that  $G_{j,t}^{[t]}\leftarrow
T_{j,t}-\sum_{g:j\in\mathcal{J}_{g}}\lambda_{g}\frac{M_{j,t}^{[t-1]}}{\|M_{\mathcal{J}_{g,t}}^{[t-1]}\|_{\mathrm{F}}},$
where $\lambda_{g}=p_{g}^{1/2}\lambda$

5 if $G^{[t]}=0$ then break

6 Compute
$\tilde{M}^{[t]}=\frac{t}{t+2}M^{[t-1]}+\frac{2}{t+2}\frac{G^{[t]}}{\|G^{[t]}\|_{\mathrm{F}}},$

7 Normalise
$\hat{M}^{[t]}\leftarrow\tilde{M}^{[t]}/\|\tilde{M}^{[t]}\|_{\mathrm{F}}$

8until  _$\|\hat{M}^{[t+1]}-\hat{M}^{[t]}\|_{\mathrm{F}}\leqslant\epsilon$_ ;

Output: $\hat{M}^{[t]}$

Algorithm 1 Frank–Wolfe algorithm for optimising (7)

Input: $X\in\mathbb{R}^{p\times n}$, $(\mathcal{J}_{g})_{g\in[G]}$, and
$\lambda>0$

1 Compute $T\leftarrow\mathcal{T}(X)$ as in (5).

2 Solve
$\hat{M}\in\operatorname*{arg\,max}_{M\in\mathcal{S}}\bigl{\\{}\langle
T,M\rangle-\lambda\|M\|_{\mathrm{grp}}\bigr{\\}}$ using either the closed-form
solution in (8) if groups are non-overlapping, or Algorithm 1.

3 Let $\hat{v}$ be the leading left singular vector of $\hat{M}$.

4 Estimate $z$ by $\hat{z}=\operatorname*{arg\,max}_{1\leqslant t\leqslant
n-1}|\hat{v}^{\top}T_{t}|$, where $T_{t}$ is the $t$th column of $T$.

Output: $\hat{z}$, $\bar{T}_{\max}=\hat{v}^{\top}T_{z}$

Algorithm 2 Single change-point estimation procedure for data with group
structure

###  3.2 Multiple change-point estimation

When the data matrix possess multiple change-points, we may combine Algorithm
2 with a top-down approach, such as the wild binary segmentation (Fryzlewicz,
2014), to recursively identify all the change-points. Specifically, we start
by drawing a large number of random intervals
$[s_{1},e_{1}],\ldots,[s_{Q},e_{Q}]$ and apply Algorithm 2 to the data matrix
$X$ restricted to each of these time intervals to obtain $Q$ candidate change-
point locations. We then aggregate $Q$ candidate change-point locations to
choose the one with the maximum projected CUSUM statistics. If the value of
the CUSUM statistic at the best candidate location is above a threshold $\xi$,
we will admit this candidate location as a change-point and repeat the above
process on the data submatrix to the left and right of this change-point. The
pseudocode for the full procedure is given in Algorithm 3.

1

Input: $X\in\mathbb{R}^{p\times n}$, $(\mathcal{J}_{g})_{g\in[G]}$,
$\lambda>0$, $\xi>0$, $Q\in\mathbb{N}$

2 Set $\hat{Z}\leftarrow\emptyset$

3 Draw $Q$ pairs of integers $(s_{1},e_{1}),\ldots,(s_{Q},e_{Q})$ uniformly at
random from the set $\\{(\ell,r)\in\mathbb{Z}^{2}:0\leqslant\ell<r\leqslant
n\\}$

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

5 Set $Q_{s,e}=\\{q:s\leqslant s_{q}<e_{q}\leqslant e\\}$

6 for  _$q\in Q_{s,e}$_ do

7 $(\hat{z}^{[q]},\bar{T}^{[q]}_{\max})\leftarrow$ output of Algorithm 2 with
inputs $(X_{j,t})_{j\in[p],t\in(s,e]}$ and $\lambda$

8 Find $q_{0}\in\operatorname*{arg\,max}_{q\in Q_{s,e}}\bar{T}^{[q]}_{\max}$
and set $b\leftarrow s_{q_{0}}+\hat{z}^{[q_{0}]}$

9 if  _$\bar{T}^{[q_{0}]}_{\max}\geqslant\xi$_ then

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

11 Run recursively $\textbf{wbs}(s,b)$ and $\textbf{wbs}(b,e)$

12

Output: $\hat{Z}$

Algorithm 3 Multiple change-point estimation procedure

##  4 Theoretical guarantees

In this section, we provide theoretical guarantees to the performance of the
groupInspect algorithm. As we have noted in Section 2, a key to the successful
change-point estimation in the current problem is a good estimator of the
oracle projection direction $v=\theta/\|\theta\|_{2}$.

The following theorem controls the sine angle risk of the estimated projection
direction $\hat{v}$ in Step 3 of Algorithm 2. We define
$\mathcal{P}_{n,p}(s,k,\tau,\vartheta,\sigma^{2},(\mathcal{J}_{g})_{g\in[G]})$
to be the set of data distributions satisfying (1), (2), (3) and (4). For any
$P\in\mathcal{P}$, we write $v(P)=\theta/\|\theta\|_{2}$ where $\theta$ is the
difference between post-change and pre-change means.

######  Theorem 1.

For a given grouping $(\mathcal{J}_{g})_{g\in[G]}$, let
$p_{*}=\min_{g\in[G]}|\mathcal{J}_{g}|$ and suppose further that there exists
a universal constant $C_{1}>0$, such that
$\max_{j\in[p]}|\\{g:j\in\mathcal{J}_{g}\\}|\leqslant C_{1}$. Let $X\sim
P\in\mathcal{P}_{n,p}(s,k,\tau,\vartheta,\sigma^{2},(\mathcal{J}_{g})_{g\in[G]})$
be a $p\times n$ data matrix, let $\theta$ be the vector of change and let
$\hat{v}$ be as in Step 3 of Algorithm 2 with input $X$,
$(\mathcal{J}_{g})_{g\in[G]}$ and
$\lambda\geqslant\sigma(1+\sqrt{4\log(Gn)/p_{*}})$. Then there exists $C>0$,
depending only on $C_{1}$, such that

$\sup_{P\in\mathcal{P}_{n,p}(s,k,\tau,\vartheta,\sigma^{2},(\mathcal{J}_{g})_{g\in[G]})}\mathbb{P}_{P}\bigg{\\{}\sin\angle{(\hat{v},v)}\leqslant\frac{C\lambda
k^{1/2}}{n^{1/2}\tau\vartheta}\bigg{\\}}\leqslant\frac{1}{nG}.$ (9)

We remark that the condition
$\max_{j\in[p]}|\\{g:j\in\mathcal{J}_{g}\\}|\leqslant C_{1}$ is to control the
extent of overlapping between different groups. Specifically, it requires that
each coordinate can belong to at most $C_{1}$ groups. In the special case when
all groups $\mathcal{J}_{g}$ are disjoint, which is often true in practical
applications, then it suffices to take $C_{1}=1$.

We note that, when $\lambda=\sigma(1+\sqrt{4\log(Gn)/p_{*}})$, with high
probability, the sine angle loss in (9) has an upper bound that is
proportional to $\sigma k^{1/2}n^{-1/2}\tau^{-1}\vartheta^{-1}$, similar to
what has been previously observed in Wang and Samworth (2018, Proposition 1).
However, Theorem 1 reveals an interesting interaction between the $\ell_{0}$
sparsity $k$ and the group sparsity $s$ when all groups are of comparable
size. Specifically, for $\lambda=\sigma(1+\sqrt{4\log(Gn)/p_{*}})$ and
assuming that $\max_{g\in[G]}p_{g}\lesssim p_{*}$, then we can simplify (9) to
obtain that

$\mathbb{E}\\{\sin\angle{(\hat{v},v)}\\}\lesssim\sqrt{\frac{\sigma^{2}\\{k+s\log(Gn)\\}}{n\tau^{2}\vartheta^{2}}}.$

In other words, the risk upper bound undergoes a phase transition as the
number of coordinates per group increases above a $\log(Gn)$ level. Similar
phase transitions have been previously observed in the context of high-
dimensional linear model where the regression coefficients satisfy a group
sparsity assumption (see, e.g. Cai et al., 2019, Theorem 3).

We now turn our attention to a minimax lower bound of the estimation risk of
the oracle projection direction. Theorem 2 below shows that the phase
transition observed in Theorem 1 is not due to the specific proof techniques
employed but rather an intrinsic feature of the problem.

######  Theorem 2.

Suppose $s>0$, $k>0$ and a grouping $(\mathcal{J}_{g})_{g\in[G]}$ satisfy that
$\mathcal{J}_{g}\cap\mathcal{J}_{g^{\prime}}=\emptyset$ for all $g\neq
g^{\prime}$ , $\min\\{k,(s-1)\log(G/s)\\}\geqslant 20$, and
$\sum_{r=1}^{s}p_{(G-r+1)}\geqslant k/2$, where $p_{(1)}\leqslant
p_{(2)}\leqslant\cdots\leqslant p_{(G)}$ are order statistics of
$p_{1},\ldots,p_{G}$. Then for some universal constant $c>0$, we have

$\inf_{\tilde{v}}\sup_{P\in\mathcal{P}_{n,p}(s,k,\tau,\vartheta,\sigma^{2},(\mathcal{J}_{g})_{g\in[G]})}\mathbb{E}_{P}L(\tilde{v}(X),v(P))\geqslant
c\sqrt{\frac{\sigma^{2}\\{k+s\log(G/s)\\}}{n\tau\vartheta^{2}}},$

where the infimum is taken over the set of all measurable functions
$\tilde{v}$ of the data $X$.

The condition that $\sum_{r=1}^{s}p_{(G-r+1)}\geqslant k/2$ is to ensure that
the upper bound $k$ on the $\ell_{0}$-sparsity is not too loose in the sense
that $k$ is not too much larger than the cardinality of the union of the
largest $s$ groups. If we assume that $\log(G/s)\asymp\log(n)$, $\tau\asymp 1$
and $\max_{g\in[G]}p_{g}\lesssim p_{*}$, then the lower bound in Theorem 2
matches the upper bound of Theorem 1 up to universal constants, when all
groups are non-overlapping.

After obtaining guarantees on the quality of the projection direction
estimator, we now provide theoretical guarantees of the overall change-point
procedure. We note that the projection direction estimator $\hat{v}$ is
dependent on the CUSUM panel $T$. While this dependence is observed to be very
weak in practice, it creates difficulties in analysing the projected CUSUM
series $\hat{v}^{\top}T$ in Step 4 of Algorithm 2. As such, for theoretical
convenience, we will instead analyse a sample-splitting version of the
algorithm. Specifically, we split the data into $X^{(1)}$ and $X^{(2)}$,
consisting of odd and even time points respectively, as described in Algorithm
4. We use $X^{(1)}$ to estimate the projected direction $\hat{v}^{(1)}$ and
then project $X^{(2)}$ along this direction to locate the change-point.
Theorem 3 below provides a performance guarantee for the estimated location of
the change-point of this sample-splitting version of our procedure.

Input: $X\in\mathbb{R}^{p\times n}$ and $\lambda>0$

1 Define $X^{(1)}$ as $X^{(1)}_{j,t}=X_{j,2t-1}$ and $X^{(2)}$ as
$X^{(2)}_{j,t}=X_{j,2t}$.

2 Compute $T^{(1)}\leftarrow\mathcal{T}(X^{(1)})$ and
$T^{(2)}\leftarrow\mathcal{T}(X^{(2)})$ as in (5).

3 Solve
$\hat{M}^{(1)}\in\operatorname*{arg\,max}_{M\in\mathcal{S}}\bigl{\\{}\langle
T^{(1)},M\rangle-\lambda\|M\|_{\mathrm{grp}}\bigr{\\}}$ using either the
closed-form solution in (8) if groups are non-overlapping, or Algorithm 1.

4 Let $\hat{v}$ be the leading left singular vector of $\hat{M}^{(1)}$.

5 Estimate $z$ by $\hat{z}=2\operatorname*{arg\,max}_{1\leqslant t\leqslant
n_{1}-1}|(\hat{v}^{(1)})^{\top}T_{t}^{(2)}|$, where $T_{t}^{(2)}$ is the $t$th
column of $T^{(2)}$.

Output: $\hat{z}$

Algorithm 4 Change-point estimation procedure: sample splitting version

######  Theorem 3.

Given data matrix $X\sim
P\in\mathcal{P}_{n,p}(s,k,\tau,\vartheta,\sigma^{2},(\mathcal{J}_{g})_{g\in[G]})$,
let $\hat{z}$ be the output from the Algorithm 4 with input $X$ and
$\lambda=\sigma(1+\sqrt{p_{*}^{-1}4\log(nG)})$. There exist universal
constants $C$, $C^{\prime}>0$ such that, if $n\geqslant 12$ is even, $z$ is
even, and

$\frac{C\sigma\sqrt{k}}{\vartheta\tau\sqrt{n}}\bigg{(}\frac{1+\sqrt{4\log(Gn)}}{p_{*}}\bigg{)}\leqslant
1,$

then,

$\mathbb{P}\bigg{\\{}\frac{1}{n}|\hat{z}-z|\leqslant\frac{C^{\prime}\sigma^{2}(1+\sqrt{4\log(n)})}{n\vartheta^{2}}\bigg{\\}}\geqslant
1-\frac{20\log n}{n}.$

##  5 Numerical studies

In this section, we provide some simulation results to demonstrate the
empirical performance of the groupInspect method. In all our numerical
studies, unless otherwise specified, we will assume that data are generated
according to (1), (2), (3) and (4), with $\sigma=1$. In all simulations, we do
not assume that $\sigma$ is known, or even equal across rows. Instead, we
estimate the variance in each row using the mean absolute deviation of
successive differences of the observations. We then standardise the data by
the estimated row standard deviation. The groupInspect procedure is then
applied to the standardised data with $\sigma=1$.

###  5.1 Theory validation

We first show that the practical performance of the groupInspect procedure is
well captured by the theoretical results in Theorems 1 and 2. There are two
related measures of the signal sparsity in our problem, which are the total
number of coordinates of change $k$ and the total number of groups with a
change $s$. We conduct two sets of simulation experiments fixing one of these
sparsity measures and varying the other. Specifically, for $n=1000$,
$p\in\\{600,1200,2400\\}$ and $\vartheta\in\\{1,2,4,8,16\\}$, we split the $p$
coordinates into disjoint groups of $p_{*}$ coordinates per group, where
$p_{*}$ is allowed to vary over all divisors of $60$. In the first set of
experiments, we fix $k=60$ so that $s=k/p_{*}$ varies with $p_{*}$, whereas in
the second set of experiments, we fix $s=3$ so that $k=sp_{*}$ varies with
$p_{*}$. The vector of change is constructed so that the magnitude of change
is equal across all coordinates of change. We will use the theoretical choice
of tuning parameter $\lambda$ for both sets of experiments here. Figure 1
shows how the $\sin\theta$ loss, averaged over $100$ Monte Carlo repetitions,
varies with $p_{*}$, for different choices of $p$ and $\vartheta$ in both
settings.

In the left panel of Figure 1, where the number of signal coordinates $k$ is
fixed, we see that the average loss decreases as $p_{*}$ increases.
Furthermore, at a log-log scale, and for relatively large signal sizes of
$\vartheta\in\\{4,8,16\\}$, we see the loss curves follow an initial linear
decreasing trend as $p_{*}$ increases before plateauing eventually. This is in
agreement with the two terms contributing to the loss described in Theorem 1.
Specifically, for small $p_{*}$, we expect the second term of (9) to dominate
and the loss decreases at a rate approximately proportional to
$1/\sqrt{p_{*}}$ initially. For large $p_{*}$, we expect the first term of (9)
to dominate and the loss will have minimal dependence on $p_{*}$. In the right
panel of Figure 1, where the number of signal groups $s$ is fixed, the average
loss increases with $p_{*}$, as expected from our theory. It appears that for
$s=3$ studied here, the first term of (9) is dominant and the average loss
increases linearly at the log-log scale with respect to $p_{*}$.

We further remark that in both panels of Figure 1, the average loss for large
$p_{*}$ shows equally spaced separation for the signal size $\vartheta$ in the
dyadic grid $\\{1,2,4,8,16\\}$. This is in good agreement with the $1/\theta$
dependence of expected loss given in Theorem 1. Finally, we note that the
ambient dimension $p$ has minimal effect on the loss curves, for all signal
strengths studied here. Again, this is predicted by our theory as the
dimension $p$ enters the mean loss in (9) only through the
$\log(Gn)=\log(pn/p_{*})$ expression in the second term.

Figure 1: Average loss (over 100 repetitions) of groupInspect for varying
elements per group $p_{*}$, plotted on a log-log scale. Left panel: $k=60$ and
$s=k/p_{*}$. Right panel: $s=3$ and $k=sp_{*}$. Other parameter: $n=1000$.

###  5.2 Practical choice of tuning parameter

The theoretical choice of $\lambda$ turns out to be conservative in practical
use. In this subsection, we will perform numerical simulations to suggest a
suitable practical tuning parameter choice. We fix $n=1000$, $z=400$, $s=3$,
$G\in\\{10,25\\}$. The signal size $\vartheta$ is varied in $\\{1,2,4,8,16\\}$
and $p$ is chosen from $\\{500,1000\\}$. All groups are set to have equal
size. For the choice of tuning parameters, we first form a logarithmic
sequence of values between $0.1$ and $3$ with length $7$ and then times each
value with the theoretically suggested value of $1+\sqrt{4p_{*}^{-1}\log(nG)}$
to form the sequence of the tuning parameter. For each setting of the signal
size, we will run algorithm with all the $\lambda$ values and record the sine
angle loss.

We plot $\sin\theta$ loss against $\lambda$ in Figure 2. The $x$-axis is the
log sequence. In most cases, the loss is minimized when tuning parameter value
is half of the theoretical value. However, when the minimum loss is achieved
by other values of $\lambda$, this lambda value can still achieve the loss
which is close to optimal value. Therefore, we suggest that using
$\lambda=1/2(1+\sqrt{4p_{*}^{-1}\log(nG)})$ in practical is less conservative.

Figure 2: Average loss (over 100 repetitions) of groupInspect for varying
tuning parameter $\lambda$. Left panel: $G=10$. Right panel: $G=25$. Other
parameter: $n=1000$, $s=3$.

###  5.3 Comparison between different methods

Now, we would like to compare our method with other existing change-point
estimation procedures. As groupInspect is a two-stage procedure that first
estimates a projection direction before localising the change-point on the
projected series, we will investigate its performance both in terms of its
accuracy in estimating the projection direction and the quality of the final
change-point location estimator. For the former, we compare the estimated
projection direction from groupInspect with that from the inspect algorithm.
We measure the accuracy in terms of the sine angle loss introduced in Section
2. We use the recommended values for tuning parameters in both methods, i.e.,
$\sqrt{2^{-1}\log\\{p\log n\\}}$ in inspect as in Wang and Samworth (2018) and
$1/2(1+\sqrt{4p_{*}^{-1}\log(nG)})$ for groupInspect as suggested in Section
5.2.

We fix $n=1000$, $p=1000$ and vary $\vartheta$ in $\\{1,2,4,8,16\\}$. We
consider settings with both non-overlapping groups and overlapping groups. For
the non-overlapping setting, we have $G=10$ groups of equal size $p_{*}=100$,
whereas for the overlapping setting, we have $G=19$ groups of size $100$ each,
where neighbouring groups overlap in exactly $50$ coordinates. Both methods
have access to exactly the same data sets and the performance is averaged over
100 Monte Carlo repetitions.

Figure 3 shows the comparison of the average sine angle loss between
groupInspect and inspect over all signal sizes on a logarithmic scale, in both
the non-overlapping and overlapping settings. In both cases, groupInspect
outperforms the inspect algorithm. From the left panel, we can see that the
estimation accuracy of the projection direction using groupInspect is
substantially better even when the signal is small.

We now turn our attention to the overall change-point localisation accuracy of
the groupInspect procedure. To this end, we compare the mean absolute
deviation of various high-dimensional change-point procedures over 100 Monte
Carlo repetitions using the same data sets. In addition to inspect, we also
compare against the $\ell_{2}$ aggregation procedures of Horváth and Hušková
(2012), the $\ell_{\infty}$ aggregation procedure of Jirak (2015) and the
double CUSUM procedure of Cho (2016). We set $n=1000$,
$p\in\\{500,1000,2000\\}$, $\vartheta\in\\{0.25,0.5,1,2,4\\}$. The simulation
results are presented in Table 1. For simplicity, we have only shown the
results for 10 equal-sized non-overlapping groups here, but qualitatively
similar results were obtained in other settings as well. We see that
groupInspect is very competitive over a wide range of dimensions and signal-
to-noise ratio settings, though the benefit of using the group sparsity
structure via groupInspect is most apparent in low signal-to-noise ratio
settings where the change-point estimation problem is more difficult.

|
---|---
Figure 3: Average loss (over 100 repetitions) comparison between groupInspect and Inspect. Left panel: non-overlap setting. Right panel: overlap setting $p$ | $\vartheta$ | groupInspect | inspect |  $\ell_{2}$-aggregate |  $\ell_{\infty}$-aggregate | double cusum
---|---|---|---|---|---|---
$500$ | $0.25$ | $\mathbf{127}$ | $143$ | $336$ | $337$ | $347$
$500$ | $0.5$ | $\mathbf{59.8}$ | $93.4$ | $231$ | $305$ | $262$
$500$ | $1$ | $\mathbf{3.83}$ | $8.83$ | $9.84$ | $94.2$ | $40.6$
$500$ | $2$ | $\mathbf{0.670}$ | $0.982$ | $0.875$ | $16.0$ | $4.16$
$500$ | $4$ | $0.045$ | $\mathbf{0.018}$ | $0.045$ | $4.04$ | $0.179$
$1000$ | $0.25$ | $\mathbf{108}$ | $138$ | $347$ | $348$ | $363$
$1000$ | $0.5$ | $\mathbf{81.8}$ | $107$ | $269$ | $326$ | $297$
$1000$ | $1$ | $\mathbf{15.6}$ | $34.6$ | $22.1$ | $204$ | $57.9$
$1000$ | $2$ | $\mathbf{0.920}$ | $1.51$ | $0.973$ | $28.3$ | $3.91$
$1000$ | $4$ | $\mathbf{0.081}$ | $0.117$ | $0.099$ | $6.70$ | $0.387$
$2000$ | $0.25$ | $\mathbf{101}$ | $139$ | $358$ | $365$ | $364$
$2000$ | $0.5$ | $\mathbf{91.2}$ | $127$ | $305$ | $353$ | $321$
$2000$ | $1$ | $\mathbf{36.3}$ | $58.1$ | $71.6$ | $305$ | $127$
$2000$ | $2$ | $\mathbf{1.88}$ | $2.76$ | $2.32$ | $52.6$ | $6.27$
$2000$ | $4$ | $\mathbf{0.134}$ | $0.161$ | $\mathbf{0.134}$ | $7.97$ | ${0.696}$
Table 1: Average mean absolute deviation (over 100 repetitions) comparison
between different methods. Other parameters used: $n=1000$ with $G=10$

###  5.4 Multiple change-points simulation

The numerical studies so far have focused mainly on the single change-point
estimation problem. In this subsection, we investigate the empirical
performance of groupInspect in multiple change-point estimation tasks. We will
compare its performance as implemented in Algorithm 3 to that of the inspect
algorithms for estimating multiple change-points under different settings. We
choose $n=1200$, $p\in\\{500,1000\\}$, $s\in\\{3,10\\}$, $G\in\\{50,100\\}$.
Each data series contains three true change-points located at $300$, $600$ and
$900$ with the $\ell_{2}$ norm of the change equal to $\vartheta$,
$1.5\vartheta$ and $2\vartheta$ respectively. We vary $\vartheta$ in
$\\{0.6,0.8,1,1.2,1.4\\}$. For simplicity, we further assume that the same $s$
coordinates undergo change in all three change-points and that all groups have
10 elements. We use the $\lambda$ tuning parameter choice suggested in Section
5.2 for the groupInspect method and that suggested in Wang and Samworth (2018)
for the inspect algorithm. For the thresholding parameter $\xi$ of the wild
binary segmentation recursion used in both groupInspect and inspect, we choose
via Monte Carlo simulation. More precisely, we randomly generate 1000 data
sets from the null model with no change-points and take the maximum absolute
CUSUM statistics from Algorithm 3 and Wang and Samworth (2018, Algorithm 4) as
$\xi_{g}$ and $\xi_{i}$ respectively. We compare the performance of two
algorithms using the Adjusted Rand index (ARI) of the estimated segmentation
against the truth (Rand, 1971; Hubert and Arabie, 1985).

From Figure 4, we see that the groupInspect algorithm generally performs much
better than the inspect algorithm in the multiple change-point localisation
tasks. The advantage of groupInspect is more pronounced when the signal is
sparser and when the dimension of the data is higher.

|
---|---
Figure 4: Average ARI comparsion between groupInspect and inspect. Left panel:
$p=500,G=50$. Right panel: $p=1000,G=100$.

To further visualise the output of the two procedures, we plot the estimated
change-point locations for one specific setting ($s=3$ and $\vartheta=1$) of
each of the two panels in Figure 4. The resulting histograms in Figure 5 shows
that when $p=500$, groupInspect was better at picking out all three change-
points with higher accuracies. When $p=1000$, inspect was only able to pick
out the change at $t=600$ in most of the trials, whereas groupInspect was
still able to identify even the weakest change signal at $t=300$ in a
substantial fraction of all trials.

|
---|---
|
Figure 5: Histograms of estimated locations by groupInspect and inspect under
two settings when $P=500,G=50$ and $p=1000,G=100$. Other parameter used:
$s=3$, $\vartheta=1$ are fixed in both settings.

###  5.5 Real data analysis

In this section, we apply groupInspect to a stock price data. The data
consists of the logarithmic daily returns (computed from the adjusted closing
prices) of S&P 500 stocks during the period 1 January 2007 to 31 December
2011. Since not all companies remained in the S&P 500 list and some companies
have missing data at a few time points, we eventually selected 256 companies
which have continuously traded throughout this this period to construct a
multivariate time series of dimension $p=256$ and length $n=1259$. We then
divide the $256$ companies into $G=11$ non-overlapping groups according to
respective Global Industry Classification Standard sector memberships. We then
rescale rows of the data matrix by their estimated standard deviation as in
Section 5. We use the same procedure in Section 5.4 to choose thresholding
parameter $\xi$.

The groupInspect algorithm identifies the following change points $t=147$,
$148$, $298$, $386$, $427$, $441$, $448$, $460$, $477$, $522$, $524$, $549$,
$559$, $1158$, $1189$, as illustrated in Figure 6. We see a large number of
changes being identified in the period between September and October 2008,
which corresponds to the period when the financial crisis reaches a climax,
and when the stock market is most volatile.

Figure 6: Estimated change point locations (red dashed lines) by groupInspect
applied to the stock return data. For ease of illustration, we have plotted
the $\ell_{2}$ norm of the returns of all stocks within each of the 11 groups
over time.

##  6 Proofs of main results

In this section, we will give the proof of our results in section 4.

###  6.1 Proof of Theorem 1

###### Proof.

By Wang and Samworth (2018, equation (9)), we can explicitly write the matrix
$A=(A_{j,t})_{j\in[p],t\in[n-1]}$ by

$A_{j,t}=\begin{cases}\sqrt{\frac{t}{n(n-t)}}(n-z)\theta_{j}&\text{if
$1\leqslant t\leqslant z$},\\\ \sqrt{\frac{n-t}{nt}}z\theta_{j}&\text{ if
$z<t\leqslant n-1$}.\end{cases}$

Thus, we have

$A=\theta\gamma^{\top},$

where
$\gamma=\frac{1}{\sqrt{n}}\bigr{(}\sqrt{\frac{1}{n-1}}(n-z),\sqrt{\frac{2}{n-2}}(n-z),\cdots\sqrt{z(n-z)},\sqrt{\frac{n-z-1}{z+1}}z,\cdots,\sqrt{\frac{1}{n-1}}z\bigr{)}^{\top}$.
In particular, by By Wang and Samworth (2018, Lemma 3), $A$ is a rank 1 matrix
with $\|A\|_{\mathrm{op}}=\|\theta\|_{2}\|\gamma\|_{2}\geqslant
n\tau\vartheta/4$. By Proposition 9 with $\delta=(nG)^{-2}$, we have

$\mathbb{P}(\|T-A\|_{\mathrm{grp}*}>\lambda)<\frac{1}{Gn}.$

By Proposition 8, on the event
$\\{\|T-A\|_{\mathrm{grp}_{*}}\leqslant\lambda\\}$, we have

$\max\bigl{\\{}\sin\angle(v,\hat{v}),\sin\angle(u,\hat{u})\bigr{\\}}\leqslant\frac{32\lambda(C_{1}k)^{1/2}}{n^{1/2}\tau\vartheta},$

as desired. ∎

###  6.2 Proof of Theorem 2

###### Proof.

We will use two different constructions to derive separate lower bounds of
order $\sqrt{\sigma^{2}s\log(G/s)/(n\tau\vartheta^{2})}$ and
$\sqrt{\sigma^{2}k/(n\tau\vartheta^{2})}$ respectively. Without loss of
generality, we may assume that $z<n/2$.

For the first bound, let $s_{0}=s-1$, $G_{0}=G-1$, then in $\mathbb{R}^{G}$.
By the Gilbert–Varshamov lemma as stated in Massart (2007, Lemma 4.10)
(applied with $\alpha=3/4$ and $\beta=1/3$), we can construct a set
$\mathcal{U}_{0}$ of $s_{0}$-sparse vectors in $\\{0,1\\}^{G_{0}}$, with
cardinality at least $({G_{0}}/{s_{0}})^{s_{0}/5}$, such that the pairwise
Hamming distance between any pair of vectors in $\mathcal{U}_{0}$ is at least
${s_{0}}/{2}$. Let $\epsilon\in(0,1)$ to be chosen later, we can define a set

$\mathcal{U}=\bigg{\\{}\biggl{(}\begin{array}[]{c}\sqrt{1-\epsilon^{2}}\\\
s_{0}^{-1/2}\epsilon u_{0}\\\
\end{array}\bigg{)}:u_{0}\in\mathcal{U}_{0}\bigg{\\}}\subseteq\mathbb{S}^{G-1}.$

We remark that for any pair of distinct $u,u^{\prime}\in\mathcal{U}$, we have
by construction that
$\epsilon/\sqrt{2}\leqslant\|u^{\prime}-u\|_{2}\leqslant\epsilon$. We then
define a map $\psi:\mathbb{R}^{G}\to\mathbb{R}^{p}$ such that for any
$u\in\mathcal{U}$ and $j\in\mathcal{J}_{g}$, we have
$\psi(u)_{j}=u_{g}p_{g}^{-1/2}$. Finally, let
$\mathcal{V}=\\{\psi(u):u\in\mathcal{U}\\}$. We note that
$\|\psi(u^{\prime})-\psi(u)\|_{2}=\|u^{\prime}-u\|_{2}$. Therefore, for
distinct $v,v^{\prime}\in\mathcal{V}$, we have

$L(v^{\prime},v)=\sqrt{1-(v^{\top}v^{\prime})^{2}}=\frac{\|v^{\prime}-v\|_{2}}{\sqrt{2}}\geqslant\frac{\epsilon}{2}.$
(10)

Now, for each $v\in\mathcal{V}$, we can define a distribution
$P_{v}\in\mathcal{P}_{n,p}(s,k,\tau,\vartheta,\sigma^{2},(\mathcal{J}_{g})_{g\in[G]})$,
such that the pre-change mean is $-\vartheta v$ and the post-change mean is
$0$ (we check that $P_{v}$ indeed satisfies the conditions of
$\mathcal{P}_{n,p}(s,k,\tau,\vartheta,\sigma^{2},(\mathcal{J}_{g})_{g\in[G]})$).
Then for any distinct $v,v^{\prime}\in\mathcal{V}$, we have

$D(P_{v}\|P_{v^{\prime}})=zD(N_{p}(-v\vartheta,\sigma^{2}I_{p})\|N_{p}(-v^{\prime}\vartheta,\sigma^{2}I_{p}))=\frac{z\vartheta^{2}}{2}\|v-v^{\prime}\|_{2}^{2}\leqslant\frac{z\vartheta^{2}\epsilon^{2}}{2\sigma^{2}}.$
(11)

By (10) and (11), we can apply Fano’s lemma (Yu, 1997, Lemma 3) to obtain that

$\displaystyle\inf_{\tilde{v}}\sup_{P\in\mathcal{P}_{n,p}(s,k,\tau,\vartheta,\sigma^{2},(\mathcal{J}_{g})_{g\in[G]})}\mathbb{E}_{P}L(\tilde{v}(X),v(P))$
$\displaystyle\geqslant\inf_{\tilde{v}}\sup_{v\in\mathcal{V}}\mathbb{E}_{P_{v}}L(\tilde{v}(X),v)$
$\displaystyle\geqslant\frac{\epsilon}{4}\biggl{\\{}1-\frac{z\vartheta^{2}\epsilon^{2}/2\sigma^{2}+\log
2}{(s_{0}/5)\log(G_{0}/s_{0})}\biggr{\\}}.$

By the condition $(s-1)\log(G/s)\geqslant 20$ in the theorem, we have
$(s_{0}/5)\log(G_{0}/s_{0})\geqslant 2\log 2$. Moreover, the choice of

$\epsilon=\sqrt{\frac{\sigma^{2}s_{0}\log(G_{0}/s_{0})}{10z\vartheta^{2}}}$

ensures that $(s_{0}/5)\log(G_{0}/s_{0})\geqslant
2z\vartheta^{2}\epsilon^{2}/\sigma^{2}$. Therefore,

$\inf_{\tilde{v}}\sup_{P\in\mathcal{P}_{n,p}(s,k,\tau,\vartheta,\sigma^{2},(\mathcal{J}_{g})_{g\in[G]})}\mathbb{E}_{P}L(\tilde{v}(X),v(P))\geqslant\frac{\epsilon}{16}\geqslant\frac{1}{72}\sqrt{\frac{\sigma^{2}s\log(G/s)}{z\vartheta^{2}}}.$
(12)

For the second lower bound, let $g_{1},\ldots,g_{s}$ be the indices of the $s$
groups with largest cardinalities. By the given condition of the Theorem, we
have that
$\tilde{k}=\sum_{r=1}^{s}p_{g_{r}}=\sum_{r=1}^{s}p_{(G-r+1)}\geqslant k/2$.
Let $S=\cup_{r=1}^{s}\mathcal{J}_{g_{r}}$, so $|S|=\tilde{k}$. By Massart
(2007, Lemma 4.7), we can construct a subset $\mathcal{V}_{0}$ of
$\\{-1,1\\}^{\tilde{k}_{0}}$ of cardinality at least $e^{\tilde{k}/8}$, such
that any two points in the set are separated in Hamming distance by at least
$\tilde{k}/4$. Construct

$\mathcal{V}=\biggl{\\{}v:v_{S}=\begin{pmatrix}\sqrt{1-\epsilon^{2}}\\\
\tilde{k}_{0}^{-1/2}\epsilon v_{0}\end{pmatrix}\text{ for some
$v_{0}\in\mathcal{V}_{0}$ and $v_{S^{\mathrm{c}}}=0$}\biggr{\\}}.$

Therefore, for distinct $v,v^{\prime}\in\mathcal{V}$, we have
$\epsilon\leqslant\|v^{\prime}-v\|_{2}\leqslant 2\epsilon$,then,

$L(v^{\prime},v)=\sqrt{1-(v^{\top}v^{\prime})^{2}}=\frac{\|v^{\prime}-v\|_{2}}{\sqrt{2}}\geqslant\frac{\epsilon}{\sqrt{2}}.$

Following the same derivation as in (11), we have that

$D(P_{v}\|P_{v^{\prime}})=zD(N_{p}(-v\vartheta,\sigma^{2}I_{p})\|N_{p}(-v^{\prime}\vartheta,\sigma^{2}I_{p}))=\frac{z\vartheta^{2}}{2\sigma^{2}}\|v-v^{\prime}\|_{2}^{2}\leqslant
2z\vartheta^{2}\epsilon^{2}/\sigma^{2}.$

Again, we can use Fano’s lemma (Yu, 1997, Lemma 3) to obtain that

$\displaystyle\inf_{\tilde{v}}\sup_{v\in\mathcal{V}}\mathbb{E}_{P_{v}}L(\tilde{v}(X),v)\geqslant\frac{\epsilon}{\sqrt{2}}\biggl{\\{}1-\frac{2z\vartheta^{2}\epsilon^{2}/\sigma^{2}+\log
2}{\tilde{k}/8}\biggr{\\}}\geqslant\frac{\epsilon}{\sqrt{2}}\biggl{\\{}1-\frac{2z\vartheta^{2}\epsilon^{2}/\sigma^{2}+\log
2}{k/16}\biggr{\\}}.$

Now, choose $\epsilon=\sigma k^{1/2}z^{-1/2}\vartheta^{-1}/4\sqrt{6}$. Since
$k\geqslant 20$, we have $k/16\geqslant 9\log(2)/5$, so that

$\displaystyle\inf_{\tilde{v}}\sup_{P\in\mathcal{P}_{n,p}(s,k,\tau,\vartheta,\sigma^{2},(\mathcal{J}_{g})_{g\in[G]})}\mathbb{E}_{P}L(\tilde{v}(X),v(P))$
$\displaystyle\geqslant\inf_{\tilde{v}}\sup_{v\in\mathcal{V}}\mathbb{E}_{P_{v}}L(\tilde{v}(X),v)$
$\displaystyle\geqslant\frac{\epsilon}{9\sqrt{2}}\geqslant\frac{1}{72\sqrt{3}}\sqrt{\frac{\sigma^{2}k}{z\theta^{2}}}.$
(13)

The desired result follows by combining (12) with (6.2), and noting that
$z\geqslant n\tau$. ∎

###  6.3 Proof of Theorem 3

###### Proof.

Recall the definition of $X^{(2)}$ and let $T^{(2)}=\mathcal{T}(X^{(2)})$.
Define similarly
$\mathbf{\mu}^{(2)}=(\mu_{1}^{(2)},\ldots,\mu_{n_{1}}^{(2)})\in\mathbb{R}^{p\times
n_{1}}$ and a random $W^{(2)}=(W_{1}^{(2)},\ldots,W_{n_{1}}^{(2)})$ taking
values in $\mathbb{R}^{p\times n_{1}}$ by $\mu_{t}^{(2)}=\mu_{2t}$ and
$W_{t}^{(2)}=W_{2t}$. Now, let $A^{(2)}=\mathcal{T}(\mathbf{\mu}^{(2)})$ and
$E^{(2)}=\mathcal{T}(W^{(2)})$. We also write
$\bar{X}=(\hat{v}^{(1)})^{\top}X^{(2)}$,$\bar{\mu}=(\hat{v}^{(1)})^{\top}\mathbf{\mu}^{(2)}$,$\bar{W}=(\hat{v}^{(1)})^{\top}W^{(2)}$,$\bar{A}=(\hat{v}^{(1)})^{\top}A^{(2)}$,
$\bar{E}=(\hat{v}^{(1)})^{\top}E^{(2)}$ and
$\bar{T}=(\hat{v}^{(1)})^{\top}T^{(2)}$ for the one-dimensional projected
images. Note that by linearity, we have $\bar{T}=\mathcal{T}(\bar{X})$,
$\bar{A}=\mathcal{T}(\bar{\mu})$ and $\bar{E}=\mathcal{T}(\bar{W})$,

Now, conditional on $\hat{v}^{(1)}$, the random variables
$\bar{X}_{1},\ldots,\bar{X}_{n_{1}}$ are independent with

$\bar{X}_{t}\mid\hat{v}^{(1)}\sim N(\bar{\mu}_{t},\sigma^{2})$

and the row vector $\bar{\mu}$ undergoes a single change at $z^{(2)}=z/2$ with
magnitude of change

$\bar{\theta}=\bar{\mu}_{z^{(2)}+1}-\bar{\mu}_{z^{(2)}}=\hat{v}^{(1)\top}\theta.$

Finally, let $\hat{z}^{(2)}\in\operatorname*{arg\,max}_{1\leqslant t\leqslant
n_{1}-1}|\bar{T_{t}}|$, so the first component of the output of the algorithm
is $\hat{z}=2\hat{z}^{(2)}$. Consider the set

$\Upsilon=\\{u\in\mathbb{S}^{p-1}:\angle(u,v)\leqslant\pi/6\\}.$

By Condition (3) and Theorem 1, we have that

$\mathbb{P}(\hat{v}^{(1)}\in\Upsilon)\geqslant 1-\frac{1}{n_{1}G}.$ (14)

Moreover, on the event $\\{\hat{v}^{(1)}\in\Upsilon\\}$, we have that
$\bar{\theta}\geqslant\sqrt{3}\vartheta/2$. Set
$\lambda_{1}=\sigma(1+\sqrt{4\log n})$, we have by Proposition 9 that

$\mathbb{P}(\|\bar{E}\|_{\infty}\geqslant\lambda_{1})=\mathbb{P}(\|\bar{E}\|_{\mathrm{grp}*}\geqslant\lambda)\leqslant\frac{1}{n_{1}}.$
(15)

Since $\bar{T}=\bar{A}+\bar{E}$ and $(\bar{A_{t}})_{t}$ and
$(\bar{T_{t}})_{t}$ are respectively maximized at $t=z^{(2)}$ and
$t=\hat{z}^{(2)}$, we have on the event
$\Omega_{0}=\\{\hat{v}_{1}\in\Upsilon,\|\bar{E}\|_{\infty}\geqslant\lambda_{1}\\}$
that

$\displaystyle\bar{A}_{z^{(2)}}-\bar{A}_{\hat{z}^{(2)}}$
$\displaystyle=(\bar{A}_{z^{(2)}}-\bar{T}_{\hat{z}^{(2)}})+(\bar{T}_{z^{(2)}}-\bar{T}_{\hat{z}^{(2)}})+(\bar{T}_{\hat{z}^{(2)}}-\bar{A}_{\hat{z}^{(2)}})$
$\displaystyle\leqslant|\bar{A}_{z^{(2)}}-\bar{A}_{\hat{z}^{(2)}}|+|\bar{T}_{\hat{z}^{(2)}}-\bar{A}_{\hat{z}^{(2)}}|\leqslant
2\lambda_{1}.$

Hence, by Wang and Samworth (2018, Lemma 7 in the online supplement), on the
event $\Omega_{0}$, we have that

$\frac{|\hat{z}^{(2)}-z^{(2)}|}{n_{1}\tau}\leqslant\frac{3\sqrt{6}\lambda_{1}}{\bar{\theta}(n_{1}\tau)^{1/2}}\leqslant\frac{6\sqrt{2}\sigma(1+\sqrt{4\log
n_{1}})}{\vartheta\sqrt{n\tau}}.$

Now we define the event

$\Omega_{1}=\biggl{\\{}\biggl{|}\sum_{r=1}^{s}\bar{W_{t}}-\sum_{r=1}^{t}\bar{W}\biggr{|}\leqslant\lambda_{1}\sqrt{|s-t|},\quad
for\text{ }all\text{ }0\leqslant t\leqslant
n_{1},s\in\\{0,z^{(2)},n_{1}\\}\biggr{\\}}.$

By Wang and Samworth (2018, Lemma 5), we have that

$\mathbb{P}(\Omega_{1}^{\mathrm{c}})\leqslant 4e^{-\lambda_{1}^{2}/4}\\{2\log
n_{1}+\log z^{(2)}+\log(n_{1}-z^{(2)}\\}\leqslant 16\log
ne^{-\lambda_{1}^{2}/4}\leqslant\frac{16\log n}{n}.$ (16)

Following the proof of Theorem 1 of Wang and Samworth (2018), we have on
$\Omega_{0}\cap\Omega_{1}$ that

$1\leqslant\frac{6\sqrt{3}\lambda_{1}}{\bar{\theta}|\hat{z}^{(2)}-z^{(2)}|^{1/2}}+\frac{12\sqrt{6}\lambda_{1}}{\bar{\theta}(n_{1}\tau)^{1/2}}\\\
\leqslant\frac{12\sqrt{2}\sigma(1+\sqrt{4\log
n})}{\vartheta\sqrt{|z-\hat{z}}|}+\frac{48\sigma(1+\sqrt{4\log
n})}{\vartheta\sqrt{n\tau}}$

From (3) for $C\geqslant 96$, we have on $\Omega_{0}\cap\Omega_{1}$ that

$|\hat{z}-z|\leqslant C^{\prime}\vartheta^{-2}\sigma^{2}(1+\sqrt{4\log n}).$

Finally, by (14), (15) and (16) we have that

$\mathbb{P}(\Omega_{0}\cap\Omega_{1})\geqslant
1-\frac{1}{n_{1}G}-\frac{1}{n_{1}}-\frac{16\log n}{n}\geqslant 1-\frac{20\log
n}{n},$

as desired. ∎

##  7 Ancillary results

We collect in this section all ancillary propositions and lemmas used in the
paper. For all results in this section, we assume that we are given a grouping
$(\mathcal{J}_{g})_{g\in[G]}$ of $[p]$ and the associated group norm
$\|\cdot\|_{\mathrm{grp}}$. It is useful to define the following counterpart
to the group norm. For any $R\in\mathbb{R}^{p\times n}$ and a grouping
$(\mathcal{J}_{g})_{g\in[G]}$ of $[p]$, we define

$\left\lVert
R\right\rVert_{\mathrm{grp*}}=p_{g}^{-1/2}\max_{g\in[G]}\max_{t\in[n]}\|R_{\mathcal{J}_{g},t}\|_{2}.$
(17)

######  Lemma 4.

The norm $\|\cdot\|_{\mathrm{grp*}}$ is a dual to $\|\cdot\|_{\mathrm{grp}}$
with respect to the inner product $\langle\cdot,\cdot\rangle$ on
$\mathbb{R}^{p\times n}$.

###### Proof.

To prove the lemma, it suffices to show that $\left\lVert
M\right\rVert_{\mathrm{grp}}=\sup\limits_{\left\lVert
R\right\rVert_{\mathrm{grp*}}\leqslant 1}\langle R,M\rangle$ for all
$M\in\mathbb{R}^{p\times(n-1)}$. First, for any
$M\in\mathbb{R}^{p\times(n-1)}$, let $M_{\mathcal{J}_{g},t}$ be the $t$th
column of $M_{\mathcal{J}_{g}}$. Define $\tilde{R}=\tilde{R}(M)$ such that

$\tilde{R}_{\mathcal{J}_{g},t}=\frac{p_{g}^{1/2}M_{\mathcal{J}_{g},t}}{\max\bigl{\\{}\|M_{\mathcal{J}_{g},t}\|_{2},1\bigr{\\}}}.$

Then,
$\|\tilde{R}\|_{\mathrm{grp*}}\leqslant\max_{g\in[G]}\max_{t\in[n-1]}p_{g}^{-1/2}p_{g}^{1/2}\frac{\|M_{\mathcal{J}_{g},t}\|_{2}}{\|M_{\mathcal{J}_{g},t}\|_{2}}=1$.
Hence,

$\displaystyle\sup_{\left\lVert R\right\rVert_{\mathrm{grp*}}\leqslant
1}\langle R,M\rangle\geqslant\langle\tilde{R},M\rangle$
$\displaystyle=\sum_{g=1}^{G}\sum_{t=1}^{n-1}p_{g}^{1/2}\frac{\langle
M_{\mathcal{J}_{g},t},M_{\mathcal{J}_{g},t}\rangle}{\|M_{\mathcal{J}_{g},t}\|_{2}}$
$\displaystyle=\sum_{g=1}^{G}\sum_{t=1}^{n-1}p_{g}^{1/2}\|M_{\mathcal{J}_{g},t}\|_{2}=\|M\|_{\mathrm{grp}}.$

On the other hand, for any $R$ such that $\|R\|_{\mathrm{grp*}}\leqslant 1$,
we have $\|R_{\mathcal{J}_{g},t}\|_{2}\leqslant p_{g}^{1/2}$ for all $g$ and
$t$. Consequently, by the Cauchy–Schwarz inequality,

$\displaystyle\langle R,M\rangle$
$\displaystyle=\sum_{g\in[G]}\sum_{t\in[n-1]}\langle
R_{\mathcal{J}_{g},t},M_{\mathcal{J}_{g},t}\rangle\leqslant\sum_{g\in[G]}\sum_{t\in[n-1]}\|R_{\mathcal{J}_{g},t}\|_{2}\|M_{\mathcal{J}_{g},t}\|_{2}$
$\displaystyle\leqslant\sum_{g\in[G]}\sum_{t\in[n-1]}p_{g}^{1/2}\|M_{\mathcal{J}_{g},t}\|_{2}=\|M\|_{\mathrm{grp}},$

thus establishing the result. ∎

######  Proposition 5.

Let $\mathcal{S}=\\{M\in\mathbb{R}^{p\times(n-1)}:\left\lVert
M\right\rVert_{\mathrm{F}}\leqslant 1\\}$. For
$T\in\mathbb{R}^{p\times(n-1)}$, $\lambda>0$, we have

$\operatorname*{arg\,max}_{M\in\mathcal{S}}\Big{\\{}\langle
T,M\rangle-\lambda\|M\|_{\mathrm{{grp}}}\Big{\\}}=\frac{T-R^{*}}{\left\lVert
T-R^{*}\right\rVert_{\mathrm{F}}},$

where $R^{*}$ satisfies
$R_{\mathcal{J}_{g},t}^{*}=T_{\mathcal{J}_{g},t}\min\big{\\{}\frac{\lambda
p_{g}^{1/2}}{\|T_{\mathcal{J}_{g},t}\|_{\mathrm{F}}},1\big{\\}}$.

###### Proof.

Define functions
$h:\mathbb{R}^{p\times(n-1)}\times\mathbb{R}^{p\times(n-1)}\to\mathbb{R}$ and
$f,g:\mathbb{R}^{p\times(n-1)}\to\mathbb{R}$ such that for
$M,R\in\mathbb{R}^{p\times(n-1)}$, $h(M,R)=\langle T-\lambda R,M\rangle$ and
$f(M)=\inf_{\|R\|_{\mathrm{grp*}}\leqslant 1}h(M,R)$ and
$g(R)=\sup_{M\in\mathcal{S}}h(M,R)$. By (17) and Lemma 4, we have that

$\langle T,M\rangle-\lambda\|M\|_{\mathrm{grp}}=\langle
T,M\rangle-\lambda\sup_{\|R\|_{\mathrm{grp*}}\leqslant 1}\langle
R,M\rangle=\inf\limits_{\left\lVert R\right\rVert_{\mathrm{grp*}}\leqslant
1}\langle T-\lambda R,M\rangle=f(M).$

By the minimax equality theorem (Fan, 1953, Theorem 1), we obtain that

$\sup_{M\in\mathcal{S}}f(M)=\sup_{M\in\mathcal{S}}\inf_{\|R\|_{\mathrm{grp*}}\leqslant
1}h(M,R)=\inf_{\|R\|_{\mathrm{grp*}}\leqslant
1}\sup_{M\in\mathcal{S}}h(M,R)=\inf_{\|R\|_{\mathrm{grp*}}\leqslant 1}g(R).$

Observe that $g(R)=\|T-\lambda R\|_{\mathrm{F}}$. To find the optimiser
$R^{*}\in\operatorname*{arg\,min}_{\|R\|_{\mathrm{grp*}}\leqslant
1}\|T-\lambda R\|_{\mathrm{F}}$, we consider the $G$ groups individually. For
each group $g$, and in the $t$th column, if
$\|T_{\mathcal{J}_{g},t}\|_{2}\leqslant\lambda p_{g}^{1/2}$, then
$R^{*}_{\mathcal{J}_{g},t}=T_{\mathcal{J}_{g},t}/\lambda$; and if
$\|T_{\mathcal{J}_{g},t}\|_{2}>\lambda p_{g}^{1/2}$, then
$R^{*}_{\mathcal{J}_{g},t}=p_{g}^{1/2}T_{\mathcal{J}_{g},t}/\|T_{\mathcal{J}_{g},t}\|_{2}$.
Since the minimizer of $g(R)$ is unique, we have that

$\operatorname*{arg\,max}_{M\in\mathcal{S}}f(M)=\operatorname*{arg\,max}_{M\in\mathcal{S}}h(M,R^{*})=\frac{T-\lambda
R^{*}}{\|T-\lambda R^{*}\|_{\mathrm{F}}},$

as desired. ∎

######  Lemma 6.

For any $A,B\in\mathbb{R}^{p\times n}$, we have $\langle
A,B\rangle\leqslant\|A\|_{\mathrm{grp}}\|B\|_{\mathrm{grp*}}$.

###### Proof.

By Cauchy–Schwarz inequality, we have that

$\displaystyle\langle A,B\rangle$ $\displaystyle=\sum_{g,t}\langle
A_{\mathcal{J}_{g},t},B_{\mathcal{J}_{g},t}\rangle\leqslant\sum_{g\in[G],t\in[n]}\|A_{\mathcal{J}_{g},t}\|_{\mathrm{F}}\|B_{\mathcal{J}_{g},t}\|_{\mathrm{F}}$
$\displaystyle\leqslant\biggl{(}\sum_{g\in[G],t\in[n]}p_{g}^{1/2}\|A_{\mathcal{J}_{g},t}\|_{\mathrm{F}}\biggr{)}\biggl{(}\max_{g\in[G],t\in[n]}p_{g}^{-1/2}\|B_{\mathcal{J}_{g},t}\|_{\mathrm{F}}\biggr{)}=\|A\|_{\mathrm{grp}}\|B\|_{\mathrm{grp*}}.$

as desired. ∎

######  Lemma 7.

Let $p_{g}=|\mathcal{J}_{g}|$ and suppose further that there exists a
universal constant $C_{1}>0$, such that
$\max_{j\in[p]}|\\{g:j\in\mathcal{J}_{g}\\}|\leqslant C_{1}$. Then, for any
$M\in\mathbb{R}^{p\times n}$, we have
$\|M\|_{\mathrm{grp}}\leqslant(C_{1}n\sum_{g}p_{g})^{1/2}\|M\|_{\mathrm{F}}$.

###### Proof.

Define $m$ with
$m_{\mathcal{J}_{g},t}=\|M_{\mathcal{J}_{g},t}\|_{\mathrm{F}}$. Then by
applying the Cauchy–Schwarz inequality twice, we have

$\displaystyle\|M\|_{\mathrm{grp}}$
$\displaystyle=\sum_{g\in[G]}p_{g}^{1/2}\sum_{t=1}^{n}\|M_{\mathcal{J}_{g},t}\|_{2}\leqslant\sum_{g\in[G]}(np_{g})^{1/2}\|M_{\mathcal{J}_{g}}\|_{\mathrm{F}}$
$\displaystyle\leqslant\sqrt{n}\biggl{(}\sum_{g\in[G]}p_{g}\biggr{)}^{1/2}\biggl{(}\sum_{g\in[G]}\|M_{\mathcal{J}_{g}}\|_{\mathrm{F}}^{2}\biggr{)}^{1/2}\leqslant\biggl{(}C_{1}n\sum_{g\in[G]}p_{g}\biggr{)}^{1/2}\|M\|_{\mathrm{F}},$

as desired. ∎

######  Proposition 8.

Let $p_{g}=|\mathcal{J}_{g}|$ and suppose further that there exists a
universal constant $C_{1}>0$, such that
$\max_{j\in[p]}|\\{g:j\in\mathcal{J}_{g}\\}|\leqslant C_{1}$. Let $A$ be a
rank one matrix with $A=\delta vu^{\top}$ for $\delta>0$,
$\|v\|_{2}=\|u\|_{2}=1$ and $\sum_{g:v_{\mathcal{J}_{g}}\neq 0}p_{g}\leqslant
k$. Suppose $T\in\mathbb{R}^{p\times(n-1)}$ satisfies
$\|T-A\|_{\mathrm{grp*}}\leqslant\lambda$ for some $\lambda>0$, and let
$\mathcal{S}=\\{M\in\mathbb{R}^{p\times(n-1)}:\|M\|_{\mathrm{F}}\leqslant
1\\}$. Then, for any

$\hat{M}\in\operatorname*{arg\,max}_{M\in\mathcal{S}}\big{\\{}\langle
T,M\rangle-\lambda\|M\|_{\mathrm{grp}}\bigr{\\}},$

we have

$\|vu^{\top}-\hat{M}\|_{\mathrm{F}}\leqslant\frac{4\lambda(C_{1}nk)^{1/2}}{\delta},$

and

$\max\\{\sin\angle(v,\hat{v}),\sin\angle(u,\hat{u})\\}\leqslant\frac{8\lambda(C_{1}nk)^{1/2}}{\delta}.$

###### Proof.

Define $\mathcal{G}_{0}=\\{g:v_{\mathcal{J}_{g}}\neq 0\\}$. Since
$vu^{\top}\in\mathcal{S}$, from the basic inequality, we have

$\langle
T,vu^{\top}\rangle-\lambda\|vu^{\top}\|_{\mathrm{grp}}\leqslant\langle
T,\hat{M}\rangle-\lambda\|\hat{M}\|_{\mathrm{grp}}.$ (18)

When $\|A-T\|_{\mathrm{grp*}}\leqslant\lambda$, or equivalently,
$p_{g}^{-1/2}\|A_{\mathcal{J}_{g},t}-T_{\mathcal{J}_{g},t}\|_{2}\leqslant\lambda$
for all $g\in[G]$ and $t\in[n-1]$, we have by Wang and Samworth (2018, Lemma
2) and (18) that

$\displaystyle\|vu^{\top}-\hat{M}\|_{\mathrm{F}}^{2}$
$\displaystyle\leqslant\frac{2}{\delta}\langle
A,vu^{\top}-\hat{M}\rangle\leqslant\frac{2}{\delta}\bigl{(}\langle
T,vu^{\top}-\hat{M}\rangle+\langle A-T,vu^{\top}-\hat{M}\rangle\bigr{)}$
$\displaystyle\leqslant\frac{2\lambda}{\delta}\bigl{(}\|vu^{\top}\|_{\mathrm{grp}}-\|\hat{M}\|_{\mathrm{grp}}+\|vu^{\top}-\hat{M}\|_{\mathrm{grp}}\bigr{)}$
$\displaystyle=\frac{4\lambda}{\delta}\sum_{g\in\mathcal{G}_{0}}\sum_{t\in[n-1]}\|(vu^{\top}-\hat{M})_{\mathcal{J}_{g},t}\|_{2}\leqslant\frac{4\lambda(C_{1}nk)^{1/2}}{\delta}\|vu^{\top}-\hat{M}\|_{\mathrm{F}},$

where we used Lemma 6 in the penultimate inequality and Lemma 7 in the final
bound. This proves the first claim of the proposition, and the second claim
follows from the first by the same argument as used in Wang and Samworth
(2018, online supplement (18) and (19)). ∎

######  Proposition 9.

Let $W$ be an $p\times n$ random matrix with independent $N(0,\sigma^{2})$
entries and set $E=\mathcal{T}(W)$. Let $p_{g}=|\mathcal{J}_{g}|$ with
$p_{*}=\min_{g\in[G]}p_{g}$. For any $\delta\in(0,1)$ and
$\lambda=\sigma(1+\sqrt{2p_{*}^{-1}\log(1/\delta)})$, we have that

$\mathbb{P}(\|E\|_{\mathrm{grp*}}>\lambda)\leqslant(n-1)G\delta.$

###### Proof.

By the definition of the CUSUM transformation $\mathcal{T}$ in (5), we have
that $E_{\mathcal{J}_{g},t}\sim N(0,\sigma^{2}I_{p_{g}})$, and
$\|E_{\mathcal{J}_{g},t}/\sigma\|_{2}^{2}\sim\chi_{p_{g}}^{2}$ for every
$g\in[G]$ and $t\in[n-1]$. Consequently, by a union bound, we have

$\displaystyle\mathbb{P}(\|E\|_{\mathrm{grp*}}>\lambda)$
$\displaystyle\leqslant\sum_{g\in[G]}\sum_{t\in[n-1]}\mathbb{P}(\|E_{\mathcal{J}_{g},t}\|_{2}^{2}>p_{g}\lambda^{2})$
$\displaystyle=\sum_{g\in[G]}\sum_{t\in[n-1]}\mathbb{P}\biggl{\\{}\frac{\|E_{\mathcal{J}_{g},t}\|_{2}^{2}}{\sigma^{2}p_{g}}>\Bigl{(}1+\sqrt{2p_{*}^{-1}\log(1/\delta)}\Bigr{)}^{2}\biggr{\\}}$
$\displaystyle\leqslant\sum_{g\in[G]}\sum_{t\in[n-1]}\mathbb{P}\biggl{\\{}\frac{\|E_{\mathcal{J}_{g},t}\|_{2}^{2}}{\sigma^{2}p_{g}}>1+2\sqrt{\frac{\log(1/\delta)}{p_{g}}}+\frac{2(\log{1/\delta})}{p_{g}}\biggr{\\}}$
$\displaystyle\leqslant(n-1)G\delta,$

as desired, where we used Laurent and Massart (2000, Lemma 1) in the final
inequality. ∎

## References

  * Aston and Kirch (2012) Aston, J. A. D. and Kirch, C. (2012) Evaluating stationarity via change-point alternatives with applications to fMRI data. _Ann. Appl. Stat._ , 6, 1906–1948.
  * Cai et al. (2019) Cai, T. T., Zhang, A. and Zhou, Y. (2019) Sparse group lasso: Optimal sample complexity, convergence rate, and statistical inference. _arXiv preprint_ , arxiv:1909.09851.
  * Cho (2016) Cho, H. (2016) Change-point detection in panel data via double CUSUM statistic. _Electron. J. Stat._ , 10, 2000-2038.
  * Cho and Fryzlewicz (2015) Cho, H. and Fryzlewicz, P. (2015) Multiple-change-point detection for high dimensional time series via sparsified binary segmentation. _J. R. Stat. Soc. Ser. B_ , 77, 475–507.
  * Davis and Kahan (1970) Davis, C. and Kahan, W. M. (1970) The rotation of eigenvectors by a perturbation. III. _SIAM J. Numer. Anal._ , 7, 1–46.
  * Enikeeva and Harchaoui (2019) Enikeeva, F. and Harchaoui, Z. (2019) High-dimensional change-point detection under sparse alternatives. _Ann. Statist._ , 47, 2051–2079.
  * Fan (1953) Fan, K. (1953) Minimax theorems. _Proc. Natl. Acad. Sci. USA_ , 39, 42–47.
  * Frank and Wolfe (1956) Frank, M. and Wolfe, P. (1956) An algorithm for quadratic programming. _Naval Res. Logist._ , 3, 95–-110.
  * Frick, Munk and Sieling (2014) Frick, K., Munk, A. and Sieling, H. (2014) Multiscale change-point inference. _J. Roy. Statist. Soc., Ser. B_ , 76, 495–580.
  * Fryzlewicz (2014) Fryzlewicz, P. (2014) Wild binary segmentation for multiple change-point detection. _Ann. Statist._ , 42, 2243–2281.
  * Hanlon and Anderson (2009) Hanlon, M. and Anderson, R. (2009) Real-time gait event detection using wearable sensors. _Gait & Posture_, 30, 523–527.
  * Horváth and Hušková (2012) Horváth, L. and Hušková, M. (2012) Change-point detection in panel data. _J. Time Ser. Anal._ , 33, 631–648.
  * Hubert and Arabie (1985) Hubert, L. and Arabie, P. (1985) Comparing partitions. _J. Classification_ , 2, 193–-218.
  * Jirak (2015) Jirak, M. (2015) Uniform change-point tests in high dimension. _Ann. Statist._ , 43, 2451–2483.
  * Killick, Fearnhead and Eckley (2012) Killick, R., Fearnhead, P. and Eckley, I. A. (2012) Optimal detection of change-points with a linear computational cost. _J. Amer. Stat. Assoc._ , 107, 1590–1598.
  * Laurent and Massart (2000) Laurent, B. and Massart, P. (2000) Adaptive estimation of a quadratic functional by model selection. _Ann. Statist._ , 28, 1302–1338.
  * Massart (2007) Massart, P. (2007) _Concentration Inequalities and Model Selection_ , Springer, Berlin.
  * Peng, Leckie and Ramamohanarao (2004) Peng, T., Leckie, C. and Ramamohanarao, K. (2004) Proactively detecting distributed denial ofservice attacks using source IP address monitoring. In Mitrou, N., Kontovasilis, K., Rouskas, G. N., Iliadis, I. and Merakos, L. eds, _Networking 2004_ , pp. 771–782. Springer-Verlag, Berlin.
  * Rand (1971) Rand, W. M. (1971) Objective criteria for the evaluation of clustering methods. _J. Amer. Statist. Assoc._ , 66, 846–-850.
  * Simon et al. (2020) Simon, N, Friedman, J, Hastie, T and Tibshirani, R (2013) A sparse-group lasso. _J. Comput. Graph. Statist._ , 22, 231–245.
  * 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_. Cambridge University Press, Cambridge. 210–268.
  * Wang and Leng (2008) Wang, H and Leng, C (2008) A note on adaptive group lasso. _Comput. Statist. Data Anal._ 52(12), 5277–5286.
  * 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.
  * Yu (1997) Yu, B. (1997) Assouad, Fano and Le Cam. In Pollard, D., Torgersen, E. and Yang G. L. (Eds.) _Festschrift for Lucien Le Cam: Research Papers in Probability and Statistics_ , 423–435. Springer, New York.
  * Yuan and Lin (2006) Yuan, M. and Lin, Y. (2006) Model selection and estimation in regression with grouped variables. _J. Roy. Statist. Soc., Ser. B_ , 68, 49–67.


######  Theorem 9 (Variance of MOP-$\alpha$ Gradient Estimate).

When $\alpha\in(0,1)$, the variance of the gradient estimate from MOP-$\alpha$
is

$\mathrm{Var}\big{(}\nabla_{\theta}\hat{\ell}^{\alpha}(\theta)\big{)}\lesssim\min_{k\leq
N}\left(\frac{k^{2}G^{\prime}(\theta)^{2}Np}{(1-\alpha)^{2}J}+\frac{\alpha^{k}}{1-\alpha}Np\,G^{\prime}(\theta)^{2}\right).$
(172)

###### Proof.

Using the derivative of the logarithm and that
$w_{n,j}^{F,\theta,\alpha}=w_{n,j}^{F,\theta,\alpha,k}=1$ when $\theta=\phi$,

$\displaystyle\left\lVert\nabla_{\theta}\log\left(\sum_{j=1}^{J}w_{n,j}^{F,\theta,\alpha}\right)-\nabla_{\theta}\log\left(\sum_{j=1}^{J}w_{n,j}^{F,\theta,\alpha,k}\right)\right\rVert_{2}^{2}$
(173)
$\displaystyle=\left\lVert\frac{\nabla_{\theta}\sum_{j=1}^{J}w_{n,j}^{F,\theta,\alpha}}{{\sum_{j=1}^{J}w_{n,j}^{F,\theta,\alpha}}}-\frac{\nabla_{\theta}\sum_{j=1}^{J}w_{n,j}^{F,\theta,\alpha,k}}{{\sum_{j=1}^{J}w_{n,j}^{F,\theta,\alpha,k}}}\right\rVert_{2}^{2}$
(174)
$\displaystyle=\left\lVert\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}w_{n,j}^{F,\theta,\alpha}-\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}w_{n,j}^{F,\theta,\alpha,k}\right\rVert_{2}^{2}$
(175)
$\displaystyle=\left\lVert\frac{1}{J}\sum_{j=1}^{J}\frac{\nabla_{\theta}w_{n,j}^{F,\theta,\alpha}}{w_{n,j}^{F,\theta,\alpha}}-\frac{1}{J}\sum_{j=1}^{J}\frac{\nabla_{\theta}w_{n,j}^{F,\theta,\alpha,k}}{w_{n,j}^{F,\theta,\alpha,k}}\right\rVert_{2}^{2}$
(176)
$\displaystyle=\left\lVert\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}\log\left(w_{n,j}^{F,\theta,\alpha}\right)-\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}\log\left(w_{n,j}^{F,\theta,\alpha,k}\right)\right\rVert_{2}^{2}.$
(177)

This lets us bound the cumulative weight discrepancies by

$\displaystyle\left\lVert\nabla_{\theta}\log\left(\sum_{j=1}^{J}w_{n,j}^{F,\theta,\alpha}\right)-\nabla_{\theta}\log\left(\sum_{j=1}^{J}w_{n,j}^{F,\theta,\alpha,k}\right)\right\rVert_{2}^{2}$
(178)
$\displaystyle=\left\lVert\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}\log\left(w_{n,j}^{F,\theta,\alpha}\right)-\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}\log\left(w_{n,j}^{F,\theta,\alpha,k}\right)\right\rVert_{2}^{2}$
(179)
$\displaystyle=\left\lVert\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}\left(\log\left(w_{n,j}^{F,\theta,\alpha}\right)-\log\left(w_{n,j}^{F,\theta,\alpha,k}\right)\right)\right\rVert_{2}^{2}$
(180)
$\displaystyle=\left\lVert\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}\left(\sum_{i=1}^{n}\alpha^{(n-i)}\log\left(\frac{g_{i,j}^{A,F,\theta}}{g_{i,j}^{A,F,\phi}}\right)-\sum_{i=n-k}^{n}{\alpha^{(n-i)}}\log\left(\frac{g_{i,j}^{A,F,\theta}}{g_{i,j}^{A,F,\phi}}\right)\right)\right\rVert_{2}^{2}$
(181)
$\displaystyle=\left\lVert\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}\left(\sum_{i=1}^{n-k}\alpha^{(n-i)}\log\left(\frac{g_{i,j}^{A,F,\theta}}{g_{i,j}^{A,F,\phi}}\right)\right)\right\rVert_{2}^{2}$
(182)
$\displaystyle\leq\frac{1}{J}\sum_{j=1}^{J}\sum_{i=1}^{n-k}\alpha^{(n-i)}\left\lVert\nabla_{\theta}\log\left(g_{i,j}^{A,F,\theta}\right)\right\rVert_{2}^{2}$
(183)
$\displaystyle\leq\frac{1}{J}\sum_{j=1}^{J}\sum_{i=1}^{n-k}\alpha^{(n-i)}p\,G^{\prime}(\theta)^{2}$
(184) $\displaystyle\leq
p\,G^{\prime}(\theta)^{2}\frac{\alpha^{k}-\alpha^{n}}{1-\alpha},$ (185)

where the second-last line follows from Assumption (A3). We can then bound
$\|\nabla_{\theta}\hat{\ell}^{\alpha}(\theta)-\nabla_{\theta}\hat{s}_{k}^{\alpha}(\theta)\|$
as follows:

$\displaystyle\left\lVert\nabla_{\theta}\hat{\ell}^{\alpha}(\theta)-\nabla_{\theta}\hat{s}_{k}^{\alpha}(\theta)\right\rVert_{2}^{2}$
(186)
$\displaystyle=\left\lVert\sum_{n=1}^{N}\nabla_{\theta}\log\left(\frac{\sum_{j=1}^{J}\prod_{i=1}^{n}\left(\frac{g_{i,j}^{A,F,\theta}}{g_{i,j}^{A,F,\phi}}\right)^{\alpha^{(n-i)}}}{\sum_{j=1}^{J}\prod_{i=1}^{n-1}\left(\frac{g_{i,j}^{A,F,\theta}}{g_{i,j}^{A,F,\phi}}\right)^{\alpha^{(n-i)}}}\right)-\sum_{n=1}^{N}\nabla_{\theta}\log\left(\frac{\sum_{j=1}^{J}\prod_{i=n-k}^{n}\left(\frac{g_{i,j}^{A,F,\theta}}{g_{i,j}^{A,F,\phi}}\right)^{\alpha^{(n-i)}}}{\sum_{j=1}^{J}\prod_{i=n-k}^{n-1}\left(\frac{g_{i,j}^{A,F,\theta}}{g_{i,j}^{A,F,\phi}}\right)^{\alpha^{(n-i)}}}\right)\right\rVert_{2}^{2}$
$\displaystyle=\Bigg{\lVert}\sum_{n=1}^{N}\nabla_{\theta}\Bigg{(}\log\left(\sum_{j=1}^{J}\prod_{i=1}^{n}\left(\frac{g_{i,j}^{A,F,\theta}}{g_{i,j}^{A,F,\phi}}\right)^{\alpha^{(n-i)}}\right)-\log\left(\sum_{j=1}^{J}\prod_{i=1}^{n-1}\left(\frac{g_{i,j}^{A,F,\theta}}{g_{i,j}^{A,F,\phi}}\right)^{\alpha^{(n-i)}}\right)$
$\displaystyle\qquad\qquad-\log\left(\sum_{j=1}^{J}\prod_{i=n-k}^{n}\left(\frac{g_{i,j}^{A,F,\theta}}{g_{i,j}^{A,F,\phi}}\right)^{\alpha^{(n-i)}}\right)+\log\left(\sum_{j=1}^{J}\prod_{i=n-k}^{n-1}\left(\frac{g_{i,j}^{A,F,\theta}}{g_{i,j}^{A,F,\phi}}\right)^{\alpha^{(n-i)}}\right)\Bigg{)}\Bigg{\rVert}_{2}^{2}$
(187)
$\displaystyle=\left\lVert\sum_{n=1}^{N}\nabla_{\theta}\Bigg{(}\\!\log\\!\left(\sum_{j=1}^{J}w_{n,j}^{F,\theta,\alpha}\right)\\!-\\!\log\\!\left(\sum_{j=1}^{J}w_{n,j}^{F,\theta,\alpha,k}\right)\\!-\\!\log\\!\left(\sum_{j=1}^{J}w_{n-1,j}^{A,F,\theta,\alpha}\right)\\!+\\!\log\\!\left(\sum_{j=1}^{J}w_{n-1,j}^{A,F,\theta,\alpha,k}\right)\\!\Bigg{)}\right\rVert_{2}^{2}$
(188)
$\displaystyle\leq\sum_{n=1}^{N}\left\lVert\nabla_{\theta}\left(\log\left(\sum_{j=1}^{J}w_{n,j}^{F,\theta,\alpha}\right)-\log\left(\sum_{j=1}^{J}w_{n,j}^{F,\theta,\alpha,k}\right)\right)\right\lVert_{\infty}$
$\displaystyle\qquad\qquad+\sum_{n=1}^{N}\left\lVert\nabla_{\theta}\left(\log\left(\sum_{j=1}^{J}w_{n-1,j}^{A,F,\theta,\alpha}\right)+\log\left(\sum_{j=1}^{J}w_{n-1,j}^{A,F,\theta,\alpha,k}\right)\right)\right\rVert_{2}^{2},$
(189)

and each of these two terms is bounded by
$\frac{\alpha^{k}}{1-\alpha}Np\,G^{\prime}(\theta)^{2}$. So we now know that

$\left\lVert\nabla_{\theta}\hat{\ell}^{\alpha}(\theta)-\nabla_{\theta}\hat{s}_{k}^{\alpha}(\theta)\right\rVert_{2}^{2}\leq\frac{2\alpha^{k}}{1-\alpha}Np\,G^{\prime}(\theta)^{2},$
(190)

which is our desired error bound. Now, we bound the variance of
$\nabla_{\theta}\hat{s}_{k}^{\alpha}(\theta)$. Recall that we defined

$\displaystyle\nabla_{\theta}\hat{s}_{k}^{\alpha}(\theta)$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{N}\nabla_{\theta}\log\left(\frac{\sum_{j=1}^{J}\prod_{i=n-k}^{n}\left(\frac{g_{i,j}^{A,F,\theta}}{g_{i,j}^{A,F,\phi}}\right)^{\alpha^{(n-i)}}}{\sum_{j=1}^{J}\prod_{i=n-k}^{n-1}\left(\frac{g_{i,j}^{A,F,\theta}}{g_{i,j}^{A,F,\phi}}\right)^{\alpha^{(n-i)}}}\right)$
(191) $\displaystyle=$
$\displaystyle\sum_{n=1}^{N}\nabla_{\theta}\left(\log\left(\sum_{j=1}^{J}w_{n,j}^{F,\theta,\alpha,k}\right)-\log\left(\sum_{j=1}^{J}w_{n-1,j}^{A,F,\theta,\alpha,k}\right)\right).$
(192)

We can then decompose this expression with the derivative of the logarithm
while noting that $w_{n,j}^{F,\theta,\alpha,k}=1$ whenever $\theta=\phi$, to
see that

$\displaystyle\nabla_{\theta}\hat{s}_{k}^{\alpha}(\theta)$
$\displaystyle=\sum_{n=1}^{N}\nabla_{\theta}\left(\log\left(\sum_{j=1}^{J}w_{n,j}^{F,\theta,\alpha,k}\right)-\log\left(\sum_{j=1}^{J}w_{n-1,j}^{A,F,\theta,\alpha,k}\right)\right)$
(193)
$\displaystyle=\sum_{n=1}^{N}\left(\frac{\sum_{j=1}^{J}\nabla_{\theta}w_{n,j}^{F,\theta,\alpha,k}}{\sum_{j=1}^{J}w_{n,j}^{F,\theta,\alpha,k}}-\frac{\sum_{j=1}^{J}\nabla_{\theta}w_{n-1,j}^{A,F,\theta,\alpha,k}}{\sum_{j=1}^{J}w_{n-1,j}^{A,F,\theta,\alpha,k}}\right)$
(194)
$\displaystyle=\sum_{n=1}^{N}\left(\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}w_{n,j}^{F,\theta,\alpha,k}-\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}w_{n-1,j}^{A,F,\theta,\alpha,k}\right)$
(195)
$\displaystyle=\frac{1}{J}\sum_{n=1}^{N}\sum_{j=1}^{J}\nabla_{\theta}\left(w_{n,j}^{F,\theta,\alpha,k}-w_{n-1,j}^{A,F,\theta,\alpha,k}\right).$
(196)

It now follows that we need to bound the variance at a single timestep, namely

$\mathrm{Var}\left(\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}\left(w_{n,j}^{F,\theta,\alpha,k}-w_{n-1,j}^{A,F,\theta,\alpha,k}\right)\right).$
(197)

We use the derivative of the logarithm yet again to find, noting that
$w_{n,j}^{F,\theta,\alpha,k}=1$ when $\theta=\phi$, that

$\displaystyle\nabla_{\theta}w_{n,j}^{F,\theta,\alpha,k}=\frac{\nabla_{\theta}w_{n,j}^{F,\theta,\alpha,k}}{w_{n,j}^{F,\theta,\alpha,k}}=\nabla_{\theta}\log(w_{n,j}^{F,\theta,\alpha,k})$
$\displaystyle=\sum_{i=n-k}^{n}\alpha^{n-i}\nabla_{\theta}\log\left(\frac{g_{i,j}^{A,F,\theta}}{g_{i,j}^{A,F,\phi}}\right)$
(198)
$\displaystyle=\sum_{i=n-k}^{n}\alpha^{n-i}\nabla_{\theta}\log\left(g_{i,j}^{A,F,\theta}\right).$
(199)

Then,

$\displaystyle\mathrm{Var}\left(\nabla_{\theta}\big{(}w_{n,j}^{F,\theta,\alpha,k}-w_{n-1,j}^{A,F,\theta,\alpha,k}\big{)}\right)$
(200)
$\displaystyle\qquad=\mathrm{Var}\left(\nabla_{\theta}\left(\sum_{i=n-k}^{n}\alpha^{n-i}\nabla_{\theta}\log\left(g_{i,j}^{A,F,\theta}\right)-\sum_{i=n-k-1}^{n-1}\alpha^{n-i}\nabla_{\theta}\log\left(g_{i,j}^{A,F,\theta}\right)\right)\right)$
(201)
$\displaystyle\qquad=\mathrm{Var}\left(\sum_{i=n-k}^{n}\nabla_{\theta}\left(\alpha^{n-i}\log\left(g_{i,j}^{A,F,\theta}\right)-\alpha^{n-i+1}\log\left(g_{i-1,j}^{A,F,\theta}\right)\right)\right).$
(202)

Note that

$\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta):=\frac{1}{J}\sum_{j=1}^{J}\sum_{i=n-k}^{n}\nabla_{\theta}\left(\alpha^{n-i}\log\left(g_{i,j}^{A,F,\theta}\right)-\alpha^{n-i+1}\nabla_{\theta}\log\left(g_{i-1,j}^{A,F,\theta}\right)\right)$
(203)

is a function bounded by
$C\frac{G^{\prime}(\theta)(1-\alpha^{k})}{1-\alpha}=:CG^{\prime}(\theta)r$ in
each coordinate for some constant $C$ (by Assumption (A3)) that depends only
on $k$ timesteps, and not $n$. Subsequently, we suppress vector and matrix
notation for the coordinates of $\theta$ and its derivatives, and their
variances and covariances, with the variance of a vector interpreted as the
sum of the variance of its components. We invoke Lemma 2 of [42] to bound the
$L_{2}$ error of each of the $k$ (from time $n-k$ to time $n$) functionals
from its expectation under the posterior by
$Cr\sqrt{p}G^{\prime}(\theta)J^{-1/2}k$. That is, we have that

$\mathbb{E}\left[\|\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta)-\mathbb{E}_{\pi}\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta)\|_{2}^{2}\right]^{1/2}\leq\frac{Cr\sqrt{p}G^{\prime}(\theta)k}{\sqrt{J}}.$
(204)

By the bias-variance decomposition, this in turn bounds the variance at a
single timestep by $C^{2}r^{2}p\,G^{\prime}(\theta)^{2}J^{-1}k^{2}$.
Concretely, we have that

$\mathrm{Var}\left(\frac{1}{J}\sum_{j=1}^{J}\sum_{i=n-k}^{n}\nabla_{\theta}\left(\alpha^{n-i}\log\left(g_{i,j}^{A,F,\theta}\right)-\alpha^{n-i+1}\log\left(g_{i-1,j}^{A,F,\theta}\right)\right)\right)\lesssim\frac{r^{2}k^{2}p2G^{\prime}(\theta)^{2}}{J}.$
(205)

It now remains to bound the variance of
$\nabla_{\theta}\hat{s}_{k}^{\alpha}(\theta)$ by considering the covariance of
each of the $N$ terms that comprise it. We decompose

$\displaystyle\mathrm{Var}\big{(}\nabla_{\theta}\hat{s}_{k}^{\alpha}(\theta)\big{)}$
$\displaystyle=\mathrm{Var}\left(\sum_{n=1}^{N}\frac{1}{J}\sum_{j=1}^{J}\sum_{i=n-k}^{n}\nabla_{\theta}\left(\alpha^{n-i}\log\left(g_{i,j}^{A,F,\theta}\right)-\alpha^{n-i+1}\log\left(g_{i-1,j}^{A,F,\theta}\right)\right)\right)$
(206)
$\displaystyle=\sum_{n=1}^{N}\mathrm{Var}\big{(}\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta)\big{)}+2\sum_{m<n}\mathrm{Cov}\big{(}\nabla_{\theta}\hat{s}_{m,k}^{\alpha}(\theta),\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta)\big{)}$
(207)
$\displaystyle\lesssim\frac{r^{2}k^{2}p\,G^{\prime}(\theta)^{2}N}{J}+2\sum_{m<n}\mathrm{Cov}\big{(}\nabla_{\theta}\hat{s}_{m,k}^{\alpha}(\theta),\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta)\big{)}.$
(208)

Similarly to the proof of the MOP-$0$ case, we use Assumptions (A2) and (A3)
that ensure strong mixing. We know from Theorem 3 of [42] that when
$\textbf{M}_{n,n+k}$ is the $k$-step Markov operator from timestep $n$ and
$\beta_{\text{TV}}(M)=\sup_{x,y\in
E}\|M(x,\cdot)-M(y,\cdot)\|_{\mathrm{TV}}=\sup_{\mu,\nu\in\mathcal{P},\mu\neq\nu}\frac{\|\mu
M-\nu M\|_{\mathrm{TV}}}{\|\mu-\nu\|_{\mathrm{TV}}}$ is the Dobrushin
contraction coefficient of a Markov operator,

$\beta_{\text{TV}}(\textbf{M}_{n,n+k})\leq(1-\epsilon)^{\lfloor
k/(c\log(J))\rfloor},$ (209)

i.e. the mixing time of the particle filter is $O(\log(J))$, where $\epsilon$
and $c$ depend on $\bar{M},\underline{M},\bar{G},\underline{G}$ in (A2) and
(A3). By Lemma 9, the particle filter itself is strongly mixing, with
$\alpha$-mixing coefficients $a(l)\leq(1-\epsilon)^{\lfloor
l/(c\log(J))\rfloor}$. Therefore, functions of particles are strongly mixing
as well, with $\alpha$-mixing coefficients bounded by the original (to see
this, observe that the $\sigma$-algebra of the functionals is contained within
the original $\sigma$-algebra).

We now derive the $\alpha$-mixing coefficients of
$(\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta))_{n=1}^{N}$. Observe that
$\nabla_{\theta}\hat{s}_{n,k}^{\alpha}$ is strong mixing at lag $k+1$, as all
weights beyond $k$ timesteps are truncated. We therefore have that the mixing
time of $\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta)$ is $O(1+k+\log(J))$,
and that the $\alpha$-mixing coefficients for
$(\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta))_{n=1}^{N}$ are given by
$a(l)\leq(1-\epsilon)^{\lfloor l/(c(1+k+\log(J)))\rfloor}$. Therefore, by
Davydov’s inequality and Lemma 2 of [42], and noting that
$\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta)\leq CrG^{\prime}(\theta)$ by
Assumption (A3), by a similar argument to the MOP-0 case, we find that

$\mathrm{Cov}\big{(}\nabla_{\theta}\hat{s}_{m,k}^{\alpha}(\theta),\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta)\big{)}\lesssim(1-\epsilon)^{\frac{1}{2}\lfloor\frac{|n-m|}{c(1+k+\log(J))}\rfloor}\frac{r^{2}p\,G^{\prime}(\theta)^{2}}{J}.$
(210)

Concretely, noting that $\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta)\leq
CrG^{\prime}(\theta)$ by Assumption (A3), and, without loss of generality,
assuming
$\mathbb{E}[(\nabla_{\theta}\hat{s}_{m,k}^{\alpha}(\theta))^{4}]^{1/4}\leq\mathbb{E}[(\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta))^{4}]^{1/4}$,
we apply Davydov’s inequality to see that

$\displaystyle\mathrm{Cov}\big{(}\nabla_{\theta}\hat{s}_{m,k}^{\alpha}(\theta),\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta)\big{)}$
$\displaystyle\leq
a(n-m)^{1/2}\,\mathbb{E}\big{[}(\nabla_{\theta}\hat{s}_{m,k}^{\alpha}(\theta))^{4}\big{]}^{1/4}\,\mathbb{E}\big{[}(\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta))^{4}\big{]}^{1/4}$
(211) $\displaystyle\leq
a(n-m)^{1/2}\,\mathbb{E}\big{[}(\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta))^{4}\big{]}^{1/2}.$
(212)

To bound this, we use the fact that

$\mathbb{E}\big{[}(\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta))^{4}\big{]}=\mathbb{E}\Big{[}\big{(}\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta)-\mathbb{E}\big{[}\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta)]\big{)}^{4}\Big{]}+\mathbb{E}\big{[}(\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta))^{2}\big{]}^{2}$
(213)

alongside Lemma 2 of [42], which alongside the fact that

$\displaystyle\mathbb{E}\Big{[}\big{(}\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta)-\mathbb{E}[\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta)]\big{)}^{4}\Big{]}$
$\displaystyle\lesssim$
$\displaystyle\frac{r^{4}p^{2}G^{\prime}(\theta)^{4}}{J^{2}},$ (214)
$\displaystyle\mathbb{E}\Big{[}\big{(}\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta)-\mathbb{E}[\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta)]\big{)}^{2}\Big{]}$
$\displaystyle\lesssim$
$\displaystyle\frac{r^{2}p\,G^{\prime}(\theta)^{2}}{J},$ (215)

allows us to find that

$\mathbb{E}\Big{[}\big{(}\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta)\big{)}^{4}\Big{]}^{1/2}\lesssim\sqrt{\frac{r^{4}p^{2}G^{\prime}(\theta)^{4}}{J^{2}}+\left(\frac{r^{2}p\,G^{\prime}(\theta)^{2}}{J}\right)^{2}}=\frac{r^{2}p\,G^{\prime}(\theta)^{2}}{J},$
(216)

and conclude that

$\mathrm{Cov}\big{(}\nabla_{\theta}\hat{s}_{m,k}^{\alpha}(\theta),\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta)\big{)}\leq(1-\epsilon)^{\frac{1}{2}\lfloor\frac{|n-m|}{c(1+k+\log(J))}\rfloor}\frac{r^{2}p\,G^{\prime}(\theta)^{2}}{J}.$
(217)

Putting it all together, we see that

$\displaystyle\mathrm{Var}\big{(}\nabla_{\theta}\hat{s}_{k}(\theta)\big{)}$
$\displaystyle=\sum_{n=1}^{N}\mathrm{Var}\big{(}\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta)\big{)}+2\sum_{m<n}\mathrm{Cov}\big{(}\nabla_{\theta}\hat{s}_{m,k}^{\alpha}(\theta),\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta)\big{)}$
(218)
$\displaystyle\leq\frac{Nr^{2}pk^{2}G^{\prime}(\theta)^{2}}{J}+2\sum_{n=1}^{N-1}(N-n)(1-\epsilon)^{\frac{1}{2}\lfloor\frac{n}{c(1+k+\log(J))}\rfloor}\frac{r^{2}p\,G^{\prime}(\theta)^{2}}{J}$
(219)
$\displaystyle\leq\frac{r^{2}pk^{2}G^{\prime}(\theta)^{2}}{J}\left(N+2\sum_{n=1}^{N-1}(N-n)(1-\epsilon)^{\frac{1}{2}\lfloor\frac{n}{c(1+k+\log(J))}\rfloor}\right)$
(220) $\displaystyle\lesssim\frac{k^{2}r^{2}p\,G^{\prime}(\theta)^{2}N}{J}.$
(221)

We will now use this result to bound
$\mathrm{Var}\big{(}\nabla_{\theta}\hat{s}_{k}^{\alpha}(\theta)\big{)}$. For
random variables $X,Y$ where $|X-Y|$ is bounded almost surely by some $M$, we
have that

$X-\mathbb{E}X-2M\leq Y-\mathbb{E}Y\leq X-\mathbb{E}X+2M$ (222)

so,

$(Y-\mathbb{E}Y)^{2}\leq(X-\mathbb{E}X-2M)^{2}+(X-\mathbb{E}X+2M)^{2}.$ (223)

Taking expectations, we get

$\mathrm{Var}(Y)\leq\mathrm{Var}(X)+8M^{2}.$ (224)

It follows from this result and the result we proved earlier, namely that

$\big{\|}\nabla_{\theta}\hat{\ell}^{\alpha}(\theta)-\nabla_{\theta}\hat{s}_{k}^{\alpha}(\theta)\big{\|}_{2}^{2}\leq\frac{2\alpha^{k}}{1-\alpha}Np\,G^{\prime}(\theta)^{2},$
(225)

that the variance of MOP-$\alpha$ proper is bounded by, for any $k\leq N$,
that the variance of MOP-$\alpha$ proper is bounded by

$\mathrm{Var}\big{(}\nabla_{\theta}\hat{\ell}^{\alpha}(\theta)\big{)}\lesssim\frac{k^{2}r^{2}G^{\prime}(\theta)^{2}Np}{J}+\frac{\alpha^{k}}{1-\alpha}Np\,G^{\prime}(\theta)^{2}.$
(226)

As the above holds for any $k\leq N$,

$\mathrm{Var}\big{(}\nabla_{\theta}\hat{\ell}^{\alpha}(\theta)\big{)}\lesssim\min_{k\leq
N}\left(\frac{k^{2}p\,G^{\prime}(\theta)^{2}Np}{(1-\alpha)^{2}J}+\frac{\alpha^{k}}{1-\alpha}Np\,G^{\prime}(\theta)^{2}\right).$
(227)

∎

###  F.3 MOP-$\alpha$ MSE Bound

######  Theorem 10 (MSE of MOP-$\alpha$).

When $\alpha\in(0,1)$, the MSE of MOP-$\alpha$ is given by

$\mathbb{E}\big{\|}\nabla_{\theta}\ell(\theta)-\nabla_{\theta}\hat{\ell}^{\alpha}(\theta)\big{\|}_{2}^{2}\;\lesssim\;\min_{k\leq
N}Np\,G^{\prime}(\theta)^{2}\left(\frac{k^{2}}{J}+(1-\epsilon)^{\lfloor
k/(c\log(J))\rfloor}+k+\frac{\alpha^{k}+\alpha^{k+1}-\alpha}{1-\alpha}\right).$
(228)

###### Proof.

The broad idea is to decompose the MSE into three terms as in Equation 230.
The first term can be controlled by a mixing argument, the second term garners
$O(Nk)$ error, and the third term is $O\left(\frac{N}{1-\alpha}\right)$. It is
unclear whether the bias can be reduced further – the gradient is a sum of $N$
terms, and is therefore $O(N)$ itself.

The first term is the problem term. We will first bound the ground truth of
conditional scores from particle approximation conditional on correct
filtering distribution at time $n-k$, and then bound that from the particle
approximation resulting from an arbitrary filtering distribution at time
$n-k$. That is, where $\hat{\pi}_{n}$ is the particle approximation at time
$n$ of the posterior $\pi_{n}$ and we write
$\nabla_{\theta}\hat{s}_{n,k|\hat{\pi}_{n-k}=\pi_{n-k}}^{1}(\theta)$ for the
truncated MOP-$1$ conditional score estimate given the correct filtering
distribution at time $n-k$, we can first observe that

$\big{\|}\nabla_{\theta}\ell_{\theta}-\nabla_{\theta}\hat{\ell}^{\alpha}(\theta)\big{\|}_{2}^{2}=\left\lVert\sum_{n=1}^{N}\nabla_{\theta}\hat{\ell}_{n}^{\alpha}(\theta)-\sum_{n=1}^{N}\nabla_{\theta}\ell_{n}(\theta)\right\rVert_{2}^{2}\leq\sum_{n=1}^{N}\left\lVert\nabla_{\theta}\hat{\ell}_{n}^{\alpha}(\theta)-\nabla_{\theta}\ell_{n}(\theta)\right\rVert_{2}^{2},$
(229)

and decompose

$\displaystyle\sum_{n=1}^{N}\left\lVert\nabla_{\theta}\hat{\ell}_{n}^{\alpha}(\theta)-\nabla_{\theta}\ell_{n}(\theta)\right\rVert_{2}^{2}\;\leq\;\sum_{n=1}^{N}\|\nabla_{\theta}\ell_{n}(\theta)-\nabla_{\theta}\hat{s}_{n,k}^{1}(\theta)\|_{2}^{2}\;+$
$\displaystyle\hskip
56.9055pt\sum_{n=1}^{N}\|\nabla_{\theta}\hat{s}_{n,k}^{1}(\theta)-\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta)\|_{2}^{2}+\sum_{n=1}^{N}\|\nabla_{\theta}\hat{s}_{n,k}^{\alpha}(\theta)-\nabla_{\theta}\hat{\ell}_{n}^{\alpha}(\theta)\|_{2}^{2}.$
(230)

We bound the first term of Equation 230 ,
$\sum_{n=1}^{N}\|\nabla_{\theta}\ell_{n}(\theta)-\nabla_{\theta}\hat{s}_{n,k}^{1}(\theta)\|_{2}^{2}$,
by decomposing it into two terms,

$\displaystyle\sum_{n=1}^{N}\|\nabla_{\theta}\ell_{n}(\theta)-\nabla_{\theta}\hat{s}_{n,k}^{1}(\theta)\|_{2}^{2}$
$\displaystyle\leq$ (231)
$\displaystyle\sum_{n=1}^{N}\|\nabla_{\theta}\ell_{n}(\theta)-\nabla_{\theta}\hat{s}_{n,k|\hat{\pi}_{n-k}=\pi_{n-k}}^{1}(\theta)\|_{2}^{2}+\sum_{n=1}^{N}\|\nabla_{\theta}\hat{s}_{n,k|\hat{\pi}_{n-k}=\pi_{n-k}}^{1}(\theta)-\nabla_{\theta}\hat{s}_{n,k}^{1}(\theta)\|_{2}^{2}.$

The first term of Equation 230 is a particle approximation dependent on $k$
timesteps, so by Lemma 2 of [42], this is bounded by

$\mathbb{E}\|\nabla_{\theta}\ell_{n}(\theta)-\nabla_{\theta}\hat{s}_{n,k|\hat{\pi}_{n-k}=\pi_{n-k}}^{1}(\theta)\|_{2}^{2}\leq\frac{Cp\,G^{\prime}(\theta)^{2}k^{2}}{J}.$
(232)

Bounding the second term of Equation 230 amounts to bounding the difference
between functionals of two different particle measures that mix under the same
Markov kernel. Here, we use Assumptions (A2) and (A3) that ensure strong
mixing. We know from Theorem 3 of [42] that when $\textbf{M}_{n,n+k}$ is the
$k$-step Markov operator from timestep $n$ and
$\beta_{\text{TV}}(M)=\sup_{x,y\in
E}\|M(x,\cdot)-M(y,\cdot)\|_{\mathrm{TV}}=\sup_{\mu,\nu\in\mathcal{P},\mu\neq\nu}\frac{\|\mu
M-\nu M\|_{\mathrm{TV}}}{\|\mu-\nu\|_{\mathrm{TV}}}$ is the Dobrushin
contraction coefficient of a Markov operator,

$\beta_{\text{TV}}(\textbf{M}_{n,n+k})\leq(1-\epsilon)^{\lfloor
k/(c\log(J))\rfloor},$ (233)

i.e. the mixing time of the particle filter is $O(\log(J))$, where $\epsilon$
and $c$ depend on $\bar{M},\underline{M},\bar{G},\underline{G}$ in (A2) and
(A3).

Then, we can bound
$\mathbb{E}\|\nabla_{\theta}\hat{s}_{n,k|\hat{\pi}_{n-k}=\pi_{n-k}}^{1}(\theta)-\nabla_{\theta}\hat{s}_{n,k}^{1}(\theta)\|_{2}^{2}$
by

$\sup_{\mu,\nu\in\mathcal{P},\mu\neq\nu}\frac{\|\mu\textbf{M}_{n,n+k}-\nu\textbf{M}_{n,n+k}\|_{\mathrm{TV}}}{\|\mu-\nu\|_{\mathrm{TV}}}=\beta_{TV}(\textbf{M}_{n,n+k})\leq(1-\epsilon)^{\lfloor
k/(c\log(J))\rfloor},$ (234)

implying that

$\displaystyle\mathbb{E}\big{\|}\nabla_{\theta}\hat{s}_{n,k|\hat{\pi}_{n-k}=\pi_{n-k}}^{1}(\theta)-\nabla_{\theta}\hat{s}_{n,k}^{1}(\theta)\big{\|}_{2}^{2}$
(235) $\displaystyle\lesssim$
$\displaystyle\sup_{\mu,\nu}\sup_{\|\psi\|_{\infty}\leq
1/2}p\,G^{\prime}(\theta)^{2}\big{|}(\mu\textbf{M}_{n,n+k})(\psi)-(\nu\textbf{M}_{n,n+k})(\psi)\big{|}$
(236) $\displaystyle\leq$
$\displaystyle\sup_{\mu,\nu}p\,G^{\prime}(\theta)^{2}\big{\|}\mu\textbf{M}_{n,n+k}-\nu\textbf{M}_{n,n+k}\big{\|}_{\mathrm{TV}}$
(237) $\displaystyle\leq$ $\displaystyle
p\,G^{\prime}(\theta)^{2}(1-\epsilon)^{\lfloor
k/(c\log(J))\rfloor}\|\hat{\pi}_{n-k}-\pi_{n-k}\|_{\text{TV}}$ (238)
$\displaystyle\leq$ $\displaystyle
p\,G^{\prime}(\theta)^{2}(1-\epsilon)^{\lfloor
k/(c\log(J))\rfloor}\sup_{\|\psi\|_{\infty}\leq
1/2}\big{|}\hat{\pi}_{n-k}(\psi)-\pi_{n-k}(\psi)\big{|}$ (239)
$\displaystyle\lesssim$ $\displaystyle
p\,G^{\prime}(\theta)^{2}(1-\epsilon)^{\lfloor k/(c\log(J))\rfloor}.$ (240)

Therefore, we have that

$\displaystyle\mathbb{E}\big{\|}\nabla_{\theta}\ell(\theta)-\nabla_{\theta}\hat{s}_{k}^{1}(\theta)\big{\|}_{2}^{2}$
(241)
$\displaystyle\leq\sum_{n=1}^{N}\mathbb{E}\big{\|}\nabla_{\theta}\ell_{n}(\theta)-\nabla_{\theta}\hat{s}_{n,k|\hat{\pi}_{n-k}=\pi_{n-k}}^{1}(\theta)\big{\|}_{2}^{2}+\sum_{n=1}^{N}\mathbb{E}\big{\|}\nabla_{\theta}\hat{s}_{n,k|\hat{\pi}_{n-k}=\pi_{n-k}}^{1}(\theta)-\nabla_{\theta}\hat{s}_{n,k}^{1}(\theta)\big{\|}_{2}^{2}$
(242) $\displaystyle\lesssim
N\frac{Cp\,G^{\prime}(\theta)^{2}k^{2}}{J}+Np\,G^{\prime}(\theta)^{2}(1-\epsilon)^{\lfloor
k/(c\log(J))\rfloor}.$ (243)

Now that the first term, $\mathbb{E}\|\nabla_{\theta}\ell(\theta)\|_{2}^{2}$,
is taken care of, it remains to bound
$\mathbb{E}\|\nabla_{\theta}\hat{s}_{k}^{1}(\theta)-\nabla_{\theta}\hat{s}_{k}^{\alpha}(\theta)\|_{2}^{2}$
and
$\mathbb{E}\|\nabla_{\theta}\hat{s}_{k}^{\alpha}(\theta)-\nabla_{\theta}\hat{\ell}^{\alpha}(\theta)\|_{2}^{2}$.
We see, by a similar argument to the variance bound, that

$\displaystyle\mathbb{E}\big{\|}\nabla_{\theta}\hat{s}_{k}^{1}(\theta)-\nabla_{\theta}\hat{s}_{k}^{\alpha}(\theta)\big{\|}_{2}^{2}$
$\displaystyle\leq\;\sum_{n=1}^{N}\mathbb{E}\left\lVert\nabla_{\theta}\Bigg{(}\left|\log\left(\sum_{j=1}^{J}w_{n,j}^{F,\theta,1,k}\right)-\log\left(\sum_{j=1}^{J}w_{n,j}^{F,\theta,\alpha,k}\right)\right|+\right.$
$\displaystyle\hskip
85.35826pt\left.\left|\log\left(\sum_{j=1}^{J}w_{n-1,j}^{A,F,\theta,1,k}\right)+\log\left(\sum_{j=1}^{J}w_{n-1,j}^{A,F,\theta,\alpha,k}\right)\right|\Bigg{)}\right\rVert_{2}^{2}.$
(244)

Each of these terms can be bounded using the derivative of the logarithm and
that $w_{n,j}^{F,\theta,1,k}=w_{n,j}^{F,\theta,\alpha,k}=1$ when
$\theta=\phi$,

$\displaystyle\left\lVert\nabla_{\theta}\log\left(\sum_{j=1}^{J}w_{n,j}^{F,\theta,1,k}\right)-\nabla_{\theta}\log\left(\sum_{j=1}^{J}w_{n,j}^{F,\theta,\alpha,k}\right)\right\rVert_{2}^{2}$
(245)
$\displaystyle\leq\frac{1}{J}\sum_{j=1}^{J}\sum_{i=n-k}^{n}(1-\alpha^{(n-i)})\left\lVert\nabla_{\theta}\log\left(g_{i,j}^{A,F,\theta}\right)\right\rVert_{2}^{2}$
(246)
$\displaystyle\leq\frac{1}{J}\sum_{j=1}^{J}\sum_{i=n-k}^{n}(1-\alpha^{(n-i)})p\,G^{\prime}(\theta)^{2}$
(247) $\displaystyle\leq
p\,G^{\prime}(\theta)^{2}\left(k-\frac{\alpha(1-\alpha^{k})}{1-\alpha}\right),$
(248)

where the second-last line follows from Assumption (A3). So,

$\displaystyle\mathbb{E}\big{\|}\nabla_{\theta}\hat{s}_{k}^{1}(\theta)-\nabla_{\theta}\hat{s}_{k}^{\alpha}(\theta)\big{\|}_{\infty}\leq
2Np\,G^{\prime}(\theta)^{2}\left(k-\frac{\alpha(1-\alpha^{k})}{1-\alpha}\right)\leq
2Np\,G^{\prime}(\theta)^{2}k.$ (249)

Now, we address the third term of Equation 230,
$\mathbb{E}\|\nabla_{\theta}\hat{s}_{k}^{\alpha}(\theta)-\nabla_{\theta}\hat{\ell}^{\alpha}(\theta)\|_{\infty}$.
From the variance argument, we know that we can bound
$\|\nabla_{\theta}\hat{\ell}^{\alpha}(\theta)-\nabla_{\theta}\hat{s}_{k}^{\alpha}(\theta)\|$:

$\displaystyle\left\lVert\nabla_{\theta}\hat{\ell}^{\alpha}(\theta)-\nabla_{\theta}\hat{s}_{k}^{\alpha}(\theta)\right\rVert_{2}^{2}$
(250)
$\displaystyle=\left\lVert\sum_{n=1}^{N}\\!\nabla_{\theta}\log\\!\left(\\!\frac{\sum_{j=1}^{J}\prod_{i=1}^{n}\left(\frac{g_{i,j}^{A,F,\theta}}{g_{i,j}^{A,F,\phi}}\right)^{\alpha^{(n-i)}}}{\sum_{j=1}^{J}\prod_{i=1}^{n-1}\left(\frac{g_{i,j}^{A,F,\theta}}{g_{i,j}^{A,F,\phi}}\right)^{\alpha^{(n-i)}}}\right)\\!-\\!\sum_{n=1}^{N}\\!\nabla_{\theta}\log\\!\left(\\!\frac{\sum_{j=1}^{J}\prod_{i=n-k}^{n}\left(\frac{g_{i,j}^{A,F,\theta}}{g_{i,j}^{A,F,\phi}}\right)^{\alpha^{(n-i)}}}{\sum_{j=1}^{J}\prod_{i=n-k}^{n-1}\left(\frac{g_{i,j}^{A,F,\theta}}{g_{i,j}^{A,F,\phi}}\right)^{\alpha^{(n-i)}}}\right)\\!\right\rVert_{2}^{2}$
(251)
$\displaystyle\leq\sum_{n=1}^{N}\left\lVert\nabla_{\theta}\left(\log\left(\sum_{j=1}^{J}w_{n,j}^{F,\theta,\alpha}\right)-\log\\!\left(\sum_{j=1}^{J}w_{n,j}^{F,\theta,\alpha,k}\right)\right)\right\lVert_{2}^{2}$
$\displaystyle\hskip
113.81102pt+\sum_{n=1}^{N}\left\lVert\nabla_{\theta}\left(\log\left(\sum_{j=1}^{J}w_{n-1,j}^{A,F,\theta,\alpha}\right)+\log\left(\sum_{j=1}^{J}w_{n-1,j}^{A,F,\theta,\alpha,k}\right)\right)\right\rVert_{2}^{2}.$
(252)

Using the derivative of the logarithm and that
$w_{n,j}^{F,\theta,\alpha}=w_{n,j}^{F,\theta,\alpha,k}=1$ when $\theta=\phi$,

$\displaystyle\left\lVert\nabla_{\theta}\log\left(\sum_{j=1}^{J}w_{n,j}^{F,\theta,\alpha}\right)-\nabla_{\theta}\log\left(\sum_{j=1}^{J}w_{n,j}^{F,\theta,\alpha,k}\right)\right\rVert_{2}^{2}$
(253)
$\displaystyle\leq\frac{1}{J}\sum_{j=1}^{J}\sum_{i=1}^{n-k}\alpha^{(n-i)}\left\lVert\nabla_{\theta}\log\left(g_{i,j}^{A,F,\theta}\right)\right\rVert_{2}^{2}$
(254)
$\displaystyle\leq\frac{1}{J}\sum_{j=1}^{J}\sum_{i=1}^{n-k}\alpha^{(n-i)}p\,G^{\prime}(\theta)^{2}$
(255) $\displaystyle\leq
p\,G^{\prime}(\theta)^{2}\frac{\alpha^{k}-\alpha^{n}}{1-\alpha},$ (256)

where the second-last line follows from Assumption (A3). Putting it together
and taking expectations on both sides, we obtain

$\mathbb{E}\left\lVert\nabla_{\theta}\hat{\ell}^{\alpha}(\theta)-\nabla_{\theta}\hat{s}_{k}^{\alpha}(\theta)\right\rVert_{2}^{2}\leq\frac{2\alpha^{k}}{1-\alpha}Np\,G^{\prime}(\theta)^{2},$
(257)

which is our desired error bound. Therefore, our decomposition yields the MSE
bound

$\displaystyle\mathbb{E}\big{\|}\nabla_{\theta}\ell(\theta)-\nabla_{\theta}\hat{\ell}^{\alpha}(\theta)\big{\|}_{2}^{2}$
(258)
$\displaystyle\leq\mathbb{E}\big{\|}\nabla_{\theta}\ell(\theta)-\nabla_{\theta}\hat{s}_{k}^{1}(\theta)\big{\|}_{2}^{2}+\mathbb{E}\big{\|}\nabla_{\theta}\hat{s}_{k}^{1}(\theta)-\nabla_{\theta}\hat{s}_{k}^{\alpha}(\theta)\big{\|}_{2}^{2}+\mathbb{E}\big{\|}\nabla_{\theta}\hat{s}_{k}^{\alpha}(\theta)-\nabla_{\theta}\hat{\ell}^{\alpha}(\theta)\big{\|}_{2}^{2}$
(259) $\displaystyle\lesssim\min_{k\leq
N}\left(N\frac{Cp\,G^{\prime}(\theta)^{2}k^{2}}{J}+Np\,G^{\prime}(\theta)^{2}(1-\epsilon)^{\lfloor
k/(c\log(J))\rfloor}\right.$ $\displaystyle\hskip
113.81102pt\left.+2Np\,G^{\prime}(\theta)^{2}\left(k-\frac{\alpha(1-\alpha^{k})}{1-\alpha}\right)+2\frac{\alpha^{k}}{1-\alpha}Np\,G^{\prime}(\theta)^{2}\right)$
(260) $\displaystyle\lesssim\min_{k\leq
N}Np\,G^{\prime}(\theta)^{2}\left(\frac{k^{2}}{J}+(1-\epsilon)^{\lfloor
k/(c\log(J))\rfloor}+k+\frac{\alpha^{k}+\alpha^{k+1}-\alpha}{1-\alpha}\right).$
(261)

∎

######  Corollary 1 (MSE of MOP-$0$).

The MSE of MOP-$0$, i.e. the estimator of [24], is

$\mathbb{E}\big{\|}\nabla_{\theta}\ell(\theta)-\nabla_{\theta}\hat{\ell}^{0}(\theta)\big{\|}_{2}^{2}\;\lesssim\;Np\,G^{\prime}(\theta)^{2}\left(\frac{1}{J}+(1-\epsilon)^{\lfloor
1/(c\log(J))\rfloor}\right).$ (262)

###### Proof.

Observe that MOP-$0$ is equivalent to MOP-$(1,1)$. Noting that $k=1$ and that
we only need to bound the difference between $\nabla_{\theta}\ell(\theta)$ and
$\nabla_{\theta}\hat{s}_{1}^{1}(\theta)$, repeating the above arguments almost
verbatim allows us to obtain this result.

As before, we bound the first term,
$\sum_{n=1}^{N}\|\nabla_{\theta}\ell_{n}(\theta)-\nabla_{\theta}\hat{s}_{n,1}^{1}(\theta)\|_{2}^{2}$,
by decomposing it into two terms,

$\displaystyle\sum_{n=1}^{N}\big{\|}\nabla_{\theta}\ell_{n}(\theta)-\nabla_{\theta}\hat{s}_{n,1}^{1}(\theta)\big{\|}_{2}^{2}$
$\displaystyle\leq$ (263)
$\displaystyle\sum_{n=1}^{N}\big{\|}\nabla_{\theta}\ell_{n}(\theta)-\nabla_{\theta}\hat{s}_{n,1|\hat{\pi}_{n-1}=\pi_{n-1}}^{1}(\theta)\big{\|}_{2}^{2}+\sum_{n=1}^{N}\big{\|}\nabla_{\theta}\hat{s}_{n,1|\hat{\pi}_{n-1}=\pi_{n-1}}^{1}(\theta)-\nabla_{\theta}\hat{s}_{n,1}^{1}(\theta)\big{\|}_{2}^{2}.$

The first term is a particle approximation dependent on $1$ timesteps, so by
Lemma 2 of [42], this is bounded by

$\mathbb{E}\big{\|}\nabla_{\theta}\ell_{n}(\theta)-\nabla_{\theta}\hat{s}_{n,1|\hat{\pi}_{n-1}=\pi_{n-1}}^{1}(\theta)\big{\|}_{2}^{2}\leq\frac{Cp\,G^{\prime}(\theta)^{2}}{J}.$
(264)

The second term amounts to bounding the difference between functionals of two
different particle measures that mix under the same Markov kernel. Here, we
use Assumptions (A2) and (A3) that ensure strong mixing. We know from Theorem
3 of [42] that when $\textbf{M}_{n,n+k}$ is the $k$-step Markov operator from
timestep $n$ and $\beta_{\text{TV}}(M)=\sup_{x,y\in
E}\|M(x,\cdot)-M(y,\cdot)\|_{\mathrm{TV}}=\sup_{\mu,\nu\in\mathcal{P},\mu\neq\nu}\frac{\|\mu
M-\nu M\|_{\mathrm{TV}}}{\|\mu-\nu\|_{\mathrm{TV}}}$ is the Dobrushin
contraction coefficient of a Markov operator,

$\beta_{\text{TV}}(\textbf{M}_{n,n+k})\leq(1-\epsilon)^{\lfloor
k/(c\log(J))\rfloor},$ (265)

i.e. the mixing time of the particle filter is $O(\log(J))$, where $\epsilon$
and $c$ depend on $\bar{M},\underline{M},\bar{G},\underline{G}$ in (A2) and
(A3).

Then, we can bound
$\mathbb{E}\|\nabla_{\theta}\hat{s}_{n,1|\hat{\pi}_{n-1}=\pi_{n-1}}^{1}(\theta)-\nabla_{\theta}\hat{s}_{n,1}^{1}(\theta)\|_{2}^{2}$
by

$\sup_{\mu,\nu\in\mathcal{P},\mu\neq\nu}\frac{\|\mu\textbf{M}_{n,n+1}-\nu\textbf{M}_{n,n+1}\|_{\mathrm{TV}}}{\|\mu-\nu\|_{\mathrm{TV}}}=\beta_{TV}(\textbf{M}_{n,n+1})\leq(1-\epsilon)^{\lfloor
1/(c\log(J))\rfloor},$ (266)

implying that

$\displaystyle\mathbb{E}\big{\|}\nabla_{\theta}\hat{s}_{n,1|\hat{\pi}_{n-1}=\pi_{n-1}}^{1}(\theta)-\nabla_{\theta}\hat{s}_{n,1}^{1}(\theta)\big{\|}_{2}^{2}$
(267) $\displaystyle\hskip
85.35826pt\lesssim\sup_{\mu,\nu}\sup_{\|\psi\|_{\infty}\leq
1/2}p\,G^{\prime}(\theta)^{2}\big{|}(\mu\textbf{M}_{n,n+1})(\psi)-(\nu\textbf{M}_{n,n+1})(\psi)\big{|}$
(268) $\displaystyle\hskip
85.35826pt\leq\sup_{\mu,\nu}p\,G^{\prime}(\theta)^{2}\big{\|}\mu\textbf{M}_{n,n+1}-\nu\textbf{M}_{n,n+1}\big{\|}_{\mathrm{TV}}$
(269) $\displaystyle\hskip 85.35826pt\leq
p\,G^{\prime}(\theta)^{2}(1-\epsilon)^{\lfloor
1/(c\log(J))\rfloor}\|\hat{\pi}_{n-1}-\pi_{n-1}\|_{\text{TV}}$ (270)
$\displaystyle\hskip 85.35826pt\leq
p\,G^{\prime}(\theta)^{2}(1-\epsilon)^{\lfloor
1/(c\log(J))\rfloor}\sup_{\|\psi\|_{\infty}\leq
1/2}\big{|}\hat{\pi}_{n-1}(\psi)-\pi_{n-1}(\psi)\big{|}$ (271)
$\displaystyle\hskip 85.35826pt\lesssim
p\,G^{\prime}(\theta)^{2}(1-\epsilon)^{\lfloor 1/(c\log(J))\rfloor}.$ (272)

Therefore, we have that

$\displaystyle\mathbb{E}\big{\|}\nabla_{\theta}\ell(\theta)-\nabla_{\theta}\hat{s}_{1}^{1}(\theta)\big{\|}_{2}^{2}$
(273)
$\displaystyle\leq\sum_{n=1}^{N}\mathbb{E}\big{\|}\nabla_{\theta}\ell_{n}(\theta)-\nabla_{\theta}\hat{s}_{n,1|\hat{\pi}_{n-1}=\pi_{n-1}}^{1}(\theta)\big{\|}_{2}^{2}+\sum_{n=1}^{N}\mathbb{E}\big{\|}\nabla_{\theta}\hat{s}_{n,1|\hat{\pi}_{n-1}=\pi_{n-1}}^{1}(\theta)-\nabla_{\theta}\hat{s}_{n,1}^{1}(\theta)\big{\|}_{2}^{2}$
(274) $\displaystyle\lesssim
N\frac{Cp\,G^{\prime}(\theta)^{2}}{J}+Np\,G^{\prime}(\theta)^{2}(1-\epsilon)^{\lfloor
1/(c\log(J))\rfloor}.$ (275)

∎

##  Appendix G Figures for Bayesian Inference

Figure 7: Convergence diagnostics for the Metropolis-Hastings variant particle
MCMC with the random walk proposal and an informative empirical prior from
IF2. Here, we display the results for the trend parameter in the Dhaka cholera
model of [18]. Figure 8: Convergence diagnostics for a No-U-Turn Sampler
(NUTS) with uniform priors on a compact set. Again, we display the results for
the trend parameter in the Dhaka cholera model of [18]. The NUTS sampler
explores more of the posterior than the Metropolis-Hastings sampler, but fails
to converge quickly.

Figure 9: Posterior estimates from NUTS powered by MOP-$\alpha$ with the
informative nonparametric empirical prior obtained from a kernel density
estimator on the IF2 parameter cloud. The posterior estimates from each chain
are largely in agreement.