Change Point Localization and Inference in Dynamic Multilayer Networks

## Abstract

Abstract We study offline change point localization and inference in dynamic multilayer random dot product graphs (D-MRDPGs), where at each time point, a multilayer network is observed with shared node latent positions and time-varying, layer-specific connectivity patterns.
We propose a novel two-stage algorithm that combines seeded binary segmentation with low-rank tensor estimation, and establish its consistency in estimating both the number and locations of change points.
Furthermore, we derive the limiting distributions of the refined estimators under both vanishing and non-vanishing jump regimes. To the best of our knowledge, this is the first result of its kind in the context of dynamic network data.
We also develop a fully data-driven procedure for constructing confidence intervals. Extensive numerical experiments demonstrate the superior performance and practical utility of our methods compared to existing alternatives.

## 1 Introduction

Statistical network analysis studies relationships among entities represented as nodes connected by edges.
While single-layer networks capture pairwise relationships efficiently, many real-world systems involve multiple types of interaction among the same set of nodes.
Multilayer networks address this complexity by organizing these varied interactions into distinct layers over a common node set.
For instance, in social networks, individuals may simultaneously interact through various relationships such as collaboration and friendship, each forming a distinct layer (e.g. Porter, 2018). Modeling multilayer structures enables the capture of heterogeneity and the identification of common structures across various interaction types.

In practice, network structures often evolve over time. For instance, transportation networks may exhibit gradual diurnal variations or sudden structural changes due to unexpected events like accidents or road closures. Detecting these sudden shifts and providing adaptive strategies, such as dynamic traffic signal control or rerouting recommendations, is crucial for efficient transportation management. These abrupt structural shifts are referred to as change points. This naturally falls in the territory of change point analysis.

Change point analysis is a well-established area in statistics concerned with detecting abrupt structural changes in ordered data. It can be broadly classified into online and offline settings, depending on whether data are analyzed sequentially as they are collected or retrospectively after the full dataset has been observed. In the context of dynamic networks, online change point detection has been studied in models such as inhomogeneous Bernoulli networks (e.g., Yu et al., 2021) and random weighted edge networks (e.g., Chen et al., 2024). Offline detection has been explored in various network models, including inhomogeneous Bernoulli networks (e.g., Wang et al., 2021), stochastic block models (e.g., Xu and Lee, 2022; Bhattacharjee et al., 2020) and random dot product graphs (e.g., Padilla et al., 2022). More recently, Wang et al. (2023) investigated online change point detection in dynamic multilayer random dot product graphs (D-MRDPGs).

In this paper, we study offline change point localization and inference for D-MRDPGs. Specifically, at each time point, we observe a realization of an $L$-layered multilayer network, where nodes are associated with fixed but latent positions, and layer-specific weight matrices capture heterogeneous interactions across layers. These weight matrices are allowed to vary over time. Our goal is to develop efficient procedures for localizing and inferring change points under this dynamic multilayer structure in the offline setting.

### 1.1 List of contributions

The main contributions of this paper are summarized as follows.

First, to the best of our knowledge, this is the first work to address offline change point detection in dynamic multilayer networks. We propose a novel two-stage procedure: the first stage generates a coarse set of candidates using seeded binary segmentation, in the spirit of Kovács et al. (2023), combined with CUSUM statistics. The second stage refines these candidates utilizing low-rank tensor estimation techniques. Under appropriate conditions, we establish consistency in both the number and locations of estimated change points.

Second, we derive the limit distributions of the refined change point estimators. Depending on whether the jump size remains fixed or vanishes as the time horizon diverges, the limiting distributions exhibit two distinct regimes. To the best of our knowledge, these are the first such results established in the network literature. We further develop a completely data-driven procedure for constructing confidence intervals for the true change points.

Lastly, we conduct extensive numerical experiments to assess the performance of our proposed methods, demonstrating substantial improvements over existing state-of-the-art algorithms.

### 1.2 Notation and organization

For any positive integer $p$, let $[p]=\{1,\ldots,p\}$.
Let $\{a_{n}\}_{n\in\mathbb{N}^{+}}$ and $\{b_{n}\}_{n\in\mathbb{N}^{+}}$ be sequences of positive real numbers. We write $a_{n}=O(b_{n})$ if
$a_{n}\leq Cb_{n}$ for some constant $C>0$ independent of $n$ and all sufficiently large $n$, and
$a_{n}=\Theta(b_{n})$ if both $a_{n}=O(b_{n})$ and $b_{n}=O(a_{n})$.
For a sequence of random variables $\{X_{n}\}_{n\in\mathbb{N}^{+}}$, we write $X_{n}=O_{p}(a_{n})$ if $\lim_{M\to\infty}\limsup_{n}\mathbb{P}(|X_{n}|\geq Ma_{n})=0$. For any two sets $\mathcal{C}$ and $\mathcal{C}^{\prime}$, define the one-sided Hausdorff distance $d(\mathcal{C}^{\prime}|\mathcal{C})=\max_{c\in\mathcal{C}}\min_{c^{\prime}\in
\mathcal{C}^{\prime}}|c^{\prime}-c|$ with the convention that $d(\mathcal{C}^{\prime}|\mathcal{C})=\infty$ if either
$\mathcal{C}^{\prime}=\emptyset$ or $\mathcal{C}=\emptyset$.

For any matrix $A\in\mathbb{R}^{p_{1}\times p_{2}}$, let $A_{i}$ and $A^{j}$ be the $i$th row and $j$th column of $A$, respectively, and let $\sigma_{1}(A)\geq\dots\geq\sigma_{p_{1}\wedge p_{2}}(A)\geq 0$ be its singular values. For any order-3 tensors $\mathbf{M},\mathbf{Q}\in\mathbb{R}^{p_{1}\times p_{2}\times p_{3}}$, define the inner product
$\langle\mathbf{M},\mathbf{Q}\rangle=\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}\sum_{
l=1}^{p_{3}}\mathbf{M}_{i,j,l}\mathbf{Q}_{i,j,l}$ and the Frobenius norm
$\|\mathbf{M}\|_{\mathrm{F}}^{2}=\langle\mathbf{M},\mathbf{M}\rangle$.
The mode-$1$ matricization of a tensor $\mathbf{M}$ is denoted by $\mathcal{M}_{1}(\mathbf{M})\in\mathbb{R}^{p_{1}\times(p_{2}p_{3})}$
with entries $\mathcal{M}_{1}(\mathbf{M})_{i_{1},(i_{2}-1)p_{3}+i_{3}}=\mathbf{M}_{i_{1},i_{
2},i_{3}}$. The mode-2 and mode-3 matricizations are analogously defined as $\mathcal{M}_{2}(\mathbf{M})\in\mathbb{R}^{p_{2}\times(p_{3}p_{1})}$ and $\mathcal{M}_{3}(\mathbf{M})\in\mathbb{R}^{p_{3}\times(p_{1}p_{2})}$, respectively.
The Tucker ranks $(r_{1},r_{2},r_{3})$ of $\mathbf{M}$ are given by $r_{s}=\mathrm{rank}(\mathcal{M}_{s}(\mathbf{M}))$ for $s\in[3]$. For any $s\in[3]$ and matrix $U_{s}\in\mathbb{R}^{q_{s}\times p_{s}}$, the marginal multiplication operator $\times_{1}$ is defined as
$\mathbf{M}\times_{1}U_{1}=\{\sum_{k=1}^{p_{1}}\mathbf{M}_{k,j,l}(U_{1})_{i,k}
\}_{i\in[q_{1}],\,j\in[p_{2}],\,l\in[p_{3}]}\in\mathbb{R}^{q_{1}\times p_{2}
\times p_{3}}$.
Marginal multiplications $\times_{2}$ and $\times_{3}$ are defined similarly.

The remainder of the paper is organized as follows. Section [2](https://arxiv.org/html/2506.21878v1#S2) introduces the D-MRDPG model and our two-stage change point localization procedure, along with its theoretical guarantees. In Section [3](https://arxiv.org/html/2506.21878v1#S3), we derive the limiting distributions of the refined estimators and propose a data-driven procedure for constructing confidence intervals. Section [4](https://arxiv.org/html/2506.21878v1#S4) presents extensive numerical experiments illustrating the empirical performance of the proposed methods. We conclude with a brief discussion in Section [5](https://arxiv.org/html/2506.21878v1#S5). All proofs and auxiliary results are deferred to the Appendix.

## 2 Change point localization

### 2.1 Problem formulation

We begin with the multilayer random dot product graph (MRDPG) model (Jones and Rubin-Delanchy, 2020), which generalizes the random dot product graph (Young and Scheinerman, 2007) to multilayer networks. Each layer is characterized by a distinct weight matrix, while all layers share a common set of latent positions. We focus on undirected edges, as the directed case is similar and thus omitted.

###### Definition 1 (Multilayer random dot product graphs, MRDPGs ) .

Given a sequence of deterministic matrices $\{W_{(l)}\}_{l=1}^{L}\subset\mathbb{R}^{d\times d}$, let $\{X_{i}\}_{i=1}^{n}\subset\mathbb{R}^{d}$ be fixed vectors satisfying $X_{i}^{\top}W_{(l)}X_{j}\in[0,1]$ for all $i,j\in[n],l\in[L]$.
An adjacency tensor $\mathbf{A}\in\{0,1\}^{n\times n\times L}$ follows an MRDPG if

$$ $\displaystyle\mathbb{P}\{\mathbf{A}\}$ $\displaystyle=\prod_{l=1}^{L}\prod_{1\leq i\leq j\leq n}\mathbf{P}_{i,j,l}^{ \mathbf{A}_{i,j,l}}(1-\mathbf{P}_{i,j,l})^{1-\mathbf{A}_{i,j,l}}$ $\displaystyle=\prod_{l=1}^{L}\prod_{1\leq i\leq j\leq n}\big{(}X_{i}^{\top}W_{ (l)}X_{j}\big{)}^{\mathbf{A}_{i,j,l}}\big{(}1-X_{i}^{\top}W_{(l)}X_{j}\big{)}^ {1-\mathbf{A}_{i,j,l}}.$ $$

We write $\mathbf{A}\sim\mathrm{MRDPG}(\{X_{i}\}_{i=1}^{n},\{W_{(l)}\}_{l\in[L]})$ and denote the probability tensor by $\mathbf{P}\in\mathbb{R}^{n\times n\times L}$.

We now extend this static model to a dynamic setting and introduce a change point framework.

###### Definition 2 (Dynamic multilayer random dot product graphs, D-MRDPGs) .

Let $\{X_{i}\}_{i=1}^{n}\subset\mathbb{R}^{d}$ be latent positions and $\{W_{(l)}(t)\}_{l\in[L],t\in[T]}\subset\mathbb{R}^{d\times d}$ be a weight matrix sequence. A sequence of mutually independent adjacency tensors $\{\mathbf{A}(t)\}_{t\in[T]}$ follows the dynamic MRDPGs if
$\mathbf{A}(t)\sim\mathrm{MRDPG}(\{X_{i}\}_{i=1}^{n},\{W_{(l)}(t)\}_{l\in[L]})$ for $t\in[T]$. We write

$$ $\{\mathbf{A}(t)\}_{t=1}^{T}\sim\mathrm{D}\mbox{-}\mathrm{MRDPGs}(\{X_{i}\}_{i= 1}^{n},\{\{W_{(l)}(t)\}_{l\in[L]}\}_{t=1}^{T}),$ $$

and write $\{\mathbf{P}(t)\}_{t=1}^{T}$ as the corresponding sequence of probability tensors.

###### Model 1 .

Let $\{\mathbf{A}(t)\}_{t\in[T]}\subset\{0,1\}^{n\times n\times L}$ follow $\mathrm{D}\mbox{-}\mathrm{MRDPGs}$ as in Definition [2](https://arxiv.org/html/2506.21878v1#Thmdefinition2).

- $(i)$
Assume that there exist change points $0=\eta_{0}<\eta_{1}<\dots<\eta_{K}<T=\eta_{K+1}$ such that for $t\in[T-1]$, $\{W_{(l)}(t)\}_{l=1}^{L}\neq\{W_{(l)}(t+1)\}_{l=1}^{L}$ if and only if $t\in\{\eta_{k}\}_{k=1}^{K}$.
Let $\Delta=\min_{k\in[K+1]}(\eta_{k}-\eta_{k-1})$ be the minimal spacing between two consecutive change points and assume $\Delta=\Theta(T)$.
- $(ii)$
For each $k\in[K]$, define the $k$-th jump size and normalized jump tensor as
$\kappa_{k}=\|\mathbf{P}(\eta_{k+1})-\mathbf{P}(\eta_{k})\|_{\mathrm{F}}\quad
\mbox{and}\quad\boldsymbol{\mathbf{\Psi}}_{k}=\kappa_{k}^{-1}\{\mathbf{P}(\eta
_{k+1})-\mathbf{P}(\eta_{k})\},$
and let $\kappa=\min_{k\in[K]}\kappa_{k}$ denote the smallest jump magnitude.

Model [1](https://arxiv.org/html/2506.21878v1#Thmmodel1) allows abrupt changes in the layers’ connectivity patterns encoded in the weight matrices, while the latent positions remain unchanged over time. This framework is motivated by a range of practical applications.
For example, in air transportation networks (see Section [4.2](https://arxiv.org/html/2506.21878v1#S4.SS2)), nodes represent airports whose intrinsic attributes, such as geographical location and logistical capacity, are relatively stable. In contrast, airline routing preferences reflected in the weight matrices may change due to factors such as seasonal demand fluctuations, route optimization strategies or policy interventions.

In Model [1](https://arxiv.org/html/2506.21878v1#Thmmodel1)$(i)$, we assume that the minimal spacing $\Delta$ between successive change points is of the same order as the total time horizon $T$, which essentially bounds the number of change points $K$. This assumption can be relaxed, as discussed further in Section [5](https://arxiv.org/html/2506.21878v1#S5). In Model [1](https://arxiv.org/html/2506.21878v1#Thmmodel1)$(ii)$, the magnitude of a change is quantified via the Frobenius norm of the difference between expected adjacency tensors. This metric is sufficiently general to accommodate both dense changes - small but widespread deviations across many layers - and sparse changes - large deviations concentrated in a few layers. Throughout, we allow all model parameters, including the number of nodes $n$, number of layers $L$, latent dimension $d$, jump size $\kappa$ and minimal spacing $\Delta$ to diverge with $T$.

###### Definition 1 (Multilayer random dot product graphs, MRDPGs ) .

###### Definition 2 (Dynamic multilayer random dot product graphs, D-MRDPGs) .

###### Model 1 .

### 2.2 Change point localization algorithm

In this section, we introduce a two-stage procedure for offline change point localization in dynamic multilayer networks, detailed in Algorithm [1](https://arxiv.org/html/2506.21878v1#alg1).
Stage I generates a coarse set of change point candidates using seeded binary segmentation and CUSUM statistics. Stage II refines them via localized scan statistics constructed using a tensor-based low-rank estimation technique. This approach builds on Wang et al. (2021) for single-layer networks and extends it to the multilayer setting.

For Stage I, we begin by defining the seeded intervals (Kovács et al., 2023) and CUSUM statistics (Page, 1954) for dynamic multilayer networks in Definitions [3](https://arxiv.org/html/2506.21878v1#Thmdefinition3) and [4](https://arxiv.org/html/2506.21878v1#Thmdefinition4).

###### Definition 3 (Seeded intervals) .

Let $J=\lceil C_{J}\log_{2}(T)\rceil$ for some sufficiently large absolute constant $C_{J}>0$.
For each $j\in[J]$, define the collection of intervals $\mathcal{J}_{j}$ as

$$ $\mathcal{J}_{j}=\{(\lfloor(i-1)T2^{-j}\rfloor,\lceil(i-1)T2^{-j}+T2^{-j+1} \rceil]\colon i\in[2^{j}-1]\}.$ $$

The full collection of seeded intervals is defined as $\mathcal{J}=\bigcup_{j=1}^{J}\mathcal{J}_{j}$.

###### Definition 4 (CUSUM statistics) .

Given a tensor sequence $\{\mathbf{B}(t)\}_{t\in[T]}$ and any $0\leq s<t<e\leq T$, define the CUSUM statistics as

$$ $\widetilde{\mathbf{B}}^{s,e}(t)=\sum_{u=s+1}^{e}\omega_{s,e}^{t}(u)\mathbf{B}( u),\quad\mbox{where}\quad\omega_{s,e}^{t}(u)=\begin{cases}\sqrt{\frac{e-t}{(e- s)(t-s)}},\quad&\mbox{for }u\in[t]\backslash[s],\\ -\sqrt{\frac{t-s}{(e-s)(e-t)}},\quad&\mbox{for }u\in[e]\backslash[t].\end{cases}$ (1) $$

Stage I implements a modified version of seeded binary segmentation (SBS), a computationally efficient algorithm introduced by Kovács et al. (2023). SBS leverages seeded intervals to construct a multiscale collection of candidate regions for detecting multiple change points. Within each interval, the algorithm computes CUSUM statistics and retains time points where the statistic is maximized and exceeds a predefined threshold, as preliminary change point estimators

We next define the refined scan statistics used in Stage II, based on tensor heteroskedastic principal component analysis (TH-PCA), a low-rank tensor estimation method proposed by Han et al. (2022) and detailed in Algorithm [2](https://arxiv.org/html/2506.21878v1#alg2) in Appendix [B](https://arxiv.org/html/2506.21878v1#A2).

###### Definition 5 (Refined scan statistics) .

Let $\{\mathbf{A}^{\prime}(t)\}_{t\in[T]}$ and $\{\mathbf{B}^{\prime}(t)\}_{t\in[T]}$ be independent sequences generated according to Definition [2](https://arxiv.org/html/2506.21878v1#Thmdefinition2). Given $\{(b_{k},s_{k},e_{k})\}_{k=1}^{\widetilde{K}}$, for any $k\in[\widetilde{K}]$ and $t\in(s_{k},e_{k})$, we define the refined scan statistic as

$$ $\widehat{D}^{s_{k},e_{k}}_{b_{k}}(t)=\big{|}\big{\langle}\widehat{\mathbf{P}}^ {s_{k},e_{k}}(b_{k})/\|\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F} },\widetilde{\mathbf{A}^{\prime}}^{s_{k},e_{k}}(t)\big{\rangle}\big{|},$ $$

where

$$ $\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})=\mathrm{TH}\text{-}\mathrm{PCA}\bigg {(}\widetilde{\mathbf{B^{\prime}}}^{s_{k},e_{k}}(b_{k}),(d,d,m^{s_{k},e_{k}}_{ b_{k}}),\sqrt{\frac{(e_{k}-b_{k})(b_{k}-s_{k})}{e_{k}-s_{k}}},\sqrt{\frac{(e_{ k}-b_{k})(b_{k}-s_{k})}{e_{k}-s_{k}}}\bigg{)}$ $$

with $\mathrm{TH}\mbox{-}\mathrm{PCA}$ detailed in Algorithm [2](https://arxiv.org/html/2506.21878v1#alg2), $\widetilde{\mathbf{B^{\prime}}}^{\cdot,\cdot}(\cdot)$ defined in Definition [4](https://arxiv.org/html/2506.21878v1#Thmdefinition4) and $m^{s,e}_{b_{k}}$ defined in Assumption [1](https://arxiv.org/html/2506.21878v1#Thmassumption1)$(ii)$.

Stage II refines each preliminary change point estimate from Stage I by locating the time point that maximizes a refined scan statistic within a local window around the initial estimate. This step leverages the TH-PCA procedure with an additional truncation step (see Algorithm [2](https://arxiv.org/html/2506.21878v1#alg2)) to more accurately estimate the local expected CUSUM adjacency tensors, leading to provably improved localization accuracy.

The assumption of mutual independence among all four sequences in Algorithm [1](https://arxiv.org/html/2506.21878v1#alg1) is imposed for theoretical convenience. In practice (and in our numerical experiments in Section [4](https://arxiv.org/html/2506.21878v1#S4)), Stage I and Stage II are implemented using the same two split tensor sequences.

Figure: Algorithm 1 Two-stage change point localization for D-MRDPGs

###### Definition 3 (Seeded intervals) .

###### Definition 4 (CUSUM statistics) .

###### Definition 5 (Refined scan statistics) .

### 2.3 Theoretical guarantees

This section establishes the theoretical guarantees of the proposed two-stage change point localization procedure (Algorithm [1](https://arxiv.org/html/2506.21878v1#alg1)).
We begin by justifying the use of low-rank tensor estimation via TH-PCA (Algorithm [2](https://arxiv.org/html/2506.21878v1#alg2)) in Stage II through an analysis of the expected CUSUM-transformed and average adjacency tensors. While the expected averaged adjacency tensors introduced below are not used in this section, they are essential for deriving the limiting distributions in Section [3](https://arxiv.org/html/2506.21878v1#S3).

For any $0\leq s<t<e\leq T$, define the expected CUSUM-transformed and average adjacency tensors as

$$ $\widetilde{\mathbf{P}}^{s,e}(t)=\mathbb{E}\big{\{}\widetilde{\mathbf{B}}^{s,e} (t)\big{\}}\quad\mbox{and}\quad{\mathbf{P}}^{s,e}=\mathbb{E}\big{\{}\mathbf{B} ^{s,e}\big{\}},\quad\mbox{where}\quad\mathbf{B}^{s,e}=(e-s)^{-1}\sum_{t=s+1}^{ e}\mathbf{B}(t),$ (2) $$

and $\widetilde{\mathbf{B}}^{\cdot,\cdot}(\cdot)$ is defined in Definition [4](https://arxiv.org/html/2506.21878v1#Thmdefinition4).
Both tensors admit Tucker representations of the form
$\widetilde{\mathbf{P}}^{s,e}(t)=\mathbf{S}\times_{1}X\times_{2}X\times_{3}
\widetilde{Q}^{s,e}(t),$ and ${\mathbf{P}}^{s,e}=\mathbf{S}\times_{1}X\times_{2}X\times_{3}Q^{s,e}$,
where $X=(X_{1},\ldots,X_{n})^{\top}\in\mathbb{R}^{n\times d}$ and
$\mathbf{S}\in\mathbb{R}^{d\times d\times d^{2}}$ with $\mathbf{S}_{i,j,l}=\mathbbm{1}\{l=(i-1)d+j\}$.
The matrices $\widetilde{Q}^{s,e}(t)$ and $Q^{s,e}$ are given by

$$ $\widetilde{Q}^{s,e}(t)=\sum_{u=s+1}^{e}\omega_{s,e}^{t}(u)Q(u),\quad Q^{s,e}=( e-s)^{-1}\sum_{t=s+1}^{e}Q(t),$ (3) $$

where $\omega_{s,e}^{t}(u)$ is define in ([1](https://arxiv.org/html/2506.21878v1#S2.E1)) and $Q(u)\in\mathbb{R}^{L\times d^{2}}$ with rows

$$ $\big{(}Q(u)\big{)}_{l}=\big{(}(W_{(l)}(u))_{1}\cdots(W_{(l)}(u))_{d}\big{)}, \quad l\in[L].$ (4) $$

To establish the low-rank structure of $\widetilde{\mathbf{P}}^{s,e}(t)$ and $\mathbf{P}^{s,e}$ (in terms of Tucker ranks, see Section [1.2](https://arxiv.org/html/2506.21878v1#S1.SS2)), and to state theoretical guarantees for Algorithm [1](https://arxiv.org/html/2506.21878v1#alg1), we state some necessary assumptions below.

###### Assumption 1 .

Consider $\mathrm{D}\mbox{-}\mathrm{MRDPGs}(\{X_{i}\}_{i=1}^{n},\{\{W_{(l)}(t)\}_{l\in[L
]}\}_{t=1}^{T})$ from Definition [2](https://arxiv.org/html/2506.21878v1#Thmdefinition2).

- $(i)$
Let $X=(X_{1},\ldots,X_{n})^{\top}\in\mathbb{R}^{n\times d}$. Assume that $\mathrm{rank}(X)=d$, $\sigma_{1}(X)/\sigma_{d}(X)\leq C_{\sigma}$ and $\sigma_{d}(X)\geq C_{\mathrm{gap}}\sqrt{n}$ with absolute constants $C_{\sigma},C_{\mathrm{gap}}>0$.
- $(ii)$
For any $0\leq s<t<e\leq T$, let $\widetilde{Q}^{s,e}(t)\in\mathbb{R}^{L\times d^{2}}$ be defined in ([3](https://arxiv.org/html/2506.21878v1#S2.E3)). Denote $m^{s,e}_{t}=\mathrm{rank}(\widetilde{Q}^{s,e}(t))$. Assume that
$\sigma_{1}\big{(}\widetilde{Q}^{s,e}(t)\big{)}/\sigma_{m^{s,e}_{t}}\big{(}
\widetilde{Q}^{s,e}(t)\big{)}\leq C_{\sigma}$ and $\sigma_{m^{s,e}_{t}}\big{(}\widetilde{Q}^{s,e}(t)\big{)}\geq C_{\mathrm{gap}}$ with absolute constants $C_{\mathrm{gap}},C_{\sigma}>0$.
- $(iii)$
For any $0\leq s<e\leq T$, let $Q^{s,e}\in\mathbb{R}^{L\times d^{2}}$ be defined in ([3](https://arxiv.org/html/2506.21878v1#S2.E3)). Denote $m^{s,e}=\mathrm{rank}(Q^{s,e})$. Assume that
$\sigma_{1}\big{(}Q^{s,e}\big{)}/\sigma_{m^{s,e}}\big{(}Q^{s,e}\big{)}\leq C_{\sigma}$ and $\sigma_{m^{s,e}}\big{(}Q^{s,e}\big{)}\geq C_{\mathrm{gap}}$
with absolute constants $C_{\mathrm{gap}},C_{\sigma}>0$.

Assumption [1](https://arxiv.org/html/2506.21878v1#Thmassumption1)$(i)$ imposes a full-rank condition on the latent position matrix $X$, requiring its smallest singular value to be at least of order $\sqrt{n}$, with all singular values of the same order. Since $X$ represents latent positions rather than observed data, the full-rankness of $X$ can be interpreted as a condition on the knowledge of the intrinsic dimension $d$, ensuring that the input dimension to TH-PCA is no smaller than the true latent dimension $d$. Further discussion on rank selection, see Wang et al. (2023).

Assumptions [1](https://arxiv.org/html/2506.21878v1#Thmassumption1)$(ii)$ and $(iii)$ - with $(iii)$ for Section [3](https://arxiv.org/html/2506.21878v1#S3) - impose low-rank conditions on the CUSUM and averaged forms of $\{Q(t)\}_{t=1}^{T}$, where each $Q(t)$ comprises the weight matrices $\{W_{(l)}(t)\}_{l=1}^{L}$.
In Appendix [C](https://arxiv.org/html/2506.21878v1#A3), we show that, with high probability, each working interval $(s_{k},e_{k}]$ or $(\tilde{s}_{k},\tilde{e}_{k}]$ contains exactly one change point $\eta_{k}$, implying $\max\{m^{s_{k},e_{k}}_{t},m^{\tilde{s}_{k},\tilde{e}_{k}}\}\leq\mathrm{rank}(Q
(\eta_{k}))+\mathrm{rank}(Q(\eta_{k+1}))$ for $t\in(s_{k},e_{k})$.
This implicitly constraints the ranks of $\{Q(\eta_{k})\}_{k=1}^{K+1}$.
While this low-rank structure may not directly or transparently reflect the explicit model structure, such ambiguity is common in tensor-based models (e.g. Jing et al., 2021).

The signal-to-noise ratio (SNR) is commonly used to characterize the inherent difficulty of change point detection. We now state the SNR condition required for our theoretical guarantees.

###### Assumption 2 (Signal-to-noise ratio condition) .

Assume that there exists a large enough absolute constant $C_{\mathrm{SNR}}>0$ such that

$$ $\kappa\sqrt{\Delta}\geq C_{\mathrm{SNR}}\log(T)\sqrt{nL^{1/2}+d^{2}m_{\max}+nd +Lm_{\max}},$ $$

where $m_{\max}=\max_{k\in[K+1]}\mathrm{rank}\big{(}Q(\eta_{k})\big{)}$ with $Q(\eta_{k})$ defined in ([4](https://arxiv.org/html/2506.21878v1#S2.E4)).

We compare Assumption [2](https://arxiv.org/html/2506.21878v1#Thmassumption2) to its counterpart in Wang et al. (2021). When the sparsity parameter $\rho=1$, their SNR condition (Assumption 3) becomes $\kappa\sqrt{\Delta}\geq C_{\mathrm{SNR}}\log^{1+\xi}(T)\sqrt{nd}$ for some $\xi>0$. Our assumption is consistent with this and extends it to the multilayer setting by accounting for the additional complexity from multilayers and the low-rank structure of layers’ weight matrices.

###### Theorem 1 .

Let $\{\widetilde{\eta}_{k}\}_{k=1}^{\widetilde{K}}$ be the output of Algorithm [1](https://arxiv.org/html/2506.21878v1#alg1).
Suppose the mutually independent adjacency tensor sequences $\{\mathbf{A}(t)\}_{t\in[T]},\{\mathbf{A}^{\prime}(t)\}_{t\in[T]},\{\mathbf{B}(
t)\}_{t\in[T]},\{\mathbf{B}^{\prime}(t)\}_{t\in[T]}\subset\{0,1\}^{n\times n
\times L}$ are generated according to Definition [2](https://arxiv.org/html/2506.21878v1#Thmdefinition2) and satisfy Model [1](https://arxiv.org/html/2506.21878v1#Thmmodel1), Assumptions [1](https://arxiv.org/html/2506.21878v1#Thmassumption1)$(i)$, $(ii)$ and [2](https://arxiv.org/html/2506.21878v1#Thmassumption2).
Assume the threshold $\tau$ is chosen such that $c_{\tau,1}n\sqrt{L}\log^{3/2}(T)<\tau<c_{\tau,2}\kappa^{2}\Delta$,
where $c_{\tau,1},c_{\tau,2}>0$ are sufficiently large and small absolute constants, respectively. We have that

$$ $\mathbb{P}\Big{\{}\widetilde{K}=K\mbox{ and }|\widetilde{\eta}_{k}-\eta_{k}| \leq\epsilon_{k},\,\forall k\in[K]\Big{\}}\geq 1-CT^{-c},\quad\mbox{where } \epsilon_{k}=C_{\epsilon}\frac{\log(T)}{\kappa_{k}^{2}},$ $$

and $C_{\epsilon},C,c>0$ are absolute constants.

Theorem [1](https://arxiv.org/html/2506.21878v1#Thmtheorem1) implies that, with probability tending to $1$ as $T\to\infty$,
the estimated number of change points satisfies $\widetilde{K}=K$ and the relative localization error vanishes:

$$ $\max_{k\in[K]}\Delta^{-1}|\widetilde{\eta}_{k}-\eta_{k}|\leq C_{\epsilon} \Delta^{-1}\kappa^{-2}\log(T)\to 0$ $$

by Assumption [2](https://arxiv.org/html/2506.21878v1#Thmassumption2). This establishes the consistency of Algorithm [1](https://arxiv.org/html/2506.21878v1#alg1) in both detecting and localizing all change points.
The localization error rate matches the minimax-optimal rate up to a logarithmic factor established for single-layer networks in Wang et al. (2021), showing that our method generalizes these guarantees to the more complex multilayer regime without sacrificing localization accuracy. Moreover, compared to the online framework for D-MRDPGs studied in Wang et al. (2023), which yields a localization error rate of order $\kappa^{-2}(d^{2}m_{\max}+nd+Lm_{\max})\log(\Delta/\alpha)$, where $\alpha$ controls the Type-I error rate, our two-stage procedure achieves a significantly sharper rate in the offline setting.

###### Assumption 1 .

###### Assumption 2 (Signal-to-noise ratio condition) .

###### Theorem 1 .

## 3 Limiting distributions

Inference on change points is generally more challenging than establishing high-probability bounds on localization errors. To address this, we introduce a final refinement step, inspired by approaches such as those in Madrid Padilla et al. (2023); Xue et al. (2024); Xu et al. (2024).

Let $\{\mathbf{A}(t)\}_{t\in[T]}$ and $\{\mathbf{B}(t)\}_{t\in[T]}$ be independent samples as defined in Definition [2](https://arxiv.org/html/2506.21878v1#Thmdefinition2). Let $\{\widetilde{\eta}_{k}\}_{k=1}^{K}$ be the output of Algorithm [1](https://arxiv.org/html/2506.21878v1#alg1) with $0=\widetilde{\eta}_{0}<\widetilde{\eta}_{1}<\cdots<\widetilde{\eta}_{
\widetilde{K}}<\widetilde{\eta}_{\widetilde{K}+1}=T$. For each $k\in[\widetilde{K}]$, define the final estimators as

$$ $\widehat{\eta}_{k}=\operatorname*{argmin}_{\tilde{s}_{k}<t<\tilde{e}_{k}} \mathcal{Q}_{k}(t)=\operatorname*{argmin}_{\tilde{s}_{k}<t<\tilde{e}_{k}}\sum_ {u=\tilde{s}_{k}+1}^{t}\|\mathbf{A}(u)-\widehat{\mathbf{P}}^{\widetilde{\eta}_ {k-1},\widetilde{\eta}_{k}}\|_{\mathrm{F}}^{2}+\sum_{u=t+1}^{\tilde{e}_{k}}\| \mathbf{A}(u)-\widehat{\mathbf{P}}^{\widetilde{\eta}_{k},\widetilde{\eta}_{k+1 }}\|_{\mathrm{F}}^{2},$ (5) $$

where $(\tilde{s}_{k},\tilde{e}_{k}]=((\widetilde{\eta}_{k-1}+\widetilde{\eta}_{k})/2
,(\widetilde{\eta}_{k}+\widetilde{\eta}_{k+1})/2]$ and

$$ $\widehat{\mathbf{P}}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}=\mathrm{TH} \text{-}\mathrm{PCA}(\mathbf{B}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}, (d,d,m^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}),1,0),$ (6) $$

with $\mathrm{TH}\mbox{-}\mathrm{PCA}$ detailed in Algorithm [2](https://arxiv.org/html/2506.21878v1#alg2), $\mathbf{B}^{\cdot,\cdot}$ defined in ([2](https://arxiv.org/html/2506.21878v1#S2.E2)) and $m_{\widetilde{\eta}^{k-1},\widetilde{\eta}_{k}}$ defined in Assumption [1](https://arxiv.org/html/2506.21878v1#Thmassumption1)$(iii)$.

###### Theorem 2 .

Let $\{\mathbf{A}(t)\}_{t\in[T]},\{\mathbf{A}^{\prime}(t)\}_{t\in[T]},\{\mathbf{B}(
t)\}_{t\in[T]},\{\mathbf{B}^{\prime}(t)\}_{t\in[T]}\subset\{0,1\}^{n\times n
\times L}$ be mutually independent adjacency tensor sequences generated according to Definition [2](https://arxiv.org/html/2506.21878v1#Thmdefinition2) and satisfying Model [1](https://arxiv.org/html/2506.21878v1#Thmmodel1), Assumptions [1](https://arxiv.org/html/2506.21878v1#Thmassumption1) and [2](https://arxiv.org/html/2506.21878v1#Thmassumption2). Let $\{\widehat{\eta}_{k}\}_{k=1}^{\widetilde{K}}$ be defined in ([5](https://arxiv.org/html/2506.21878v1#S3.E5)) with $\{\widetilde{\eta}_{k}\}_{k=1}^{\widetilde{K}}$ obtained from Algorithm [1](https://arxiv.org/html/2506.21878v1#alg1), using a threshold $\tau$ satisfying condition stated in Theorem [1](https://arxiv.org/html/2506.21878v1#Thmtheorem1).

For $k\in[K]$, if $\kappa_{k}\to 0$, as $T\to\infty$, then when $T\to\infty$, we have $|\widehat{\eta}_{k}-\eta_{k}|=O_{p}(\kappa_{k}^{-2})$ and

$$ $\kappa_{k}^{2}(\widehat{\eta}_{k}-\eta_{k})\xrightarrow{\mathcal{D}} \operatorname*{argmin}_{r\in\mathbb{R}}\mathcal{P}^{\prime}_{k}(r),\quad\mbox{ where}\quad\mathcal{P}^{\prime}_{k}(r)=\begin{cases}-r+2\sigma_{k,k}\mathbb{B} _{1}(-r),&r<0,\\ 0,&r=0,\\ r+2\sigma_{k,k+1}\mathbb{B}_{2}(r),&r>0,\end{cases}$ $$

for $r\in\mathbb{Z}$. Here, $\mathbb{B}_{1}(r)$ and $\mathbb{B}_{2}(r)$ are independent standard Brownian motions, and for any $k^{\prime}\in\{k,k+1\}$, $\sigma_{k,k^{\prime}}^{2}=\operatorname{Var}\big{(}\langle\boldsymbol{\mathbf{
\Psi}}_{k},\mathbf{E}_{k^{\prime}}(1)\rangle\big{)}$, where $\boldsymbol{\mathbf{\Psi}}_{k}$ is the normalized jump tensor (Model [1](https://arxiv.org/html/2506.21878v1#Thmmodel1)), and
$\mathbf{E}_{k^{\prime}}(t)=\mathbf{A}_{k^{\prime}}(t)-\mathbf{P}(\eta_{k^{
\prime}})$ with $\{\mathbf{A}_{k^{\prime}}(t)\}_{t\in\mathbb{Z}}\overset{\mathrm{i.i.d.}}{\sim}
\mathrm{MRDPG}(\{X_{i}\}_{i=1}^{n},\{W_{(l)}(\eta_{k^{\prime}})\}_{l\in[L]})$.

Theorem [2](https://arxiv.org/html/2506.21878v1#Thmtheorem2) establishes the localization error bounds and limiting distributions for the refined change point estimators in the vanishing jump regime ($\kappa_{k}\to 0$).
In particular, it shows the uniform tightness
$\kappa_{k}^{2}|\widehat{\eta}_{k}-\eta_{k}|=O_{p}(1),$
which improves upon Theorem [1](https://arxiv.org/html/2506.21878v1#Thmtheorem1) by a logarithmic factor and guarantees the existence of limiting distributions.
To the best of our knowledge, Theorem [2](https://arxiv.org/html/2506.21878v1#Thmtheorem2) is the first to derive limiting distributions for change point estimators in network data. These limiting distributions are associated with a two-sided Brownian motion. Results for the non-vanishing jump regime ($\kappa_{k}\to\rho_{k}>0$) are deferred to Appendix [A](https://arxiv.org/html/2506.21878v1#A1).

###### Theorem 2 .

### 3.1 Confidence interval construction

Using Theorem [2](https://arxiv.org/html/2506.21878v1#Thmtheorem2), we construct $(1-\alpha)$ confidence intervals for $\eta_{k}$, $k\in[K]$, in the vanishing regime, for a user-specified confidence level $\alpha\in(0,1)$.

Step 1: Estimate the jump size and normalized jump tensor.
Compute the estimated jump size

$$ $\hat{\kappa}_{k}=\|\widehat{\mathbf{P}}^{\widetilde{\eta}_{k},\widetilde{\eta} _{k+1}}-\widehat{\mathbf{P}}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}\|_{ \mathrm{F}},$ $$

and the estimated normalized jump tensor

$$ $\widehat{\boldsymbol{\mathbf{\Psi}}}_{k}=\hat{\kappa}_{k}^{-1}(\widehat{ \mathbf{P}}^{\widetilde{\eta}_{k},\widetilde{\eta}_{k+1}}-\widehat{\mathbf{P}} ^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}),$ $$

where $\widehat{\mathbf{P}}^{\cdot,\cdot}$ is defined in ([6](https://arxiv.org/html/2506.21878v1#S3.E6)).

Step 2: Estimate the variances.
For each $k^{\prime}\in\{k,k+1\}$, compute

$$ $\hat{\sigma}_{k,k^{\prime}}^{2}=\frac{1}{\widetilde{\eta}_{k^{\prime}}- \widetilde{\eta}_{k^{\prime}-1}-1}\sum_{t=\widetilde{\eta}_{k^{\prime}-1}+1}^{ \widetilde{\eta}_{k^{\prime}}}\big{(}\langle\boldsymbol{\widehat{\mathbf{\Psi} }}_{k},\mathbf{A}(t)-\widehat{\mathbf{P}}^{\widetilde{\eta}_{k^{\prime}-1}, \widetilde{\eta}_{k^{\prime}}}\rangle\big{)}^{2}.$ $$

Step 3: Simulate limiting distributions.
Let $B\in\mathbb{N}^{+}$ and $M\in\mathbb{R}^{+}$. For each $b\in[B]$, let

$$ $\hat{u}^{(b)}_{k}=\operatorname*{argmin}_{r\in(-M,M)}\widehat{\mathcal{P}}^{ \prime}_{k}(r),\quad\mbox{where}\quad\widehat{\mathcal{P}}^{\prime}_{k}(r)= \begin{cases}-r+\frac{2\hat{\sigma}_{k,k}}{\sqrt{T}}\sum_{i=\lceil Tr\rceil}^{ -1}z_{i}^{(b)},&r<0,\\ 0,&r=0,\\ r+\frac{2\hat{\sigma}_{k,k+1}}{\sqrt{T}}\sum_{i=1}^{\lfloor Tr\rfloor}z_{i}^{( b)},&r>0,\end{cases}$ $$

with independent standard Gaussian random variables $\{z_{i}^{(b)}\}_{i=-\lfloor TM\rfloor}^{\lceil TM\rceil}$.

Step 4: Construct the confidence interval. Let $\hat{q}_{\alpha/2},\hat{q}_{1-\alpha/2}$ be empirical quantiles of $\{\hat{u}_{k}^{(b)}\}_{b=1}^{B}$. The $(1-\alpha)$ confidence interval for $\eta_{k}$ is given by

$$ $\bigg{[}\widehat{\eta}_{k}-\frac{\hat{q}_{1-\alpha/2}}{\hat{\kappa}_{k}^{2}} \mathbbm{1}\{\hat{\kappa}_{k}\neq 0\},\,\widehat{\eta}_{k}-\frac{\hat{q}_{ \alpha/2}}{\hat{\kappa}_{k}^{2}}\mathbbm{1}\{\hat{\kappa}_{k}\neq 0\}\bigg{]}.$ $$

The empirical performance of this procedure is evaluated in Section [4.1](https://arxiv.org/html/2506.21878v1#S4.SS1).

## 4 Numerical experiments

We evaluate the proposed method on synthetic and real data in Sections [4.1](https://arxiv.org/html/2506.21878v1#S4.SS1) and [4.2](https://arxiv.org/html/2506.21878v1#S4.SS2), respectively.

### 4.1 Simulation studies

**Table 1: Means of evaluation metrics for networks simulated from Scenarios 1 and 2.**
|  |  | Scenario 1 | Scenario 2 |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $n$ | Method | $|\widehat{K}-K|\downarrow$ | $d(\widehat{\mathcal{C}}|\mathcal{C})\downarrow$ | $d(\mathcal{C}|\widehat{\mathcal{C}})\downarrow$ | $C(\mathcal{G},\mathcal{G^{\prime}})\uparrow$ | $|\widehat{K}-K|\downarrow$ | $d(\widehat{\mathcal{C}}|\mathcal{C})\downarrow$ | $d(\mathcal{C}|\widehat{\mathcal{C}})\downarrow$ | $C(\mathcal{G},\mathcal{G^{\prime}})\uparrow$ |
| $50$ | CPDmrdpg | $0.01$ | $0.00$ | $0.42$ | $99.86\%$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |
| gSeg (nets.) | $1.09$ | Inf | Inf | $52.82\%$ | $1.60$ | Inf | Inf | $67.68\%$ |  |
| kerSeg (nets.) | $0.10$ | $0.00$ | $3.12$ | $99.13\%$ | $0.15$ | $0.00$ | $1.53$ | $99.32\%$ |  |
| gSeg (frob.) | $0.52$ | Inf | Inf | $90.12\%$ | $0.23$ | Inf | Inf | $97.71\%$ |  |
| kerSeg (frob.) | $0.26$ | $0.00$ | $5.76$ | $98.35\%$ | $0.35$ | $0.11$ | $3.43$ | $98.37\%$ |  |
| $100$ | CPDmrdpg | $0.00$ | $0.00$ | $0.00$ | $100\%$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |
| gSeg (nets.) | $1.12$ | Inf | Inf | $52.62\%$ | $1.58$ | Inf | Inf | $69.24\%$ |  |
| kerSeg (nets.) | $0.12$ | $0.00$ | $2.82$ | $99.17\%$ | $0.16$ | $0.00$ | $1.81$ | $99.31\%$ |  |
| gSeg (frob.) | $0.47$ | Inf | Inf | $88.71\%$ | $0.16$ | $0.04$ | $1.65$ | $99.17\%$ |  |
| kerSeg (frob.) | $0.30$ | $0.00$ | $6.07$ | $98.11\%$ | $0.40$ | $0.02$ | $4.42$ | $97.81\%$ |  |

To evaluate the performance of our method (Algorithm [1](https://arxiv.org/html/2506.21878v1#alg1)) for change point detection and localization, we compare it to gSeg (Chen and Zhang, 2015) and kerSeg (Song and Chen, 2024). For the competitors, we consider two input types: networks (nets.) and their layer-wise Frobenius norms (frob.). For gSeg, we construct the similarity graph using the minimum spanning tree and apply the original edge-count scan statistics. For kerSeg, we use the kernel-based scan statistics fGKCP1. For both methods, we set the significance level $\alpha=0.05$. Our proposed method is referred to as CPDmrdpg.
Following Wang et al. (2023), we use relatively large Tucker ranks as inputs to TH-PCA (Algorithm [2](https://arxiv.org/html/2506.21878v1#alg2)) for robustness, setting $r_{1}=r_{2}=15$ and $r_{3}=L$ to compute the refined scan statistics (Definition [5](https://arxiv.org/html/2506.21878v1#Thmdefinition5)).
Based on Theorem [1](https://arxiv.org/html/2506.21878v1#Thmtheorem1), we set the threshold $\tau=c_{\tau,1}n\sqrt{L}\log^{3/2}(T)$ with $c_{\tau,1}=0.1$; justification and sensitivity analysis are provided in Appendix [E.1](https://arxiv.org/html/2506.21878v1#A5.SS1). We also assess the confidence intervals constructed utilizing the procedure in Section [3.1](https://arxiv.org/html/2506.21878v1#S3.SS1), a capability not supported by the competitors. We set $B=500$ and $M=T$ as suggested by Xu et al. (2024).

Performance is quantified using the following metrics:

- $(i)$
Absolute error: $|\widehat{K}-K|$ where $\widehat{K}$ and $K$ denote the numbers of estimated and true change points, respectively;
- $(ii)$
One-sided Hausdorff distances (see Section [1.2](https://arxiv.org/html/2506.21878v1#S1.SS2)):
$d(\widehat{\mathcal{C}}|\mathcal{C})$ and $d(\mathcal{C}|\widehat{\mathcal{C}})$
where $\widehat{\mathcal{C}}$ and $\mathcal{C}$ denote the sets of estimated and true change points, respectively;
- $(iii)$
Time segment coverage:
$C(\mathcal{G},\mathcal{G^{\prime}})=T^{-1}\sum_{\mathcal{A}\in\mathcal{G}}|
\mathcal{A}|\cdot\max_{\mathcal{A^{\prime}}\in\mathcal{G^{\prime}}}|\mathcal{A
}\cap\mathcal{A^{\prime}}|/|\mathcal{A}\cup\mathcal{A^{\prime}}|$
where $\mathcal{G}$ and $\mathcal{G}^{\prime}$ denote the partitions of the time span $[1,T]$ into intervals between consecutive true and estimated change points, respectively.

Throughout, we set the time horizon to $T=200$ and the number of layers to $L=4$, and consider node sizes $n\in\{50,100\}$. Each setting is evaluated over 100 Monte Carlo trials. We consider two network models: the Dirichlet distribution model (DDM) and the multilayer stochastic block model (MSBM), with structural changes specified in each scenario. In the DDM, we generate latent positions

$$ $\{X_{i}\}_{i=1}^{n}\cup\{Y_{i}\}_{i=1}^{n}\overset{\mathrm{i.i.d.}}{\sim}\text {Dirichlet}(\mathbf{1}_{d})$ $$

with $d=5$ and $\mathbf{1}_{d}\in\mathbb{R}^{d}$ denoting the all-one vector. For each time $t$, we sample weight matrices $\{W_{(l)}(t)\}_{l=1}^{L}\subset\mathbb{R}^{d\times d}$ with entries

$$ $(W_{(l)}(t))_{u,v}\sim\text{Uniform}((\rho_{t}L+l)/(4L),(\rho_{t}L+l+1)/(4L)).$ $$

The edge probabilities are given by $\mathbf{P}_{i,j,l}(t)=X_{i}^{\top}W_{(l)}(t)Y_{j}$ and the adjacency entries are sampled as $\mathbf{A}_{i,j,l}(t)\sim\text{Bernoulli}(\mathbf{P}_{i,j,l}(t))$. In the MSBM, the edge probability tensor $\mathbf{P}_{i,j,l}(t)\in[0,1]^{n\times n\times L}$ is defined as $\mathbf{P}_{i,j,l}(t)=p_{1,l}$ if nodes $i,j\in\mathcal{B}_{c}$ for some $c\in[C_{t}]$, and $p_{2,l}$ otherwise, where $\{\mathcal{B}_{c}\}_{c\in[C_{t}]}$ partitions the nodes into $C_{t}$ communities. The connection probabilities are drawn from

$$ $p_{1,l}\sim\text{Uniform}((3L+l-1)/(4L),(3L+l)/(4L))$ $$

and

$$ $p_{2,l}\sim\text{Uniform}((2L+l-1)/(4L),(2L+l)/(4L)).$ $$

The adjacency tensor $\mathbf{A}(t)\in\{0,1\}^{n\times n\times L}$ is then sampled $\mathbf{A}_{i,j,l}(t)\overset{\mathrm{ind.}}{\sim}\text{Bernoulli}(\mathbf{P}_
{i,j,l}(t))$.

Scenario 1. We consider the DDM with $K=2$ change points at $t\in\{70,140\}$, yielding $3$ time segments $\{\mathcal{A}_{i}\}_{i=1}^{3}$. We set $\rho_{t}=2$ for $t\in\mathcal{A}_{1}\cup\mathcal{A}_{3}$, and $\rho_{t}=3$ with reversed layer order for $t\in\mathcal{A}_{2}$.

Scenario 2. We consider the MSBM with $K=5$ change points at $t\in\{20,60,80,160,180\}$, resulting in $6$ time segments $\{\mathcal{A}_{i}\}_{i=1}^{6}$. We let $\{\mathcal{B}_{c}(t)\}_{c\in[C_{t}]}$ be evenly-sized communities and specify the changes as follows: $C_{t}=4$ for $t\in\mathcal{A}_{1}$, $C_{t}=2$ for $t\in\mathcal{A}_{2}$, $C_{t}=4$ for $t\in\mathcal{A}_{3}$, $C_{t}=4$ with reversed layer order for $t\in\mathcal{A}_{4}$, $C_{t}=3$ for $t\in\mathcal{A}_{5}$ and $C_{t}=4$ for $t\in\mathcal{A}_{6}$.

Scenario 3.
We consider the MSBM with $K=3$ change points at $t\in\{50,100,150\}$, yielding $4$ time segments $\{\mathcal{A}_{i}\}_{i=1}^{4}$.
The number of communities is fixed at $C_{t}=3$ but in the first layer, the the community sizes vary across segments $(0.3n,0.4n,0.3n)$ in $\mathcal{A}_{1}\cup\mathcal{A}_{4}$, $(0.4n,0.3n,0.3n)$ in $\mathcal{A}_{2}$ and $(0.5n,0.3n,0.2n)$ in $\mathcal{A}_{3}$. The remaining layers retain equal-sized communities.

Scenario 4.
We consider the MSBM with $K=5$ change points at $t\in\{20,60,80,160,180\}$, resulting in $6$ time segments $\{\mathcal{A}_{i}\}_{i=1}^{6}$. The number of communities is fixed at $C_{t}=4$ with equal-sized partitions, while the connection probabilities vary across segments. Specifically, for $\epsilon=0.1$, we let

$$ $p_{1,l}\sim\text{Uniform}\left(0.5\cdot[0.21+\delta_{t}\cdot\epsilon],0.5\cdot [0.25+\delta_{t}\cdot\epsilon]\right)$ $$

and

$$ $p_{2,l}\sim\text{Uniform}\left(0.21+\delta_{t}\cdot\epsilon,0.25+\delta_{t} \cdot\epsilon\right),$ $$

where $\delta_{t}=0$ for $t\in\mathcal{A}_{1}\cup\mathcal{A}_{5}$, $\delta_{t}=1$ for $t\in\mathcal{A}_{2}\cup\mathcal{A}_{4}\cup\mathcal{A}_{6}$ and $\delta_{t}=2$ for $t\in\mathcal{A}_{3}$.

The changes in Scenarios 1 and 4 follow Model [1](https://arxiv.org/html/2506.21878v1#Thmmodel1), while those in Scenarios 2 and 3 do not, allowing us to assess the robustness of our methods. Table [1](https://arxiv.org/html/2506.21878v1#S4.T1) presents results for Scenarios 1 and 2, and Table [3](https://arxiv.org/html/2506.21878v1#A5.T3) in Appendix [E.1](https://arxiv.org/html/2506.21878v1#A5.SS1) for Scenarios 3 and 4.
Across all scenarios, our method achieves the best overall performance, nearly accurately estimating both the number and locations of change points, and remains robust even when Model [1](https://arxiv.org/html/2506.21878v1#Thmmodel1) is violated.
For gSeg, Frobenius norm (frob.) inputs yield better results than networks (nets.), while kerSeg performs better with networks, benefiting from its high-dimensional kernel-based design. Although both competitors exhibit low Hausdorff distances $d(\widehat{\mathcal{C}}|\mathcal{C})$, their higher reverse distances $d(\mathcal{C}|\widehat{\mathcal{C}})$ and frequent errors in estimating the number of change points suggest they often detect spurious change points.

Table [2](https://arxiv.org/html/2506.21878v1#S4.T2) reports the coverage and average lengths of the confidence intervals constructed via the procedure in Section [3.1](https://arxiv.org/html/2506.21878v1#S3.SS1) for node size $n\in\{100,150\}$. The proposed method generally achieves strong coverage with reasonably narrow intervals. Coverage is lower in Scenario 3, where violations of Model [1](https://arxiv.org/html/2506.21878v1#Thmmodel1) and relatively small, layer-specific changes pose greater challenges. The performance improves with larger $n$ as the change magnitudes $\kappa_{k}$ (Model [1](https://arxiv.org/html/2506.21878v1#Thmmodel1)) increase.

**Table 2: The $95\%$ confidence interval coverage (average length) for change points across all scenarios.**
| $n$ | Scenario 1 | Scenario 2 | Scenario 3 | Scenario 4 |
| --- | --- | --- | --- | --- |
| $100$ | $100\%$ $(0.003)$ | $100\%$ $(0.106)$ | $76.67\%$ $(1.528)$ | $100\%$ $(0.605)$ |
| $150$ | $100\%$ $(0.001)$ | $100\%$ $(0.029)$ | $95.33\%$ $(0.653)$ | $100\%$ $(0.294)$ |

### 4.2 Real data experiments

Our analysis incorporates two real data sets, the worldwide agricultural trade network data set presented here and the U.S. air transport network data set in Appendix [E.2](https://arxiv.org/html/2506.21878v1#A5.SS2).

The worldwide agricultural trade network data are available from Food and Agricultural Organization of the United Nations (2022). The dataset comprises annual multilayer networks from 1986 to 2020 ($T=35$), with nodes representing countries and layers representing agricultural products. A directed edge within a layer indicates the trade relation between two countries of a specific agricultural product. We use the top $L=4$ agricultural products by the trade volume and the $n=75$ most active countries based on import/export volume.
Tuning parameters follow the setup in Section [4.1](https://arxiv.org/html/2506.21878v1#S4.SS1), and our method detects change points in 1991, 1999, 2005, and 2013. Results for competing methods and confidence intervals (Section [3.1](https://arxiv.org/html/2506.21878v1#S3.SS1)) are provided in Appendix [E.2](https://arxiv.org/html/2506.21878v1#A5.SS2).

The 1991 change point aligns with the German reunification and the dissolution of the Soviet Union, both of which triggered major political shifts that significantly affected the trade dynamics. The 1999 change point corresponds to the World Trade Organization’s (WTO) Third Ministerial Conference, a key moment in debates on globalization, particularly regarding agricultural subsidies and tariff reductions, with developing nations demanding fairer trade terms. The 2005 change point marks a WTO agreement to eliminate agricultural export subsidies, promoting greater equity in global markets. Finally, the 2013 change point corresponds to the adoption of the WTO’s Bali Package, the first fully endorsed multilateral agreement, which introduced the Trade Facilitation Agreement and key provisions on food security and tariff quota administration, significantly impacting agricultural trade.

## 5 Conclusion

In this paper, we study offline change point localization and inference in dynamic multilayer networks — a setting that, to the best of our knowledge, has not has not been previously addressed.
We propose a two-stage algorithm and establish its consistency in estimating both the number and locations of change points.
We further develop local refinement procedures, derive the limiting distributions and introduce a data-driven method for constructing confidence intervals for the true change points.

Several limitations of this work remain open for future research. First, the assumption $\Delta=\Theta(T)$ precludes frequent change points. This could be relaxed using alternative selection strategies such as the narrowest-over-threshold approach (Baranowski et al., 2019) instead of greedy selection in this paper. Second, our inference procedure is limited to vanishing jumps. It would be interesting to explore practical procedures for the non-vanishing regime, potentially building on bootstrap methods (e.g. Cho and Kirch, 2022). Lastly, the current framework assumes temporal independence; a natural extension would be to incorporate temporal dependence structures (e.g. Padilla et al., 2022; Cho and Owens, 2023).

## References

- Baranowski et al. (2019)
Rafal Baranowski, Yining Chen, and Piotr Fryzlewicz.
Narrowest-over-threshold detection of multiple change points and change-point-like features.
*Journal of the Royal Statistical Society Series B: Statistical Methodology*, 81(3):649–672, May 2019.
ISSN 1467-9868.
doi: 10.1111/rssb.12322.
URL [http://dx.doi.org/10.1111/rssb.12322](http://dx.doi.org/10.1111/rssb.12322).
- Bhattacharjee et al. (2020)
Monika Bhattacharjee, Moulinath Banerjee, and George Michailidis.
Change point estimation in a dynamic stochastic block model.
*Journal of machine learning research*, 21(107):1–59, 2020.
- Bureau of Transportation Statistics (2022)
Bureau of Transportation Statistics.
T-100 domestic market (u.s. carriers), 2022.
URL [https://www.transtats.bts.gov/DL_SelectFields.aspx?gnoyr_VQ=GDL&QO_fu146_anzr=Nv4](https://www.transtats.bts.gov/DL_SelectFields.aspx?gnoyr_VQ=GDL&QO_fu146_anzr=Nv4).
- Chen et al. (2024)
Cathy Yi-Hsuan Chen, Yarema Okhrin, and Tengyao Wang.
Monitoring network changes in social media.
*Journal of Business & Economic Statistics*, 42(2):391–406, 2024.
- Chen and Zhang (2015)
Hao Chen and Nancy Zhang.
Graph-based change-point detection.
*The Annals of Statistics*, 43(1):139 – 176, 2015.
doi: 10.1214/14-AOS1269.
URL [https://doi.org/10.1214/14-AOS1269](https://doi.org/10.1214/14-AOS1269).
- Cho and Kirch (2022)
Haeran Cho and Claudia Kirch.
Bootstrap confidence intervals for multiple change points based on moving sum procedures.
*Computational Statistics & Data Analysis*, 175:107552, 2022.
- Cho and Owens (2023)
Haeran Cho and Dom Owens.
High-dimensional data segmentation in regression settings permitting temporal dependence and non-gaussianity, 2023.
URL [https://arxiv.org/abs/2209.08892](https://arxiv.org/abs/2209.08892).
- Food and Agricultural Organization of the United Nations (2022)
Food and Agricultural Organization of the United Nations.
Food and agriculture data, 2022.
URL [https://www.fao.org/faostat/en/#data/TM](https://www.fao.org/faostat/en/#data/TM).
- Han et al. (2022)
Rungang Han, Rebecca Willett, and Anru R Zhang.
An optimal statistical and computational framework for generalized tensor estimation.
*The Annals of Statistics*, 50(1):1–29, 2022.
- Jing et al. (2021)
Bing-Yi Jing, Ting Li, Zhongyuan Lyu, and Dong Xia.
Community detection on mixture multilayer networks via regularized tensor decomposition.
*The Annals of Statistics*, 49(6):3181–3205, 2021.
- Jones and Rubin-Delanchy (2020)
Andrew Jones and Patrick Rubin-Delanchy.
The multilayer random dot product graph.
*arXiv preprint arXiv:2007.10455*, 2020.
- Kovács et al. (2023)
Solt Kovács, Peter Bühlmann, Housen Li, and Axel Munk.
Seeded binary segmentation: a general methodology for fast and optimal changepoint detection.
*Biometrika*, 110(1):249–256, 2023.
- Madrid Padilla et al. (2022)
Carlos Misael Madrid Padilla, Daren Wang, Zifeng Zhao, and Yi Yu.
Change-point detection for sparse and dense functional data in general dimensions.
*Advances in Neural Information Processing Systems*, 35:37121–37133, 2022.
- Madrid Padilla et al. (2023)
Carlos Misael Madrid Padilla, Haotian Xu, Daren Wang, OSCAR HERNAN MADRID PADILLA, and Yi Yu.
Change point detection and inference in multivariate non-parametric models under mixing conditions.
*Advances in Neural Information Processing Systems*, 36:21081–21134, 2023.
- Padilla et al. (2022)
Oscar Hernan Madrid Padilla, Yi Yu, and Carey E Priebe.
Change point localization in dependent dynamic nonparametric random dot product graphs.
*Journal of Machine Learning Research*, 23(234):1–59, 2022.
- Page (1954)
Ewan S Page.
Continuous inspection schemes.
*Biometrika*, 41(1/2):100–115, 1954.
- Porter (2018)
Mason A. Porter.
What is… a multilayer network?
*Notices of the American Mathematical Society*, 2018.
URL [https://api.semanticscholar.org/CorpusID:53702798](https://api.semanticscholar.org/CorpusID:53702798).
- Song and Chen (2024)
Hoseung Song and Hao Chen.
Practical and powerful kernel-based change-point detection.
*IEEE Transactions on Signal Processing*, 2024.
- Vershynin (2018)
Roman Vershynin.
*High-dimensional probability: An introduction with applications in data science*, volume 47.
Cambridge university press, 2018.
- Wang et al. (2017)
Daren Wang, Yi Yu, and Alessandro Rinaldo.
Optimal covariance change point localization in high dimension.
*arXiv preprint arXiv:1712.09912*, 2017.
- Wang et al. (2021)
Daren Wang, Yi Yu, and Alessandro Rinaldo.
Optimal change point detection and localization in sparse dynamic networks.
*The Annals of Statistics*, 49(1):203–232, 2021.
- Wang et al. (2023)
Fan Wang, Wanshan Li, Oscar Hernan Madrid Padilla, Yi Yu, and Alessandro Rinaldo.
Multilayer random dot product graphs: Estimation and online change point detection.
*arXiv preprint arXiv:2306.15286*, 2023.
- Wellner et al. (2013)
Jon Wellner et al.
*Weak convergence and empirical processes: with applications to statistics*.
Springer Science & Business Media, 2013.
- Xu and Lee (2022)
Cong Xu and Thomas CM Lee.
Statistical consistency for change point detection and community estimation in time-evolving dynamic networks.
*IEEE Transactions on Signal and Information Processing over Networks*, 8:215–227, 2022.
- Xu et al. (2024)
Haotian Xu, Daren Wang, Zifeng Zhao, and Yi Yu.
Change-point inference in high-dimensional regression models under temporal dependence.
*The Annals of Statistics*, 52(3):999–1026, 2024.
- Xue et al. (2024)
Gengyu Xue, Haotian Xu, and Yi Yu.
Change point localisation and inference in fragmented functional data.
*arXiv preprint arXiv:2405.05730*, 2024.
- Young and Scheinerman (2007)
Stephen J Young and Edward R Scheinerman.
Random dot product graph models for social networks.
In *International Workshop on Algorithms and Models for the Web-Graph*, pages 138–149. Springer, 2007.
- Yu et al. (2021)
Yi Yu, Oscar Hernan Madrid Padilla, Daren Wang, and Alessandro Rinaldo.
Optimal network online change point localisation.
*arXiv preprint arXiv:2101.05477*, 2021.
- Zhang et al. (2022)
Anru R Zhang, T Tony Cai, and Yihong Wu.
Heteroskedastic pca: Algorithm, optimality, and applications.
*The Annals of Statistics*, 50(1):53–80, 2022.

## Appendix

All technical details of this paper are deferred to the Appendix. Limiting distributions for the non-vanishing regime are established in Appendix [A](https://arxiv.org/html/2506.21878v1#A1). Additional algorithms used in our procedures are provided in Appendix [B](https://arxiv.org/html/2506.21878v1#A2).
The proof of Theorem [1](https://arxiv.org/html/2506.21878v1#Thmtheorem1) is given in Appendix [C](https://arxiv.org/html/2506.21878v1#A3), while the proofs of the limiting distribution results, including Theorem [2](https://arxiv.org/html/2506.21878v1#Thmtheorem2) in the main text,
and Theorem [3](https://arxiv.org/html/2506.21878v1#Thmtheorem3) in Appendix [A](https://arxiv.org/html/2506.21878v1#A1) are presented in Appendix [D](https://arxiv.org/html/2506.21878v1#A4). Further details and results for Section [4](https://arxiv.org/html/2506.21878v1#S4) are collected in Appendix [E](https://arxiv.org/html/2506.21878v1#A5).

## Appendix A Limiting distributions in the non-vanishing regime

###### Theorem 3 .

Let $\{\mathbf{A}(t)\}_{t\in[T]},\{\mathbf{A}^{\prime}(t)\}_{t\in[T]},\{\mathbf{B}(
t)\}_{t\in[T]},\{\mathbf{B}^{\prime}(t)\}_{t\in[T]}\subset\{0,1\}^{n\times n
\times L}$ be mutually independent adjacency tensor sequences generated according to Definition [2](https://arxiv.org/html/2506.21878v1#Thmdefinition2) and satisfying Model [1](https://arxiv.org/html/2506.21878v1#Thmmodel1), Assumptions [1](https://arxiv.org/html/2506.21878v1#Thmassumption1) and [2](https://arxiv.org/html/2506.21878v1#Thmassumption2). Let $\{\widehat{\eta}_{k}\}_{k=1}^{\widetilde{K}}$ be defined in ([5](https://arxiv.org/html/2506.21878v1#S3.E5)) with $\{\widetilde{\eta}_{k}\}_{k=1}^{\widetilde{K}}$ obtained from Algorithm [1](https://arxiv.org/html/2506.21878v1#alg1), using a threshold $\tau$ satisfying condition stated in Theorem [1](https://arxiv.org/html/2506.21878v1#Thmtheorem1).

For $k\in[K]$, if $\kappa_{k}\to\rho_{k}$, as $T\to\infty$, with $\rho_{k}>0$ being an absolute constant, then when $T\to\infty$, we have $|\widehat{\eta}_{k}-\eta_{k}|=O_{p}(1)$ and

$$ $\widehat{\eta}_{k}-\eta_{k}\xrightarrow{\mathcal{D}}\operatorname*{argmin}_{r \in\mathbb{Z}}\mathcal{P}_{k}(r),\quad\mbox{where}\quad\mathcal{P}_{k}(r)= \begin{cases}-r\rho_{k}^{2}-2\rho_{k}\sum_{t=r+1}^{0}\langle\boldsymbol{ \mathbf{\Psi}}_{k},\mathbf{E}_{k}(t)\rangle,&r<0,\\ 0,&r=0,\\ r\rho_{k}^{2}+2\rho_{k}\sum_{t=1}^{r}\langle\boldsymbol{\mathbf{\Psi}}_{k}, \mathbf{E}_{k+1}(t)\rangle,&r>0,\end{cases}$ $$

for $r\in\mathbb{Z}$. Here, the normalized jump tensor
$\boldsymbol{\mathbf{\Psi}}_{k}$ is defined in Model [1](https://arxiv.org/html/2506.21878v1#Thmmodel1), and for any $k\in[K+1]$ and $t\in\mathbb{Z}$,
$\mathbf{E}_{k}(t)=\mathbf{A}_{k}(t)-\mathbf{P}(\eta_{k})$ with $\{\mathbf{A}_{k}(t)\}_{t\in\mathbb{Z}}\overset{\mathrm{i.i.d.}}{\sim}\mathrm{
MRDPG}(\{X_{i}\}_{i=1}^{n},\{W_{(l)}(\eta_{k})\}_{l\in[L]})$.

The proof of Theorem [3](https://arxiv.org/html/2506.21878v1#Thmtheorem3) is given in Appendix [D](https://arxiv.org/html/2506.21878v1#A4).

Similar to Theorem [2](https://arxiv.org/html/2506.21878v1#Thmtheorem2), Theorem [3](https://arxiv.org/html/2506.21878v1#Thmtheorem3) establishes the uniform tightness $\kappa_{k}^{2}|\widehat{\eta}_{k}-\eta_{k}|=O_{p}(1)$ and further derives the limiting distributions of the refined change point estimators defined in ([5](https://arxiv.org/html/2506.21878v1#S3.E5)), which are associated with a two-sided random walk.

###### Theorem 3 .

## Appendix B Additional algorithms

We present the tensor heteroskedastic principal component analysis (TH-PCA) algorithm introduced in Han et al. (2022), incorporating an additional truncation step, in Algorithm [2](https://arxiv.org/html/2506.21878v1#alg2). Its subroutine,
the heteroskedastic principal component analysis (H-PCA) algorithm proposed by Zhang et al. (2022), is provided in Algorithm [3](https://arxiv.org/html/2506.21878v1#alg3).

Figure: Algorithm 2 Tensor heteroskedastic principal component analysis, TH-PCA$(\mathbf{A},(r_{1},r_{2},r_{3}),\tau_{1},\tau_{2})$

Figure: Algorithm 3 Heteroskedastic principal component analysis, H-PCA$(\Sigma,r)$

## Appendix C Proof of Theorem 1

The proof of Theorem [1](https://arxiv.org/html/2506.21878v1#Thmtheorem1) is in Appendix [C.1](https://arxiv.org/html/2506.21878v1#A3.SS1) with all necessary auxiliary results in Appendix [C.2](https://arxiv.org/html/2506.21878v1#A3.SS2).

### C.1 Proof of Theorem 1

###### Proof.

We first define the event

$$ $\mathcal{A}=\Big{\{}\widetilde{K}=K\mbox{ and }|b_{k}-\eta_{k}|\leq\tilde{ \epsilon},\,\forall k\in[K]\Big{\}},\quad\mbox{where }\tilde{\epsilon}=C_{ \tilde{\epsilon}}\log(T)\bigg{\{}\frac{n\sqrt{L}\log^{1/2}(T)}{\kappa^{2}}+ \frac{\sqrt{\Delta}}{\kappa}\bigg{\}},$ $$

where $\{b_{k}\}_{k=1}^{\widetilde{K}}$ are preliminary change point estimates obtained from Stage I in Algorithm [1](https://arxiv.org/html/2506.21878v1#alg1).
Then by Proposition [4](https://arxiv.org/html/2506.21878v1#Thmtheorem4), it holds that

$$ $\mathbb{P}\{\mathcal{A}\}\geq 1-C_{0}T^{-c_{0}}$ $$

and $C_{\tilde{\epsilon}},C_{0},c_{0}>0$ are absolute constants. Since $\{\mathbf{A}^{\prime}(t)\}_{t=1}^{T}\cup\{\mathbf{B}^{\prime}(t)\}_{t=1}^{T}$ are independent of $\{\mathbf{A}(t)\}_{t=1}^{T}\cup\{\mathbf{B}(t)\}_{t=1}^{T}$, the distribution of $\{\mathbf{A}^{\prime}(t)\}_{t=1}^{T}\cup\{\mathbf{B}^{\prime}(t)\}_{t=1}^{T}$ remains unaffected under the conditioning on the event $\mathcal{A}$. All subsequent analysis in this proof is carried out under the event $\mathcal{A}$. Consequently, we can derive that

$$ $|b_{k}-\eta_{k}|\leq\tilde{\epsilon}\leq\Delta/6,\quad\forall k\in[K],$ (7) $$

where the last inequality follows from Assumption [2](https://arxiv.org/html/2506.21878v1#Thmassumption2) and the fact that $C_{\mathrm{SNR}}$ is a sufficiently large constant.

Step 1. We first establish that for any $k\in[K]$, each working interval $(s_{k},e_{k})$ contains exactly one true change point, namely $\eta_{k}$, and the two endpoints are well separated.

From ([7](https://arxiv.org/html/2506.21878v1#A3.E7)), we obtain that
$\eta_{k}\in[b_{k-1},b_{k+1}]$,

$$ $\eta_{k}-b_{k-1}\geq\eta_{k}-\eta_{k-1}-|\eta_{k-1}-b_{k-1}|\geq\Delta-\Delta/ 6\geq 5\Delta/6,$ $$

and

$$ $b_{k+1}-\eta_{k}\geq\eta_{k+1}-\eta_{k}-|\eta_{k+1}-b_{k+1}|\geq\Delta-\Delta/ 6\geq 5\Delta/6.$ $$

Similarly, we can derive that

$$ $\min\{b_{k}-b_{k-1},\,b_{k+1}-b_{k}\}\geq 2\Delta/3.$ $$

As a result, the working interval

$$ $(s_{k},e_{k}]=\big{(}b_{k-1}+(b_{k}-b_{k-1})/2,\;b_{k+1}-(b_{k+1}-b_{k})/2\big {]},$ $$

contains exactly one change point $\eta_{k}$.
For any $t\in(s_{k},e_{k})$, denote $m^{s_{k},e_{k}}_{t}=\mathrm{rank}(\widetilde{Q}^{s_{k},e_{k}}(t))$ with $\widetilde{Q}^{s_{k},e_{k}}(t)$ defined in ([3](https://arxiv.org/html/2506.21878v1#S2.E3)). Then we have that

$$ $m^{s_{k},e_{k}}_{t}\leq m_{k}+m_{k+1}\leq 2m_{\max},$ (8) $$

where $m_{k}=\mathrm{rank}\big{(}Q(\eta_{k})\big{)}$ with $Q(\eta_{k})$ defined in ([4](https://arxiv.org/html/2506.21878v1#S2.E4)),
and $m_{\max}=\max_{k\in[K+1]}m_{k}$.

In addition, we have that

$$ $b_{k}-s_{k}=(b_{k}-b_{k-1})/2\geq\Delta/3,$ $$

and

$$ $e_{k}-b_{k}=(b_{k+1}-b_{k})/2\geq\Delta/3.$ $$

Therefore,

$$ $\min\{e_{k}-b_{k},b_{k}-s_{k}\}\geq\Delta/3.$ (9) $$

Step 2. We now show that the population statistics $\widetilde{\mathbf{P}}^{s_{k},e_{k}}(b_{k})$ provide a sufficiently strong signal within each working interval $(s_{k},e_{k}]$.

By Lemma [6](https://arxiv.org/html/2506.21878v1#Thmtheorem6), it holds that

$$ $\big{\|}\widetilde{\mathbf{P}}^{s_{k},e_{k}}(t)\big{\|}_{\mathrm{F}}^{2}= \begin{cases}\frac{t-s_{k}}{(e_{k}-s_{k})(e_{k}-t)}(e_{k}-\eta_{k})^{2}\kappa_ {k}^{2},&s_{k}<t\leq\eta_{k},\\ \frac{e_{k}-t}{(e_{k}-s_{k})(t-s_{k})}(\eta_{k}-s_{k})^{2}\kappa_{k}^{2},&\eta _{k}<t<e_{k}.\end{cases}$ $$

Define the scaling factor

$$ $\widetilde{\Delta}_{k}=\sqrt{\frac{(b_{k}-s_{k})\,(e_{k}-b_{k})}{e_{k}-s_{k}}}.$ $$

Assuming without loss of generality that $b_{k}\leq\eta_{k}$ and using ([9](https://arxiv.org/html/2506.21878v1#A3.E9)), we obtain

$$ $\widetilde{\Delta}_{k}^{2}\geq\frac{\min\{b_{k}-s_{k},\,e_{k}-b_{k}\}}{2}\geq \frac{\Delta}{6}.$ (10) $$

Thus, we have that

$$ $\displaystyle\big{\|}\widetilde{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\big{\|}_{ \mathrm{F}}^{2}=$ $\displaystyle\frac{b_{k}-s_{k}}{(e_{k}-s_{k})(e_{k}-b_{k})}(e_{k}-\eta_{k})^{2 }\ \kappa_{k}^{2}=\widetilde{\Delta}_{k}^{2}\bigg{(}\frac{e_{k}-\eta_{k}}{e_{k }-b_{k}}\bigg{)}^{2}\kappa_{k}^{2}$ $\displaystyle=$ $\displaystyle\widetilde{\Delta}_{k}^{2}\bigg{(}1-\frac{\eta_{k}-b_{k}}{e_{k}-b _{k}}\bigg{)}^{2}\kappa_{k}^{2}\geq\frac{\Delta}{6}(1-\frac{\Delta/6}{\Delta/3 })^{2}\kappa_{k}^{2}=\Delta\kappa_{k}^{2}/24,$ (11) $$

where the first inequality follows from ([7](https://arxiv.org/html/2506.21878v1#A3.E7)), ([9](https://arxiv.org/html/2506.21878v1#A3.E9)) and ([10](https://arxiv.org/html/2506.21878v1#A3.E10)).

Step 3.
Note that each entry of the tensor $\widetilde{\mathbf{B}^{\prime}}^{s,e}(b_{k})$ is independently $c_{\sigma}$-sub-Gaussian distributed with mean tensor $\mathbb{E}\{\widetilde{\mathbf{B}^{\prime}}^{s,e}(b_{k})\}=\widetilde{\mathbf{
P}}^{s_{k},e_{k}}(\eta_{k})$ and an absolute constant $c_{\sigma}>0$. By Theorem 4 and Lemma 5 in Wang et al. (2023), and Assumption [1](https://arxiv.org/html/2506.21878v1#Thmassumption1) $(i)$ and $(ii)$, for any $t\in(s,e)$, it holds that

$$ $\mathbb{P}\bigg{\{}\big{\|}\widehat{\mathbf{P}}^{s_{k},e_{k}}(t)-\widetilde{ \mathbf{P}}^{s_{k},e_{k}}(t)\big{\|}_{F}\leq C_{1}\sqrt{(d^{2}m^{s,e}_{t}+nd+ Lm^{s,e}_{t})\log(T)}\bigg{\}}\geq 1-T^{-c_{1}},$ $$

for some constants $C_{1}>0$ and $c_{1}>3$. By ([8](https://arxiv.org/html/2506.21878v1#A3.E8)), we can derive that

$$ $\mathbb{P}\bigg{\{}\big{\|}\widehat{\mathbf{P}}^{s_{k},e_{k}}(t)-\widetilde{ \mathbf{P}}^{s_{k},e_{k}}(t)\big{\|}_{F}\leq C_{2}\sqrt{(d^{2}m_{\max}+nd+Lm_{ \max})\log(T)}\bigg{\}}\geq 1-T^{-c_{1}},$ $$

where $C_{2}>0$ is a constant. Define the event

$$ $\displaystyle\mathcal{B}=\Bigg{\{}\sup_{\begin{subarray}{c}0\leq s_{k}<t<e_{k} \leq T\\ (s_{k},e_{k})\mbox{ contains only one change point }\eta_{k}\end{subarray}}$ $\displaystyle\big{\|}\widehat{\mathbf{P}}^{s_{k},e_{k}}(t)-\widetilde{\mathbf{ P}}^{s_{k},e_{k}}(t)\big{\|}_{F}\leq$ $\displaystyle C_{\mathcal{B}}\sqrt{(d^{2}m_{\max}+nd+Lm_{\max})\log(T)}\Bigg{ \}}.$ (12) $$

with a constant $C_{\mathcal{B}}>0$.
By the union bound argument, it holds that

$$ $\mathbb{P}\big{\{}\mathcal{B}\big{\}}\geq 1-T^{-c_{\mathcal{B}}},$ $$

with a constant $c_{\mathcal{B}}>0$.
By the event $\mathcal{B}$ and the triangle equality, we have that

$$ $\big{\|}\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\big{\|}_{\mathrm{F}}\geq\big {\|}\widetilde{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\big{\|}_{\mathrm{F}}-C_{ \mathcal{B}}\sqrt{(d^{2}m_{\max}+nd+Lm_{\max})\log(T)}\geq\kappa_{k}\sqrt{ \Delta}/48,$ (13) $$

where the last inequality follows from ([C.1](https://arxiv.org/html/2506.21878v1#A3.Ex34)), Assumption [2](https://arxiv.org/html/2506.21878v1#Thmassumption2) and the fact that $C_{\mathrm{SNR}}$ is a sufficiently large constant.
As a consequence,

$$ $\displaystyle 2\Big{\langle}\widetilde{\mathbf{P}}^{s_{k},e_{k}}(b_{k})/\| \widetilde{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F}},\widehat{\mathbf{P} }^{s_{k},e_{k}}(b_{k})/\|\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{ F}}\Big{\rangle}$ $\displaystyle=$ $\displaystyle 2-\Bigg{\|}\frac{\widetilde{\mathbf{P}}^{s_{k},e_{k}}(b_{k})}{\| \widetilde{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F}}}-\frac{\widehat{ \mathbf{P}}^{s_{k},e_{k}}(b_{k})}{\|\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k}) \|_{\mathrm{F}}}\Bigg{\|}_{\mathrm{F}}^{2}$ $\displaystyle=$ $\displaystyle 2-\Bigg{\|}\frac{\big{(}\widetilde{\mathbf{P}}^{s_{k},e_{k}}(b_{ k})-\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\big{)}\|\widetilde{\mathbf{P}}^{ s_{k},e_{k}}(b_{k})\|_{\mathrm{F}}}{\|\widetilde{\mathbf{P}}^{s_{k},e_{k}}(b_{ k})\|_{\mathrm{F}}\|\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F}}}$ $\displaystyle\hskip 14.22636pt+\frac{\widetilde{\mathbf{P}}^{s_{k},e_{k}}(b_{k })\big{(}\|\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F}}-\| \widetilde{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F}}\big{)}}{\| \widetilde{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F}}\|\widehat{\mathbf{P }}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F}}}\Bigg{\|}_{\mathrm{F}}^{2}$ $\displaystyle\geq$ $\displaystyle 2-\Bigg{(}\frac{\|\widetilde{\mathbf{P}}^{s_{k},e_{k}}(b_{k})- \widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F}}\|\widetilde{\mathbf{P }}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F}}}{\|\widetilde{\mathbf{P}}^{s_{k},e_{k}} (b_{k})\|_{\mathrm{F}}\|\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F }}}$ $\displaystyle\hskip 14.22636pt+\frac{\|\widetilde{\mathbf{P}}^{s_{k},e_{k}}(b_ {k})\|_{\mathrm{F}}\|\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})-\widetilde{ \mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F}}}{\|\widetilde{\mathbf{P}}^{s_{ k},e_{k}}(b_{k})\|_{\mathrm{F}}\|\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{ \mathrm{F}}}\Bigg{)}^{2}$ $\displaystyle=$ $\displaystyle 2-\Bigg{(}\frac{\big{(}\|\widetilde{\mathbf{P}}^{s_{k},e_{k}}(b_ {k})\|_{\mathrm{F}}+\|\widetilde{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F }}\big{)}\|\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})-\widetilde{\mathbf{P}}^{s _{k},e_{k}}(b_{k})\|_{\mathrm{F}}}{\|\widetilde{\mathbf{P}}^{s_{k},e_{k}}(b_{k })\|_{\mathrm{F}}\|\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F}}} \Bigg{)}^{2}$ $\displaystyle\geq$ $\displaystyle 2-4\Bigg{(}\frac{\|\widetilde{\mathbf{P}}^{s_{k},e_{k}}(b_{k})- \widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F}}}{\min\{\|\widetilde{ \mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F}},\|\widehat{\mathbf{P}}^{s_{k}, e_{k}}(b_{k})\|_{\mathrm{F}}\}}\Bigg{)}^{2}$ $\displaystyle\geq$ $\displaystyle 2-4\frac{48^{2}C_{\mathcal{B}}^{2}(d^{2}m^{s,e}_{t}+nd+Lm^{s,e}_ {t})\log(T)}{\kappa_{k}^{2}\Delta}\geq 1,$ $$

where the first inequality follows from the reverse triangle inequality, the third inequality follows from
the definition of the event $\mathcal{B}$, ([C.1](https://arxiv.org/html/2506.21878v1#A3.Ex34)) and ([13](https://arxiv.org/html/2506.21878v1#A3.E13)),
and the final inequality follows from Assumption [2](https://arxiv.org/html/2506.21878v1#Thmassumption2) and the fact that $C_{\mathrm{SNR}}$ is a sufficiently large constant. Therefore,

$$ $\Big{\langle}\widetilde{\mathbf{P}}^{s_{k},e_{k}}(b_{k}),\widehat{\mathbf{P}}^ {s_{k},e_{k}}(b_{k})/\|\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F} }\Big{\rangle}\geq\|\widetilde{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F}} /2\geq\kappa_{k}\sqrt{\Delta/96},$ (14) $$

where the last inequality follows from ([C.1](https://arxiv.org/html/2506.21878v1#A3.Ex34)).

Step 4. Since $\{\mathbf{A}^{\prime}(t)\}_{t=1}^{T}$ is independent of $\{\mathbf{B}^{\prime}(t)\}_{t=1}^{T}$, the distribution of $\{\mathbf{A}^{\prime}(t)\}_{t=1}^{T}$ remain unaffected under the conditioning on the event $\mathcal{B}$. By the truncation in the construction of $\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})$ stated in Algorithm [2](https://arxiv.org/html/2506.21878v1#alg2),
we have that

$$ $\|\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\infty}\leq\sqrt{\frac{(e_{k}-b _{k})(b_{k}-s_{k})}{(e_{k}-s_{k})}}$ $$

Combined with ([13](https://arxiv.org/html/2506.21878v1#A3.E13)), it follows that

$$ $(e_{k}-s_{k})^{-1/2}\|\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\infty}/\| \widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F}}\leq\frac{1}{c_{3} \kappa_{k}\sqrt{\Delta}},$ $$

for some constant $c_{3}>0$. Applying Lemma [5](https://arxiv.org/html/2506.21878v1#Thmtheorem5), we obtain for any $\varepsilon>0$

$$ $\displaystyle\mathbb{P}\Big{\{}\Big{|}\frac{1}{\sqrt{e_{k}-s_{k}}}\sum_{t=s_{k }+1}^{e_{k}}\Big{\langle}\mathbf{P}(t)-\mathbf{A}^{\prime}(t),\widehat{\mathbf {P}}^{s_{k},e_{k}}(b_{k})/\|\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{ \mathrm{F}}\Big{\rangle}\Big{|}\geq\varepsilon\Big{\}}$ $\displaystyle\leq$ $\displaystyle 2\exp\Bigg{\{}-c_{4}\frac{\varepsilon^{2}}{1+\varepsilon/(c_{3} \kappa_{k}\sqrt{\Delta})}\Bigg{\}}.$ $$

where $c_{4}>0$ is a constant.
Choosing $\varepsilon=C_{3}\sqrt{\log(T)}$ for a large enough constant $C_{3}>0$, and applying Assumption [2](https://arxiv.org/html/2506.21878v1#Thmassumption2) and the fact that $C_{\mathrm{SNR}}$ is a sufficiently large constant, we finally derive that

$$ $\mathbb{P}\Big{\{}\Big{|}\frac{1}{\sqrt{e_{k}-s_{k}}}\sum_{t=s_{k}+1}^{e_{k}} \Big{\langle}\mathbf{P}(t)-\mathbf{A}^{\prime}(t),\widehat{\mathbf{P}}^{s_{k}, e_{k}}(b_{k})/\|\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F}}\Big{ \rangle}\Big{|}\geq C_{3}\sqrt{\log(T)}\Big{\}}\leq 2T^{-c_{5}},$ (15) $$

where $c_{5}>3$ is a constant.
A similar argument also demonstrates that

$$ $\mathbb{P}\Big{\{}\Big{|}\Big{\langle}\widetilde{\mathbf{P}}^{s_{k},e_{k}}(t)- \widetilde{\mathbf{A}^{\prime}}^{s_{k},e_{k}}(t),\widehat{\mathbf{P}}^{s_{k},e _{k}}(b_{k})/\|\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F}}\Big{ \rangle}\Big{|}\geq C_{3}\sqrt{\log(T)}\Bigg{\}}\leq 2T^{-c_{5}}.$ (16) $$

Step 5. We now consider the univariate time series defined for all $t\in(s_{k},e_{k})$ as

$$ $y(t)=\Big{\langle}\mathbf{A}^{\prime}(t),\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_ {k})/\|\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F}}\Big{\rangle}$ $$

and

$$ $y^{s_{k},e_{k}}(t)=\Big{\langle}\widetilde{\mathbf{A}^{\prime}}^{s_{k},e_{k}}( t),\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})/\|\widehat{\mathbf{P}}^{s_{k},e_{ k}}(b_{k})\|_{\mathrm{F}}\Big{\rangle}.$ $$

Conditional on the event $\mathcal{B}$, define the corresponding mean functions

$$ $f(t)=\mathbb{E}(y(t))=\Big{\langle}\mathbf{P}(t),\widehat{\mathbf{P}}^{s_{k},e _{k}}(b_{k})/\|\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F}}\Big{ \rangle},$ $$

and

$$ $f^{s_{k},e_{k}}(t)=\mathbb{E}(y^{s_{k},e_{k}}(t))=\Big{\langle}\widetilde{ \mathbf{P}}^{s_{k},e_{k}}(t),\widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})/\| \widehat{\mathbf{P}}^{s_{k},e_{k}}(b_{k})\|_{\mathrm{F}}\Big{\rangle}.$ $$

The function $f(t)$ is a piecewise constant on $(s_{k},e_{k}]$ with a single change point at $\eta_{k}$. Using ([14](https://arxiv.org/html/2506.21878v1#A3.E14)), we obtain that

$$ $|f^{s_{k},e_{k}}(\eta_{k})|\Big{|}\geq\kappa_{k}\sqrt{\Delta/96},$ $$

Moreover, from ([15](https://arxiv.org/html/2506.21878v1#A3.E15)), ([16](https://arxiv.org/html/2506.21878v1#A3.E16)) and an union bound argument, we have that

$$ $\displaystyle\mathbb{P}\Bigg{\{}\sup_{\begin{subarray}{c}0\leq s_{k}<t<e_{k} \leq T\\ (s_{k},e_{k})\mbox{ contains only one change point }\eta_{k}\end{subarray}}$ $\displaystyle\bigg{|}\frac{1}{\sqrt{e_{k}-s_{k}}}\sum_{t=s_{k}+1}^{e_{k}}\left (y(t)-f(t)\right)\bigg{|}$ $\displaystyle\hskip 14.22636pt\geq C_{3}\sqrt{\log(T)}\Bigg{\}}\leq 2T^{-c_{6}}$ $$

and

$$ $\mathbb{P}\Bigg{\{}\sup_{\begin{subarray}{c}0\leq s_{k}<t<e_{k}\leq T\\ (s_{k},e_{k})\mbox{ contains only one change point }\eta_{k}\end{subarray}} \big{|}y^{s_{k},e_{k}}(t)-f^{s_{k},e_{k}}(t)\big{|}\geq C_{3}\sqrt{\log(T)} \Bigg{\}}\leq 2T^{-c_{6}},$ $$

for some constant $c_{6}>0$. Applying Lemma 12 from Wang et al. (2017) with $\lambda=C\sqrt{\log(T)}$, it follows that the estimated change point $\widetilde{\eta}_{k}=\operatorname*{argmax}_{s_{k}<t<e_{k}}|y^{s_{k},e_{k}}(t)|$ is an undetected change point and satisfies for a large enough constant $C_{5}>0$,

$$ $|\widetilde{\eta}_{k}-\eta_{k}|\leq C_{5}\frac{\log(T)}{\kappa_{k}^{2}},$ $$

which completes the proof.

∎

###### Proof.

### C.2 Additional results

###### Proposition 4 .

Let $\{b_{k}\}_{k=1}^{\widetilde{K}}$ denote the output of Stage I in Algorithm [1](https://arxiv.org/html/2506.21878v1#alg1) applied to two independent adjacency tensor sequences $\{\mathbf{A}(t)\}_{t\in[T]},\{\mathbf{B}(t)\}_{t\in[T]}\subset\{0,1\}^{n\times
n
\times L}$, generated according to Definition [2](https://arxiv.org/html/2506.21878v1#Thmdefinition2) and satisfying Model [1](https://arxiv.org/html/2506.21878v1#Thmmodel1) and Assumption [2](https://arxiv.org/html/2506.21878v1#Thmassumption2). Suppose the threshold $\tau$ is chosen such that

$$ $c_{\tau,1}n\sqrt{L}\log^{3/2}(T)<\tau<c_{\tau,2}\kappa^{2}\Delta,$ (17) $$

where $c_{\tau,1},c_{\tau,2}>0$ are sufficiently large and small absolute constants, respectively.

Then, it holds that

$$ $\mathbb{P}\Big{\{}\widetilde{K}=K\mbox{ and }|b_{k}-\eta_{k}|\leq\tilde{ \epsilon},\,\forall k\in[K]\Big{\}}\geq 1-CT^{-c},$ $$

where

$$ $\tilde{\epsilon}=C_{\tilde{\epsilon}}\log(T)\bigg{\{}\frac{n\sqrt{L}\log^{1/2} (T)}{\kappa^{2}}+\frac{\sqrt{\Delta}}{\kappa}\bigg{\}},$ $$

and $C_{\tilde{\epsilon}},C,c>0$ are absolute constants.

###### Proof.

The proof presented here is a minor modification of Theorem 1 in Wang et al. (2021). For completeness, we include the full details below.

For $0\leq s<t<e\leq T$, we define the event

$$ $\mathcal{A}(s,t,e)=\Big{\{}\big{|}\big{\langle}\widetilde{\mathbf{A}}^{s,e}(t) ,\widetilde{\mathbf{B}}^{s,e}(t)\big{\rangle}-\|\widetilde{\mathbf{P}}^{s,e}(t )\|_{\mathrm{F}}^{2}\big{|}\leq C_{\mathcal{A}}\log(T)\big{(}\|\widetilde{ \mathbf{P}}^{s,e}(t)\|_{\mathrm{F}}+\log^{1/2}(T)n\sqrt{L}\big{)}\Big{\}},$ $$

where $\widetilde{\mathbf{P}}^{s,e}(t)$ is defined in ([2](https://arxiv.org/html/2506.21878v1#S2.E2)) and $C_{\mathcal{A}}>0$ is a constant.
Due to Lemma S.4 in Wang et al. (2021), it holds that $\mathbb{P}(\mathcal{A}(s,t,e)^{c})\leq C_{1}T^{-c_{1}}$ for some constants $C_{1}>0$ and $c_{1}>3$. By an union bound argument, it holds that

$$ $\mathbb{P}(\mathcal{A})=\mathbb{P}\bigg{\{}\bigcup_{0\leq s<t<e\leq T}\mathcal {A}(s,t,e)\bigg{\}}\geq 1-C_{1}T^{3-c_{1}}.$ $$

All subsequent analysis in this proof is carried out under the event $\mathcal{A}$.

We aim to demonstrate that, conditioned on the event $\mathcal{A}$ and assuming that the algorithm has accurately detected and localized change points so far, the procedure will also successfully identify any remaining undetected change point, if one exists, and estimate its location within an error of $\tilde{\epsilon}$. To this end, it suffices to consider an arbitrary time interval $0\leq s<e\leq T$ that satisfies

$$ $\eta_{r-1}\leq s<\eta_{r}<\dots<\eta_{r+q}<e\leq\eta_{r+q+1},\quad q\geq-1,$ $$

and

$$ $\max\big{\{}\min\{\eta_{r}-s,s-\eta_{r-1}\},\min\{\eta_{r+q+1}-e,e-\eta_{r+q} \}\big{\}}\leq\tilde{\epsilon},$ $$

where $q=-1$ indicates that there is no change point contained in $(s,e)$ and

$$ $\tilde{\epsilon}=C_{\tilde{\epsilon}}\log(T)\bigg{\{}\frac{n\sqrt{L}\log^{1/2} (T)}{\kappa^{2}}+\frac{\sqrt{\Delta}}{\kappa}\bigg{\}}.$ $$

with an absolute constant $C_{\tilde{\epsilon}}>0$.
By Assumption [2](https://arxiv.org/html/2506.21878v1#Thmassumption2), we have that

$$ $\tilde{\epsilon}\leq C_{\tilde{\epsilon}}\Delta\bigg{(}\frac{1}{C_{\mathrm{SNR }}^{2}\log^{1/2}(T)}+\frac{1}{C_{\mathrm{SNR}}\sqrt{nL^{1/2}}}\bigg{)}\leq \Delta/64,$ (18) $$

where the final inequality follows that $C_{\mathrm{SNR}}$ is large enough. Consequently, for any change point $\eta_{k}$, it must be that either $|\eta_{k}-s|\leq\tilde{\epsilon}$ or $|\eta_{k}-s|\geq\Delta-\tilde{\epsilon}\geq 3\Delta/4$. This implies that if
$\min\{|e-\eta_{k}|,|\eta_{k}-s|\}\leq\tilde{\epsilon}$, then $\eta_{k}$ corresponds to a change point that has already been detected and estimated within an error of at most $\tilde{\epsilon}$ during the previous induction step. We refer to a change point $\eta_{k}$ in $(s,e)$ as undetected if

$$ $\min\{\eta_{k}-s,\eta_{k}-e\}\geq 3\Delta/4.$ $$

To complete the induction step, it suffices to show that $\mathrm{SBS}\big{(}(s,e),\tau,\mathcal{J}\big{)}$ satisfies the following tow properties:
$(i)$ It does not detect any new change point within $(s,e)$ if all change points in that interval have already been detected; and
$(ii)$ It detects a point $b$ in $(s,e)$ such that $|\eta_{k}-b|\leq\tilde{\epsilon}$ if there exists at least one previously undetected change point $\eta_{k}$ in $(s,e)$.

Step 1. Assume that there are no undetected change points within the interval $(s,e)$. Then, for any interval $(\alpha^{\prime},\beta^{\prime}]\in\mathcal{J}$ with $(\alpha^{\prime},\beta^{\prime}]\subseteq(s,e]$, one of the following scenarios must hold:
$(i)$ The interval $(\alpha^{\prime},\beta^{\prime})$ contains no change points;
$(ii)$ The interval $(\alpha^{\prime},\beta^{\prime})$ contains exactly one change point $\eta_{k}$ and $\min\{\eta_{k}-\alpha^{\prime},\beta^{\prime}-\eta_{k}\}\leq\tilde{\epsilon}$;
$(iii)$ The interval $(\alpha^{\prime},\beta^{\prime})$ contains two change points $\eta_{k}$ and $\eta_{k+1}$, and $\max\{\eta_{k}-\alpha^{\prime},\beta^{\prime}-\eta_{k+1}\}\leq\tilde{\epsilon}$.

We focus on analyzing the scenario $(iii)$, as the other two scenarios are similar and more straightforward.
If scenario $(iii)$ holds, then by ([18](https://arxiv.org/html/2506.21878v1#A3.E18)), we have

$$ $\tilde{\epsilon}\leq\Delta/64\leq(\beta^{\prime}-\alpha^{\prime})/64,$ $$

This implies that the interval

$$ $(\alpha,\beta]=\big{(}\alpha^{\prime}+64^{-1}(\beta^{\prime}-\alpha^{\prime}), \beta^{\prime}-64^{-1}(\beta^{\prime}-\alpha^{\prime})\big{]},$ $$

contains no change points. Note that $\widetilde{\mathbf{P}}^{\alpha,\beta}(t)=0$ for all $t\in(\alpha,\beta)$, since there are no true change points within $(\alpha,\beta)$. Moreover, by the event $\mathcal{A}$,

$$ $\max_{\alpha<t<\beta}\big{\langle}\widetilde{\mathbf{A}}^{\alpha,\beta}(t), \widetilde{\mathbf{B}}^{\alpha,\beta}(t)\big{\rangle}\leq C_{\mathcal{A}}n \sqrt{L}\log^{3/2}(T).$ $$

Therefore, if the input parameter $\tau$ satisfies

$$ $\tau>C_{\mathcal{A}}n\sqrt{L}\log^{3/2}(T),$ $$

Algorithm [1](https://arxiv.org/html/2506.21878v1#alg1) will correctly reject the existence of undetected change points.

Step 2. Now suppose there exists a change point $\eta_{k}\in(s,e)$ such that

$$ $\min\{\eta_{k}-s,\eta_{k}-e\}\geq 3\Delta/4.$ $$

Let $a_{\mathcal{I}},b_{\mathcal{I}}$ and $\mathcal{I}^{*}$ be as defined in the procedure $\mathrm{SBS}\big{(}(s,e),\tau,\mathcal{J}\big{)}$. Denote $\mathcal{I}^{*}=(\alpha^{*\prime},\beta^{*\prime}]$. By Lemma 8 in Madrid Padilla et al. (2022), for any change point $\eta_{k}\in(s,e)$ satisfying $\min\{\eta_{k}-s,e-\eta_{k}\}\geq 3\Delta/4$,
there exists an interval $(\alpha^{\prime},\beta^{\prime}]\subseteq(s,e]$ containing only one $\eta_{k}$ such that

$$ $\eta_{k}-3\Delta/4\leq\alpha^{\prime}\leq\eta_{k}-3\Delta/16\quad\mbox{and} \quad\eta_{k}+3\Delta/16\leq\beta^{\prime}\leq\eta_{k}+3\Delta/4.$ $$

Since $(\alpha,\beta]=[\alpha^{\prime}+(\beta^{\prime}-\alpha^{\prime})/64,\beta^{
\prime}-(\beta^{\prime}-\alpha^{\prime})/64]$, we have

$$ $\eta_{k}-\Delta 3/4\leq\alpha\leq\eta_{k}-\Delta/8\quad\mbox{and}\quad\eta_{k} +\Delta/8\leq\beta\leq\eta_{k}+\Delta 3/4.$ $$

On the event $\mathcal{A}$, it holds that

$$ $\big{\langle}\widetilde{\mathbf{A}}^{\alpha,\beta}(\eta_{k}),\widetilde{ \mathbf{B}}^{\alpha,\beta}(\eta_{k})\big{\rangle}\geq\|\widetilde{\mathbf{P}}^ {\alpha,\beta}(\eta_{k})\|^{2}_{\mathrm{F}}-C_{\mathcal{A}}\log(T)\big{\{}\log ^{1/2}(T)n\sqrt{L}+\|\widetilde{\mathbf{P}}^{\alpha,\beta}(\eta_{k})\|_{ \mathrm{F}}\big{\}}.$ (19) $$

Furthermore, by Lemma [6](https://arxiv.org/html/2506.21878v1#Thmtheorem6), it hold that

$$ $\|\widetilde{\mathbf{P}}^{\alpha,\beta}(\eta_{k})\|^{2}_{\mathrm{F}}=\frac{( \eta_{k}-\alpha)(\beta-\eta_{k})}{\beta-\alpha}\kappa_{k}^{2}.$ $$

Then we can derive that

$$ $\|\widetilde{\mathbf{P}}^{\alpha,\beta}(\eta_{k})\|^{2}_{\mathrm{F}}\geq\kappa _{k}^{2}\Delta/16\quad\mbox{and}\quad\|\widetilde{\mathbf{P}}^{\alpha,\beta}( \eta_{k})\|^{2}_{\mathrm{F}}\leq 3\kappa_{k}^{2}\Delta/4.$ (20) $$

Combining ([19](https://arxiv.org/html/2506.21878v1#A3.E19)) and ([20](https://arxiv.org/html/2506.21878v1#A3.E20)), and using Assumption [2](https://arxiv.org/html/2506.21878v1#Thmassumption2) along with fact that $C_{\mathrm{SNR}}$ is a sufficiently large constant, we obtain

$$ $\displaystyle\big{\langle}\widetilde{\mathbf{A}}^{\alpha,\beta}(\eta_{k}), \widetilde{\mathbf{B}}^{\alpha,\beta}(\eta_{k})\big{\rangle}\geq\kappa_{k}^{2} \Delta/16-\kappa_{k}^{2}\Delta/64-\kappa_{k}^{2}\Delta/64\geq\kappa_{k}^{2} \Delta/32.$ $$

By the definition of $\mathcal{I}^{*}$, it follows that

$$ $a_{\mathcal{I}^{*}}=\big{\langle}\tilde{A}^{\alpha^{*},\beta^{*}}(b_{\mathcal{ I}^{*}}),\tilde{B}^{\alpha^{*},\beta^{*}}(b_{\mathcal{I}^{*}})\big{\rangle} \geq(\kappa_{\max}^{s,e})^{2}\Delta/32,$ (21) $$

where

$$ $\kappa_{\max}^{s,e}=\max\big{\{}\kappa_{k}\colon\min\{\eta_{k}-s,e-\eta_{k}\} \geq 3\Delta/4\big{\}}.$ $$

Therefore, if the threshold $\tau$ satisfies

$$ $\tau<\kappa^{2}\Delta/32,$ $$

Algorithm [1](https://arxiv.org/html/2506.21878v1#alg1) will consistently detect the existence of any previously undetected change points within the interval.

Step 3. Suppose there exists at least one undetected change point $\eta_{k}\in(s,e)$ such that

$$ $\min\{\eta_{k}-s,\eta_{k}-e\}\geq 3\Delta/4.$ $$

Let $a_{\mathcal{I}},b_{\mathcal{I}}$ and $\mathcal{I}^{*}$ be defined according to the procedure $\mathrm{SBS}\big{(}(s,e),\tau,\mathcal{J}\big{)}$, and denote $\mathcal{I}^{*}=(\alpha^{*\prime},\beta^{*\prime}]$. To complete the induction step, it suffices to establish that there exists an undetected change point $\eta_{k}\in(\alpha^{*\prime},\beta^{*\prime})$ satisfying

$$ $\min\{\eta_{k}-\alpha^{*\prime},\beta^{*\prime}-\eta_{k}\}\geq 3\Delta/4,$ (22) $$

and that

$$ $|b_{\mathcal{I}^{*}}-\eta_{k}|\leq\tilde{\epsilon}.$ (23) $$

Step 3.1. Proof of ([22](https://arxiv.org/html/2506.21878v1#A3.E22)).
Let

$$ $(\alpha^{*},\beta^{*}]=(\alpha^{*\prime}+(\beta^{*\prime}-\alpha^{*\prime})/64 ,\beta^{*\prime}-(\beta^{*\prime}-\alpha^{*\prime})/64].$ (24) $$

Assume for contradiction that

$$ $\max_{\alpha^{*}<t<\beta^{*}}\big{\|}\widetilde{\mathbf{P}}^{\alpha^{*},\beta^ {*}}(t)\big{\|}_{\mathrm{F}}^{2}<(\kappa_{\max}^{s,e})^{2}\Delta/64.$ (25) $$

Then on the event $\mathcal{A}$, we obtain

$$ $\displaystyle\max_{\alpha^{*}<t<\beta^{*}}\big{\langle}\widetilde{A}^{\alpha^{ *},\beta^{*}}(t),\widetilde{B}^{\alpha^{*},\beta^{*}}(t)\big{\rangle}$ $\displaystyle\leq$ $\displaystyle\max_{\alpha^{*}<t<\beta^{*}}\big{\|}\widetilde{\mathbf{P}}^{ \alpha^{*},\beta^{*}}(t)\big{\|}_{\mathrm{F}}^{2}+C_{\mathcal{A}}\log(T)\max_{ \alpha^{*}<t<\beta^{*}}\big{\|}\widetilde{\mathbf{P}}^{\alpha^{*},\beta^{*}}(t )\big{\|}_{\mathrm{F}}+C_{\mathcal{A}}\log^{3/2}(T)n\sqrt{L}$ $\displaystyle<$ $\displaystyle(\kappa_{\max}^{s,e})^{2}\Delta/64+C_{\mathcal{A}}\log(T)\kappa_{ \max}^{s,e}\sqrt{\Delta}/8+C_{\mathcal{A}}\log^{3/2}(T)n\sqrt{L}$ $\displaystyle\leq$ $\displaystyle(\kappa_{\max}^{s,e})^{2}\Delta/64+(\kappa_{\max}^{s,e})^{2} \Delta/128+(\kappa_{\max}^{s,e})^{2}\Delta/128=(\kappa_{\max}^{s,e})^{2}\Delta /32,$ $$

where the second inequality follows from ([25](https://arxiv.org/html/2506.21878v1#A3.E25)), the third inequality follows from Assumption [2](https://arxiv.org/html/2506.21878v1#Thmassumption2) and the fact that $C_{\mathrm{SNR}}$ is a large enough constant. This contradicts the inequality ([21](https://arxiv.org/html/2506.21878v1#A3.E21)).
Thus,

$$ $\max_{\alpha^{*}<t<\beta^{*}}\big{\|}\widetilde{\mathbf{P}}^{\alpha^{*},\beta^ {*}}(t)\big{\|}_{\mathrm{F}}^{2}\geq(\kappa_{\max}^{s,e})^{2}\Delta/64.$ (26) $$

Now, observe
$(i)$ If $(\alpha^{*},\beta^{*})$ contains at least two change points, then $\beta^{*}-\alpha^{*}\geq\Delta$.
$(ii)$ If it contains exactly one change point $\eta_{k}$, but
$\min\{\eta_{k}-\alpha^{*},\beta^{*}-\eta_{k}\}<\Delta/64$, then by Lemma [6](https://arxiv.org/html/2506.21878v1#Thmtheorem6), we would have

$$ $\displaystyle\max_{\alpha^{*}<t<\beta^{*}}\big{\|}\widetilde{P}^{\alpha^{*}, \beta^{*}}(t)\big{\|}_{\mathrm{F}}^{2}=$ $\displaystyle\big{\|}\widetilde{P}^{\alpha^{*},\beta^{*}}(\eta_{k})\big{\|}_{ \mathrm{F}}^{2}=\frac{(\eta_{k}-\alpha^{*})(\beta^{*}-\eta_{k})}{\beta^{*}- \alpha^{*}}\kappa_{k}^{2}$ $\displaystyle\leq$ $\displaystyle\min\{\beta^{*}-\eta_{k},\eta_{k}-\alpha^{*}\}\kappa_{k}^{2}<( \kappa_{\max}^{s,e})^{2}\Delta/64,$ $$

contradicting ([26](https://arxiv.org/html/2506.21878v1#A3.E26)). Therefore, it has to be the case that $\min\{\eta_{k}-\alpha^{*},\beta^{*}-\eta_{k}\}\geq\Delta/64$.
Moreover, by ([24](https://arxiv.org/html/2506.21878v1#A3.E24)), it holds that

$$ $\beta^{*,\prime}-\alpha^{*,\prime}\geq\beta^{*}-\alpha^{*}\geq\Delta/64.$ $$

Then, using Assumption [2](https://arxiv.org/html/2506.21878v1#Thmassumption2) and the the fact that $C_{\mathrm{SNR}}$ is large enough, we have that

$$ $\tilde{\epsilon}\leq C_{\tilde{\epsilon}}\Delta\bigg{(}\frac{1}{C_{\mathrm{SNR }}^{2}\log^{1/2+2\xi}(T)}+\frac{1}{C_{\mathrm{SNR}}\sqrt{nL^{1/2}}\log^{1+\xi} (T)}\bigg{)}\leq(\beta^{*,\prime}-\alpha^{*,\prime})/64,$ $$

Hence, by a similar argument as in Step 1, no previously detected change point lies in $(\alpha^{*},\beta^{*})$. Note that by ([21](https://arxiv.org/html/2506.21878v1#A3.E21)), there is at least one undetected change point in $(\alpha^{*},\beta^{*})$.

Step 3.2. Proof of ([23](https://arxiv.org/html/2506.21878v1#A3.E23)). To this end, we apply Lemma S.5 in Wang et al. (2021). Define

$$ $\lambda=\max_{\alpha^{*}<t<\beta^{*}}\Big{|}\big{\langle}\widetilde{\mathbf{A} }^{\alpha^{*},\beta^{*}}(t),\widetilde{\mathbf{B}}^{\alpha^{*},\beta^{*}}(t) \big{\rangle}-\big{\|}\widetilde{\mathbf{P}}^{\alpha^{*},\beta^{*}}(t)\big{\|} _{\mathrm{F}}^{2}\Big{|}.$ (27) $$

From ([26](https://arxiv.org/html/2506.21878v1#A3.E26)), Assumption [2](https://arxiv.org/html/2506.21878v1#Thmassumption2) and the fact that $C_{\mathrm{SNR}}$ is a sufficiently large constant, it follows that

$$ $\displaystyle c_{2}\max_{\alpha^{*}<t<\beta^{*}}\big{\|}\widetilde{\mathbf{P}} ^{\alpha^{*},\beta^{*}}(t)\big{\|}_{\mathrm{F}}^{2}/2\geq$ $\displaystyle\max\big{\{}C_{\mathcal{A}}\log(T)\max_{\alpha^{*}<t<\beta^{*}} \big{\|}\widetilde{\mathbf{P}}^{\alpha^{*},\beta^{*}}(t)\big{\|}_{\mathrm{F}},$ $\displaystyle\hskip 28.45274ptC_{\mathcal{A}}\log^{3/2}(T)n\sqrt{L}\Big{\}},$ (28) $$

where $c_{2}>0$ is a small enough constant. By the definition of the event $\mathcal{A}$, we obtain

$$ $\displaystyle\lambda\leq$ $\displaystyle C_{\mathcal{A}}\log(T)\big{\{}\max_{\alpha^{*}<t<\beta^{*}}\big{ \|}\widetilde{\mathbf{P}}^{\alpha^{*},\beta^{*}}(t)\big{\|}_{\mathrm{F}}+\log^ {1/2}(T)n\sqrt{L}\big{)}\big{\}}\leq c_{2}\max_{\alpha^{*}<t<\beta^{*}}\big{\| }\widetilde{\mathbf{P}}^{\alpha^{*},\beta^{*}}(t)\big{\|}_{\mathrm{F}}^{2},$ (29) $$

where the last inequality follows from ([C.2](https://arxiv.org/html/2506.21878v1#A3.Ex89)).
Note that ([21](https://arxiv.org/html/2506.21878v1#A3.E21)), ([27](https://arxiv.org/html/2506.21878v1#A3.E27)) and ([29](https://arxiv.org/html/2506.21878v1#A3.E29)) verify conditions (2), (3), (4) of Lemma S.5 in Wang et al. (2021), respectively. Therefore, Lemma S.5 in Wang et al. (2021) implies that there exists an undetected change point $\eta_{k}$ within $(s,e)$ such that

$$ $\displaystyle|\eta_{k}-b_{\mathcal{I}^{*}}|\leq\frac{C_{3}\Delta\lambda}{\| \widetilde{P}^{\alpha^{*},\beta^{*}}(\eta_{k})\|^{2}_{\mathrm{F}}}\quad\mbox{ and}\quad\|\widetilde{P}^{\alpha^{*},\beta^{*}}(\eta_{k})\|^{2}_{\mathrm{F}} \geq c_{4}\max_{\alpha^{*}\leq t\leq\beta^{*}}\|\widetilde{P}^{\alpha^{*}, \beta^{*}}(t)\|^{2}_{\mathrm{F}},$ (30) $$

where $C_{3},c_{4}>0$ are constants. Then combining ([29](https://arxiv.org/html/2506.21878v1#A3.E29)) and ([30](https://arxiv.org/html/2506.21878v1#A3.E30)), we can derive that

$$ $\displaystyle|\eta_{k}-b_{\mathcal{I}^{*}}|\leq$ $\displaystyle\frac{C_{3}C_{\mathcal{A}}\Delta\log(T)\big{\{}\max_{\alpha^{*}<t <\beta^{*}}\big{\|}\widetilde{\mathbf{P}}^{\alpha^{*},\beta^{*}}(t)\big{\|}_{ \mathrm{F}}+\log^{1/2}(T)n\sqrt{L}\big{)}\big{\}}\ }{c_{4}\max_{\alpha^{*}\leq t \leq\beta^{*}}\|\widetilde{\mathbf{P}}^{\alpha^{*},\beta^{*}}(t)\|^{2}_{ \mathrm{F}}}$ $\displaystyle=$ $\displaystyle\frac{C_{3}C_{\mathcal{A}}\log(T)}{c_{4}}\bigg{\{}\frac{\Delta}{ \max_{\alpha^{*}<t<\beta^{*}}\big{\|}\widetilde{\mathbf{P}}^{\alpha^{*},\beta^ {*}}(t)\big{\|}_{\mathrm{F}}}+\frac{\Delta\log^{1/2}(T)n\sqrt{L}}{\max_{\alpha ^{*}\leq t\leq\beta^{*}}\|\widetilde{P}^{\alpha^{*},\beta^{*}}(t)\|^{2}_{ \mathrm{F}}}\bigg{\}}$ $\displaystyle\leq$ $\displaystyle\frac{C_{3}C_{\mathcal{A}}\log(T)}{c_{4}}\bigg{\{}\frac{8\sqrt{ \Delta}}{\kappa}+\frac{64\log^{1/2}(T)n\sqrt{L}}{\kappa^{2}}\bigg{\}}$ $\displaystyle\leq$ $\displaystyle C_{5}\log(T)\bigg{\{}\frac{\sqrt{\Delta}}{\kappa}+\frac{\log^{1/ 2}(T)n\sqrt{L}}{\kappa^{2}}\bigg{\}},$ $$

where the second inequality follows form ([26](https://arxiv.org/html/2506.21878v1#A3.E26)) and $C_{5}>0$ is an constant. This completes the induction step and therefore, the proof.

∎

###### Lemma 5 .

Let $\{\mathbf{A}(t)\}_{t\in[T]}$ follow $\mathrm{D}\mbox{-}\mathrm{MRDPGs}$ as in Definition [2](https://arxiv.org/html/2506.21878v1#Thmdefinition2). Let $\mathbf{V}\in\mathbb{R}^{n\times n\times L}$ and $\{w_{t}\}_{t=1}^{T}\subset\mathbb{R}$ satisfy $\sum_{t=1}^{T}w_{t}^{2}=1$. Then for any $\varepsilon>0$, it holds that

$$ $\mathbb{P}\bigg{(}\bigg{|}\sum_{t=1}^{T}w_{t}\big{\langle}\mathbf{A}(t)- \mathbf{P}(t),\mathbf{V}/\|\mathbf{V}\|_{\mathrm{F}}\big{\rangle}\bigg{|}\geq \varepsilon\bigg{)}\leq 2\exp\bigg{(}-c\frac{\varepsilon^{2}}{1+\varepsilon\| \mathbf{V}\|_{\mathrm{F}}^{-1}\|\mathbf{V}\|_{\infty}\max_{1\leq t\leq T}|w_{t }|}\bigg{)},$ $$

where $c>0$ is an absolute constant and $\|\mathbf{V}\|_{\infty}=\max_{i,j\in[n],l\in[L]}|V_{i,j,l}|$.

###### Proof.

By definition of the tensor inner product, we have that

$$ $\sum_{t=1}^{T}w_{t}\big{\langle}\mathbf{A}(t)-\mathbf{P}(t),\mathbf{V}/\| \mathbf{V}\|_{\mathrm{F}}\big{\rangle}=\sum_{t=1}^{T}\sum_{i=1}^{n}\sum_{j=1}^ {n}\sum_{l=1}^{L}\|\mathbf{V}\|_{\mathrm{F}}^{-1}w_{t}\mathbf{V}_{i,j,l}\Big{ \{}\big{(}\mathbf{A}(t)\big{)}_{i,j,l}-\big{(}\mathbf{P}(t)\big{)}_{i,j,l}\Big {\}}.$ $$

We can derive that

$$ $\sum_{t=1}^{T}\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{l=1}^{L}\|\mathbf{V}\|_{ \mathrm{F}}^{-2}w_{t}^{2}\mathbf{V}_{i,j,l}^{2}=1,$ $$

and

$$ $\max_{t\in[T],i,j\in[n],l\in[L]}\Big{\{}\|\mathbf{V}\|_{\mathrm{F}}^{-1}w_{t} \mathbf{V}_{i,j,l}\Big{\}}\leq\|\mathbf{V}\|_{\mathrm{F}}^{-1}\|\mathbf{V}\|_{ \infty}\max_{1\leq t\leq T}|w_{t}|.$ $$

Since

$$ $\{\big{(}\mathbf{A}(t)\big{)}_{i,j,l}-\big{(}\mathbf{P}(t)\big{)}_{i,j,l}\}_{i ,j\in[n],l\in[L],t\in[T]}$ $$

are mutually independent centered Bernoulli random variables By Bernstein inequality (e.g. Theorem 2.8.2 in Vershynin, 2018), it holds with an absolute constant $c_{0}>0$ that

$$ $\mathbb{P}\bigg{(}\bigg{|}\sum_{t=1}^{T}w_{t}\big{\langle}\mathbf{A}(t)- \mathbf{P}(t),\mathbf{V}/\|\mathbf{V}\|_{\mathrm{F}}\big{\rangle}\bigg{|}\geq \varepsilon\bigg{)}\leq 2\exp\bigg{(}-c_{0}\frac{\varepsilon^{2}}{1+ \varepsilon\|\mathbf{V}\|_{\mathrm{F}}^{-1}\|\mathbf{V}\|_{\infty}\max_{1\leq t \leq T}|w_{t}|}\bigg{)},$ $$

which completes the proof.

∎

###### Lemma 6 .

Suppose the adjacency tensor sequence $\{\mathbf{B}(t)\}_{t\in[T]}\subset\{0,1\}^{n\times n\times L}$ is generated according to Definition [2](https://arxiv.org/html/2506.21878v1#Thmdefinition2) and satisfy Model [1](https://arxiv.org/html/2506.21878v1#Thmmodel1). For any $0\leq s<t<e\leq T$, let $\widetilde{\mathbf{P}}^{s,e}(t)$ be defined
as in ([2](https://arxiv.org/html/2506.21878v1#S2.E2)).
If $(s,e)$ contains exactly one change point $\eta_{k}$, then for any $t\in(s,e)$

$$ $\big{\|}\widetilde{\mathbf{P}}^{s,e}(t)\big{\|}_{\mathrm{F}}^{2}=\begin{cases} \frac{t-s}{(e-s)(e-t)}(e-\eta_{k})^{2}\kappa_{k}^{2},&s<t\leq\eta_{k},\\ \frac{e-t}{(e-s)(t-s)}(\eta_{k}-s)^{2}\kappa_{k}^{2},&\eta_{k}<t<e.\end{cases}$ $$

###### Proof.

This follows directly the definition of $\widetilde{\mathbf{P}}^{s,e}(t)$ in ([2](https://arxiv.org/html/2506.21878v1#S2.E2)).
∎

###### Proposition 4 .

###### Proof.

###### Lemma 5 .

###### Proof.

###### Lemma 6 .

###### Proof.

## Appendix D Proofs of Theorems 2 and 3

###### Proof.

Step 1. Preliminary bounds.
We first define the event

$$ $\mathcal{A}=\Big{\{}\widetilde{K}=K\mbox{ and }|\widetilde{\eta}_{k}-\eta_{k}| \leq\epsilon_{k},\,\forall k\in[K]\Big{\}},\quad\mbox{where}\quad\epsilon_{k}= C_{\epsilon}\frac{\log(T)}{\kappa_{k}^{2}}.$ $$

Then by Theorem [1](https://arxiv.org/html/2506.21878v1#Thmtheorem1), it holds that

$$ $\mathbb{P}\{\mathcal{A}\}\geq 1-C_{0}T^{-c_{0}},$ $$

and $C_{\tilde{\epsilon}},C_{0},c_{0}>0$ are absolute constants. Since $\mathcal{A}$ holds with probability tending to $1$ as $T\to\infty$,we
condition the remainder of the proof on $\mathcal{A}$.

From $\mathcal{A}$, Assumption [2](https://arxiv.org/html/2506.21878v1#Thmassumption2) and the fact that $C_{\mathrm{SNR}}$ is a sufficiently large constant, we have for all $k\in[K]$ that
$\eta_{k}\in[\widetilde{\eta}_{k-1},\widetilde{\eta}_{k+1}]$,

$$ $\eta_{k}-\widetilde{\eta}_{k-1}\geq\eta_{k}-\eta_{k-1}-|\eta_{k-1}-\widetilde{ \eta}_{k-1}|\geq\Delta-C_{\epsilon}\frac{\log(T)}{\kappa^{2}}\geq\Delta-\Delta /6=5\Delta/6,$ (31) $$

and

$$ $\widetilde{\eta}_{k+1}-\eta_{k}\geq\eta_{k+1}-\eta_{k}-|\eta_{k+1}-\widetilde{ \eta}_{k+1}|\geq\Delta-C_{\epsilon}\frac{\log(T)}{\kappa^{2}}\geq\Delta-\Delta /6=5\Delta/6.$ (32) $$

Similarly, we can derive that for any $k\in[K+1]$

$$ $\widetilde{\eta}_{k}-\widetilde{\eta}_{k-1}\geq 2\Delta/3.$ (33) $$

As a result, the working interval

$$ $(\tilde{s}_{k},\tilde{e}_{k}]=(\widetilde{\eta}_{k-1}+(\widetilde{\eta}_{k}- \widetilde{\eta}_{k-1})/2,\;\widetilde{\eta}_{k+1}-(\widetilde{\eta}_{k+1}- \widetilde{\eta}_{k})/2],$ $$

contains exactly one change point $\eta_{k}$.

Next, by Theorem 4 and Lemma 5 in Wang et al. (2023), and Assumption [1](https://arxiv.org/html/2506.21878v1#Thmassumption1) $(i)$ and $(iii)$, for any $k\in[K+1]$, we have

$$ $\big{\|}\widehat{\mathbf{P}}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}- \mathbf{P}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}\big{\|}_{\mathrm{F}}= O_{p}\bigg{(}\sqrt{\frac{(d^{2}m_{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}} +nd+Lm_{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}})\log(T)}{\widetilde{\eta} _{k}-\widetilde{\eta}_{k-1}}}\bigg{)},$ $$

for some absolute constant $C_{1}>0$. For any $k\in[K+1]$, by ([31](https://arxiv.org/html/2506.21878v1#A4.E31)) and ([32](https://arxiv.org/html/2506.21878v1#A4.E32)), each interval $(\widetilde{\eta}_{k-1},\widetilde{\eta}_{k})$
contains at most two true change points $\eta_{k}-1$ and $\eta_{k}$. Consequently, we have that

$$ $m_{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}\leq m_{k-1}+m_{k}+m_{k+1}\leq 3 m_{\max}.$ $$

Thus, it holds that

$$ $\big{\|}\widehat{\mathbf{P}}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}- \mathbf{P}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}\big{\|}_{\mathrm{F}} \leq C_{2}\sqrt{\frac{(d^{2}m_{\max}+nd+Lm_{\max})\log(T)}{\widetilde{\eta}_{k }-\widetilde{\eta}_{k-1}}},\quad\forall k\in[K+1],$ (34) $$

with an absolute constant $C_{2}>0$.

Step 2. Characterization of bias.
From ([31](https://arxiv.org/html/2506.21878v1#A4.E31)) and ([32](https://arxiv.org/html/2506.21878v1#A4.E32)), for any $k\in[K]$, the interval $(\widetilde{\eta}_{k-1},\widetilde{\eta}_{k+1})$ may contain one, two or three change points.
One example that contains three change points is illustrated in the figure below. We analyze the biases in this case. The analyses for the other scenarios are similar but simpler and therefore omitted.

In the following, we analyze three types of bias terms. Denote $\alpha_{T}=\log(T)$, then $\alpha_{T}\to\infty$ as $T\to\infty$. Observe that

$$ $\displaystyle\big{\|}\mathbf{P}(\eta_{k})-\mathbf{P}^{\widetilde{\eta}_{k-1}, \widetilde{\eta}_{k}}\big{\|}_{\mathrm{F}}$ $\displaystyle=$ $\displaystyle\bigg{\|}\mathbf{P}(\eta_{k})-\frac{\eta_{k-1}-\widetilde{\eta}_{ k-1}}{\widetilde{\eta}_{k}-\widetilde{\eta}_{k-1}}\mathbf{P}(\eta_{k-1})-\frac {\eta_{k}-\eta_{k-1}}{\widetilde{\eta}_{k}-\widetilde{\eta}_{k-1}}\mathbf{P}( \eta_{k})-\frac{\widetilde{\eta}_{k}-\eta_{k}}{\widetilde{\eta}_{k}-\widetilde {\eta}_{k-1}}\mathbf{P}(\eta_{k+1})\bigg{\|}_{\mathrm{F}}$ $\displaystyle\leq$ $\displaystyle\frac{\eta_{k-1}-\widetilde{\eta}_{k-1}}{\widetilde{\eta}_{k}- \widetilde{\eta}_{k-1}}\big{\|}\mathbf{P}(\eta_{k})-\mathbf{P}(\eta_{k-1})\big {\|}_{\mathrm{F}}+\frac{\widetilde{\eta}_{k}-\eta_{k}}{\widetilde{\eta}_{k}- \widetilde{\eta}_{k-1}}\big{\|}\mathbf{P}(\eta_{k+1})-\mathbf{P}(\eta_{k})\big {\|}_{\mathrm{F}}$ $\displaystyle\leq$ $\displaystyle\frac{3C_{\epsilon}\log(T)}{2\Delta\kappa_{k-1}^{2}}\kappa_{k-1}+ \frac{3C_{\epsilon}\log(T)}{2\Delta\kappa_{k}^{2}}\kappa_{k}\leq\alpha_{T}^{-1 }\kappa_{k},$ (35) $$

where the second inequity follows from the event $\mathcal{A}$ and ([33](https://arxiv.org/html/2506.21878v1#A4.E33)), and the last inequality follows from Assumption [2](https://arxiv.org/html/2506.21878v1#Thmassumption2) and the fact that $C_{\mathrm{SNR}}$ is a large enough constant.
Similarly, we have that

$$ $\displaystyle\big{\|}\mathbf{P}(\eta_{k+1})-\mathbf{P}^{\widetilde{\eta}_{k}, \widetilde{\eta}_{k+1}}\big{\|}_{\mathrm{F}}=$ $\displaystyle\bigg{\|}\mathbf{P}(\eta_{k+1})-\frac{\eta_{k+1}-\widetilde{\eta} _{k}}{\widetilde{\eta}_{k+1}-\widetilde{\eta}_{k}}\mathbf{P}(\eta_{k+1})-\frac {\widetilde{\eta}_{k+1}-\eta_{k+1}}{\widetilde{\eta}_{k+1}-\widetilde{\eta}_{k }}\mathbf{P}(\eta_{k+2})\bigg{\|}_{\mathrm{F}}$ $\displaystyle=$ $\displaystyle\frac{\widetilde{\eta}_{k+1}-\eta_{k+1}}{\widetilde{\eta}_{k+1}- \widetilde{\eta}_{k}}\big{\|}\mathbf{P}(\eta_{k+1})-\mathbf{P}(\eta_{k+2})\big {\|}_{\mathrm{F}}$ $\displaystyle\leq$ $\displaystyle\frac{3C_{\epsilon}\log(T)}{2\Delta\kappa_{k+1}^{2}}\kappa_{k+1} \leq\alpha_{T}^{-1}\kappa_{k},$ (36) $$

and

$$ $\displaystyle\big{\|}\mathbf{P}(\eta_{k+1})-\mathbf{P}^{\widetilde{\eta}_{k-1} ,\widetilde{\eta}_{k}}\big{\|}_{\mathrm{F}}$ $\displaystyle=$ $\displaystyle\bigg{\|}\mathbf{P}(\eta_{k+1})-\frac{\eta_{k-1}-\widetilde{\eta} _{k-1}}{\widetilde{\eta}_{k}-\widetilde{\eta}_{k-1}}\mathbf{P}(\eta_{k-1})- \frac{\eta_{k}-\eta_{k-1}}{\widetilde{\eta}_{k}-\widetilde{\eta}_{k-1}}\mathbf {P}(\eta_{k})-\frac{\widetilde{\eta}_{k}-\eta_{k}}{\widetilde{\eta}_{k}- \widetilde{\eta}_{k-1}}\mathbf{P}(\eta_{k+1})\bigg{\|}_{\mathrm{F}}$ $\displaystyle\leq$ $\displaystyle\frac{\eta_{k-1}-\widetilde{\eta}_{k-1}}{\widetilde{\eta}_{k}- \widetilde{\eta}_{k-1}}\big{\|}\mathbf{P}(\eta_{k+1})-\mathbf{P}(\eta_{k-1}) \big{\|}_{\mathrm{F}}+\frac{\eta_{k}-\eta_{k-1}}{\widetilde{\eta}_{k}- \widetilde{\eta}_{k-1}}\big{\|}\mathbf{P}(\eta_{k+1})-\mathbf{P}(\eta_{k})\big {\|}_{\mathrm{F}}$ $\displaystyle\leq$ $\displaystyle\frac{3C_{\epsilon}\log(T)}{2\Delta\kappa_{k-1}^{2}}(\kappa_{k-1} +\kappa_{k})+\kappa_{k}\leq\alpha_{T}^{-1}\kappa_{k}+\kappa_{k}\leq C_{3} \kappa_{k},$ (37) $$

for some constant $C_{3}>0$.

Step 3. Uniform tightness of $\kappa_{k}^{2}|\widehat{\eta}_{k}-\eta_{k}|$. In this step, we show that $\kappa_{k}^{2}|\widehat{\eta}_{k}-\eta_{k}|=O_{p}(1)$. Let $r=\widehat{\eta}_{k}-\eta_{k}$ and without loss of generality, assume $r\geq 0$. Our goal is to establish that

$$ $r\kappa_{k}^{2}=O_{p}(1)$ $$

If $r\kappa_{k}^{2}<1$, the conclusion holds trivially. Thus, for the remainder of the argument, we assume that $r\kappa_{k}^{2}\geq 1$.
Since $\widehat{\eta}_{k}=\eta_{k}+r$, it follows that

$$ $\mathcal{Q}_{k}(\eta_{k}+r)-\mathcal{Q}_{k}(\eta_{k})\leq 0.$ $$

Now observe that

$$ $\displaystyle\mathcal{Q}_{k}(\eta_{k}+r)-\mathcal{Q}_{k}(\eta_{k})=$ $\displaystyle\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\big{\|}\mathbf{A}(t)-\widehat{ \mathbf{P}}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}\big{\|}_{\mathrm{F}} ^{2}-\big{\|}\mathbf{A}(t)-\widehat{\mathbf{P}}^{\widetilde{\eta}_{k}, \widetilde{\eta}_{k+1}}\big{\|}_{\mathrm{F}}^{2}$ $\displaystyle=$ $\displaystyle\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\Big{\{}\big{\|}\mathbf{A}(t)- \widehat{\mathbf{P}}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}\big{\|}_{ \mathrm{F}}^{2}-\big{\|}\mathbf{A}(t)-\mathbf{P}^{\widetilde{\eta}_{k-1}, \widetilde{\eta}_{k}}\big{\|}_{\mathrm{F}}^{2}\Big{\}}$ $\displaystyle\hskip 14.22636pt-\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\Big{\{}\big{\| }\mathbf{A}(t)-\widehat{\mathbf{P}}^{\widetilde{\eta}_{k},\widetilde{\eta}_{k+ 1}}\big{\|}_{\mathrm{F}}^{2}-\big{\|}\mathbf{A}(t)-\mathbf{P}^{\widetilde{\eta }_{k},\widetilde{\eta}_{k+1}}\big{\|}_{\mathrm{F}}^{2}\Big{\}}$ $\displaystyle\hskip 14.22636pt+\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\Big{\{}\big{\| }\mathbf{A}(t)-\mathbf{P}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}\big{\| }_{\mathrm{F}}^{2}-\big{\|}\mathbf{A}(t)-\mathbf{P}(\eta_{k})\big{\|}_{\mathrm {F}}^{2}\Big{\}}$ $\displaystyle\hskip 14.22636pt-\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\Big{\{}\big{\| }\mathbf{A}(t)-\mathbf{P}^{\widetilde{\eta}_{k},\widetilde{\eta}_{k+1}}\big{\| }_{\mathrm{F}}^{2}-\big{\|}\mathbf{A}(t)-\mathbf{P}(\eta_{k+1})\big{\|}_{ \mathrm{F}}^{2}\Big{\}}$ $\displaystyle\hskip 14.22636pt+\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\Big{\{}\big{\| }\mathbf{A}(t)-\mathbf{P}(\eta_{k})\big{\|}_{\mathrm{F}}^{2}-\big{\|}\mathbf{A }(t)-\mathbf{P}(\eta_{k+1})\big{\|}_{\mathrm{F}}^{2}\Big{\}}$ $\displaystyle=$ $\displaystyle I-II+III-IV+V.$ $$

Therefore, we have that

$$ $V\leq-I+II-III+IV\leq|I|+|II|+|III|+|IV|.$ (38) $$

Step 3.1. Order of magnitude of $I$. We start by analyzing the term

$$ $\displaystyle I=$ $\displaystyle\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\Big{\{}\big{\|}\mathbf{A}(t)- \widehat{\mathbf{P}}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}\big{\|}_{ \mathrm{F}}^{2}-\big{\|}\mathbf{A}(t)-\mathbf{P}^{\widetilde{\eta}_{k-1}, \widetilde{\eta}_{k}}\big{\|}_{\mathrm{F}}^{2}\Big{\}}$ $\displaystyle=$ $\displaystyle\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\big{\|}\widehat{\mathbf{P}}^{ \widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}-\mathbf{P}^{\widetilde{\eta}_{k-1 },\widetilde{\eta}_{k}}\big{\|}_{\mathrm{F}}^{2}-2\sum_{t=\eta_{k}+1}^{\eta_{k }+r}\big{\langle}\mathbf{A}(t)-\mathbf{P}^{\widetilde{\eta}_{k-1},\widetilde{ \eta}_{k}},\widehat{\mathbf{P}}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}- \mathbf{P}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}\big{\rangle}$ $\displaystyle=$ $\displaystyle\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\big{\|}\widehat{\mathbf{P}}^{ \widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}-\mathbf{P}^{\widetilde{\eta}_{k-1 },\widetilde{\eta}_{k}}\big{\|}_{\mathrm{F}}^{2}-2\sum_{t=\eta_{k}+1}^{\eta_{k }+r}\big{\langle}\mathbf{A}(t)-\mathbf{P}(\eta_{k+1}),\widehat{\mathbf{P}}^{ \widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}-\mathbf{P}^{\widetilde{\eta}_{k-1 },\widetilde{\eta}_{k}}\big{\rangle}$ $\displaystyle\hskip 14.22636pt-2\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\big{\langle} \mathbf{P}(\eta_{k+1})-\mathbf{P}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k} },\widehat{\mathbf{P}}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}-\mathbf{P }^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}\big{\rangle}$ $\displaystyle=$ $\displaystyle I.1-2I.2-2I.3.$ (39) $$

By ([34](https://arxiv.org/html/2506.21878v1#A4.E34)), we have that

$$ $\displaystyle\big{\|}\widehat{\mathbf{P}}^{\widetilde{\eta}_{k-1},\widetilde{ \eta}_{k}}-\mathbf{P}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}\big{\|}_{ \mathrm{F}}^{2}\leq$ $\displaystyle C_{2}^{2}\frac{(d^{2}m_{\max}+nd+Lm_{\max})\log(T)}{\widetilde{ \eta}_{k}-\widetilde{\eta}_{k-1}}$ $\displaystyle\leq$ $\displaystyle 3C_{2}^{2}\frac{(d^{2}m_{\max}+nd+Lm_{\max})\log(T)}{2\Delta} \leq\alpha_{T}^{-1}\kappa_{k}^{2},$ (40) $$

where the second inequality is by ([33](https://arxiv.org/html/2506.21878v1#A4.E33)) and the last inequality follows from Assumption [2](https://arxiv.org/html/2506.21878v1#Thmassumption2) and the fact that $C_{\mathrm{SNR}}$ is a sufficiently large constant.
This yields that

$$ $\displaystyle|I.1|=$ $\displaystyle\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\big{\|}\widehat{\mathbf{P}}^{ \widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}-\mathbf{P}^{\widetilde{\eta}_{k-1 },\widetilde{\eta}_{k}}\big{\|}_{\mathrm{F}}^{2}=O_{p}(r\alpha_{T}^{-1}\kappa_ {k}^{2}),$ (41) $$

We now turn to the term $I.2$ in ([D](https://arxiv.org/html/2506.21878v1#A4.Ex123)). By Lemma [5](https://arxiv.org/html/2506.21878v1#Thmtheorem5) and ([D](https://arxiv.org/html/2506.21878v1#A4.Ex127)), we obtain that

$$ $\displaystyle|I.2|=O_{p}\Big{(}r^{1/2}\big{\|}\widehat{\mathbf{P}}^{\widetilde {\eta}_{k-1},\widetilde{\eta}_{k}}-\mathbf{P}^{\widetilde{\eta}_{k-1}, \widetilde{\eta}_{k}}\big{\|}_{\mathrm{F}}\Big{)}=O_{p}(r^{1/2}\alpha_{T}^{-1/ 2}\kappa_{k}).$ (42) $$

Next, by the Cauchy–Schwarz inequality, we derive that

$$ $\displaystyle|I.3|\leq r\big{\|}\mathbf{P}(\eta_{k+1})-\mathbf{P}^{\widetilde{ \eta}_{k-1},\widetilde{\eta}_{k}}\big{\|}_{\mathrm{F}}\big{\|}\widehat{\mathbf {P}}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}-\mathbf{P}^{\widetilde{\eta }_{k-1},\widetilde{\eta}_{k}}\big{\|}_{\mathrm{F}}=O_{p}(r\alpha_{T}^{-1/2} \kappa_{k}^{2}),$ (43) $$

where the last inequality follows from ([D](https://arxiv.org/html/2506.21878v1#A4.Ex111)) and ([D](https://arxiv.org/html/2506.21878v1#A4.Ex127)).

Combining ([D](https://arxiv.org/html/2506.21878v1#A4.Ex123)), ([41](https://arxiv.org/html/2506.21878v1#A4.E41)), ([42](https://arxiv.org/html/2506.21878v1#A4.E42)) and ([43](https://arxiv.org/html/2506.21878v1#A4.E43)), we conclude that

$$ $\displaystyle(I)=o_{p}(r\kappa_{k}^{2}+r^{1/2}\kappa_{k}).$ (44) $$

Step 3.2. Order of magnitude of $II$.
We now analyze the term

$$ $\displaystyle II=$ $\displaystyle\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\Big{\{}\big{\|}\mathbf{A}(t)- \widehat{\mathbf{P}}^{\widetilde{\eta}_{k},\widetilde{\eta}_{k+1}}\big{\|}_{ \mathrm{F}}^{2}-\big{\|}\mathbf{A}(t)-\mathbf{P}^{\widetilde{\eta}_{k}, \widetilde{\eta}_{k+1}}\big{\|}_{\mathrm{F}}^{2}\Big{\}}$ $\displaystyle=$ $\displaystyle\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\big{\|}\widehat{\mathbf{P}}^{ \widetilde{\eta}_{k},\widetilde{\eta}_{k+1}}-\mathbf{P}^{\widetilde{\eta}_{k}, \widetilde{\eta}_{k+1}}\big{\|}_{\mathrm{F}}^{2}-2\sum_{t=\eta_{k}+1}^{\eta_{k }+r}\big{\langle}\mathbf{A}(t)-\mathbf{P}^{\widetilde{\eta}_{k},\widetilde{ \eta}_{k+1}},\widehat{\mathbf{P}}^{\widetilde{\eta}_{k},\widetilde{\eta}_{k+1} }-\mathbf{P}^{\widetilde{\eta}_{k},\widetilde{\eta}_{k+1}}\big{\rangle}$ $\displaystyle=$ $\displaystyle\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\big{\|}\widehat{\mathbf{P}}^{ \widetilde{\eta}_{k},\widetilde{\eta}_{k+1}}-\mathbf{P}^{\widetilde{\eta}_{k}, \widetilde{\eta}_{k+1}}\big{\|}_{\mathrm{F}}^{2}-2\sum_{t=\eta_{k}+1}^{\eta_{k }+r}\big{\langle}\mathbf{A}(t)-\mathbf{P}(\eta_{k+1}),\widehat{\mathbf{P}}^{ \widetilde{\eta}_{k},\widetilde{\eta}_{k+1}}-\mathbf{P}^{\widetilde{\eta}_{k}, \widetilde{\eta}_{k+1}}\big{\rangle}$ $\displaystyle\hskip 14.22636pt-2\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\big{\langle} \mathbf{P}(\eta_{k+1})-\mathbf{P}^{\widetilde{\eta}_{k},\widetilde{\eta}_{k+1} },\widehat{\mathbf{P}}^{\widetilde{\eta}_{k},\widetilde{\eta}_{k+1}}-\mathbf{P }^{\widetilde{\eta}_{k},\widetilde{\eta}_{k+1}}\big{\rangle}$ $\displaystyle=$ $\displaystyle II.1-2II.2-2II.3.$ (45) $$

By ([34](https://arxiv.org/html/2506.21878v1#A4.E34)), we have that

$$ $\displaystyle\big{\|}\widehat{\mathbf{P}}^{\widetilde{\eta}_{k},\widetilde{ \eta}_{k+1}}-\mathbf{P}^{\widetilde{\eta}_{k},\widetilde{\eta}_{k+1}}\big{\|}_ {\mathrm{F}}^{2}\leq$ $\displaystyle C_{2}^{2}\frac{(d^{2}m_{\max}+nd+Lm_{\max})\log(T)}{\widetilde{ \eta}_{k+1}-\widetilde{\eta}_{k}}$ $\displaystyle\leq$ $\displaystyle 3C_{2}^{2}\frac{(d^{2}m_{\max}+nd+Lm_{\max})\log(T)}{2\Delta} \leq\alpha_{T}^{-1}\kappa_{k}^{2},$ (46) $$

where the second inequality follows from ([33](https://arxiv.org/html/2506.21878v1#A4.E33)) and the last inequality follows from Assumption [2](https://arxiv.org/html/2506.21878v1#Thmassumption2) and the fact that $C_{\mathrm{SNR}}$ is a sufficiently large constant.
It then follows that

$$ $\displaystyle|II.1|=$ $\displaystyle\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\big{\|}\widehat{\mathbf{P}}^{ \widetilde{\eta}_{k},\widetilde{\eta}_{k+1}}-\mathbf{P}^{\widetilde{\eta}_{k}, \widetilde{\eta}_{k+1}}\big{\|}_{\mathrm{F}}^{2}=O_{p}(r\alpha_{T}^{-1}\kappa_ {k}^{2}),$ (47) $$

To control $II.2$, by Lemma [5](https://arxiv.org/html/2506.21878v1#Thmtheorem5) and ([D](https://arxiv.org/html/2506.21878v1#A4.Ex132)), we obtain that

$$ $\displaystyle|II.2|=O_{p}\Big{(}r^{1/2}\big{\|}\widehat{\mathbf{P}}^{ \widetilde{\eta}_{k},\widetilde{\eta}_{k+1}}-\mathbf{P}^{\widetilde{\eta}_{k}, \widetilde{\eta}_{k+1}}\big{\|}_{\mathrm{F}}\Big{)}=O_{p}(r^{1/2}\alpha_{T}^{- 1/2}\kappa_{k}).$ (48) $$

Next, by the Cauchy–Schwarz inequality, we derive that

$$ $\displaystyle|II.3|\leq r\big{\|}\mathbf{P}(\eta_{k+1})-\mathbf{P}^{\widetilde {\eta}_{k},\widetilde{\eta}_{k+1}}\big{\|}_{\mathrm{F}}\big{\|}\widehat{ \mathbf{P}}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}-\mathbf{P}^{ \widetilde{\eta}_{k},\widetilde{\eta}_{k+1}}\big{\|}_{\mathrm{F}}=O_{p}(r \alpha_{T}^{-3/2}\kappa_{k}^{2}),$ (49) $$

where the last inequality follows from ([D](https://arxiv.org/html/2506.21878v1#A4.Ex109)) and ([D](https://arxiv.org/html/2506.21878v1#A4.Ex132)).

Combining ([D](https://arxiv.org/html/2506.21878v1#A4.Ex128)), ([47](https://arxiv.org/html/2506.21878v1#A4.E47)), ([48](https://arxiv.org/html/2506.21878v1#A4.E48)) and ([49](https://arxiv.org/html/2506.21878v1#A4.E49)), we conclude that

$$ $\displaystyle|II|=o_{p}(r\kappa_{k}^{2}+r^{1/2}\kappa_{k}).$ (50) $$

Step 3.3. Order of magnitude of $III$. We now analyze the term

$$ $\displaystyle III=$ $\displaystyle\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\Big{\{}\big{\|}\mathbf{A}(t)- \mathbf{P}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}\big{\|}_{\mathrm{F}}^ {2}-\big{\|}\mathbf{A}(t)-\mathbf{P}(\eta_{k})\big{\|}_{\mathrm{F}}^{2}\Big{\}}$ $\displaystyle=$ $\displaystyle\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\big{\|}\mathbf{P}(\eta_{k})- \mathbf{P}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}\big{\|}_{\mathrm{F}}^ {2}-2\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\big{\langle}\mathbf{A}(t)-\mathbf{P}( \eta_{k}),\mathbf{P}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}-\mathbf{P}( \eta_{k})\big{\rangle}$ $\displaystyle=$ $\displaystyle\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\big{\|}\mathbf{P}(\eta_{k})- \mathbf{P}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}\big{\|}_{\mathrm{F}}^ {2}-2\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\big{\langle}\mathbf{P}(\eta_{k+1})- \mathbf{P}(\eta_{k}),\mathbf{P}^{\widetilde{\eta}_{k-1},\widetilde{\eta}_{k}}- \mathbf{P}(\eta_{k})\big{\rangle}$ $\displaystyle\hskip 14.22636pt-2\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\big{\langle} \mathbf{A}(t)-\mathbf{P}(\eta_{k+1}),\mathbf{P}^{\widetilde{\eta}_{k-1}, \widetilde{\eta}_{k}}-\mathbf{P}(\eta_{k})\big{\rangle}$ $\displaystyle=$ $\displaystyle III.1-2III.2-2III.3.$ (51) $$

From ([D](https://arxiv.org/html/2506.21878v1#A4.Ex106)), we obtain that

$$ $\displaystyle|III.1|=O_{p}(r\alpha_{T}^{-2}\kappa_{k}^{2}).$ (52) $$

Using the Cauchy–Schwarz inequality and again ([D](https://arxiv.org/html/2506.21878v1#A4.Ex106)), we have that

$$ $\displaystyle|III.2|\leq r\big{\|}\mathbf{P}(\eta_{k+1})-\mathbf{P}(\eta_{k}) \big{\|}_{\mathrm{F}}\big{\|}\mathbf{P}^{\widetilde{\eta}_{k-1},\widetilde{ \eta}_{k}}-\mathbf{P}(\eta_{k})\big{\|}_{\mathrm{F}}=O_{p}(r\alpha_{T}^{-1} \kappa_{k}^{2}).$ (53) $$

To bound $III.3$, by Lemma [5](https://arxiv.org/html/2506.21878v1#Thmtheorem5) and ([D](https://arxiv.org/html/2506.21878v1#A4.Ex106)), we get that

$$ $\displaystyle|III.3|=O_{p}\Big{(}r^{1/2}\big{\|}\mathbf{P}^{\widetilde{\eta}_{ k-1},\widetilde{\eta}_{k}}-\mathbf{P}(\eta_{k})\big{\|}_{\mathrm{F}}\Big{)}=O_ {p}(r^{1/2}\alpha_{T}^{-1/2}\kappa_{k}).$ (54) $$

Combining ([D](https://arxiv.org/html/2506.21878v1#A4.Ex133)), ([52](https://arxiv.org/html/2506.21878v1#A4.E52)), ([53](https://arxiv.org/html/2506.21878v1#A4.E53)) and ([54](https://arxiv.org/html/2506.21878v1#A4.E54)), we conclude that

$$ $\displaystyle|III|=o_{p}(r\kappa_{k}^{2}+r^{1/2}\kappa_{k}).$ (55) $$

Step 3.4. Order of magnitude of $IV$. Consider the term

$$ $\displaystyle IV=$ $\displaystyle\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\Big{\{}\big{\|}\mathbf{A}(t)- \mathbf{P}^{\widetilde{\eta}_{k},\widetilde{\eta}_{k+1}}\big{\|}_{\mathrm{F}}^ {2}-\big{\|}\mathbf{A}(t)-\mathbf{P}(\eta_{k+1})\big{\|}_{\mathrm{F}}^{2}\Big{\}}$ $\displaystyle=$ $\displaystyle\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\big{\|}\mathbf{P}(\eta_{k+1})- \mathbf{P}^{\widetilde{\eta}_{k},\widetilde{\eta}_{k+1}}\big{\|}_{\mathrm{F}}^ {2}-2\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\big{\langle}\mathbf{A}(t)-\mathbf{P}( \eta_{k+1}),\mathbf{P}^{\widetilde{\eta}_{k},\widetilde{\eta}_{k+1}}-\mathbf{P }(\eta_{k+1})\big{\rangle}$ $\displaystyle=$ $\displaystyle IV.1-2IV.2.$ (56) $$

By ([D](https://arxiv.org/html/2506.21878v1#A4.Ex109)), we derive that

$$ $\displaystyle|IV.1|=O_{p}(r\alpha_{T}^{-2}\kappa_{k}^{2}).$ (57) $$

By Lemma [5](https://arxiv.org/html/2506.21878v1#Thmtheorem5) and ([D](https://arxiv.org/html/2506.21878v1#A4.Ex109)), we have that

$$ $\displaystyle|IV.2|=O_{p}\Big{(}r^{1/2}\big{\|}\mathbf{P}(\eta_{k+1})-\mathbf{ P}^{\widetilde{\eta}_{k},\widetilde{\eta}_{k+1}}\big{\|}_{\mathrm{F}}\Big{)}=O _{p}(r^{1/2}\alpha_{T}^{-1/2}\kappa_{k}).$ (58) $$

Combining ([D](https://arxiv.org/html/2506.21878v1#A4.Ex137)), ([57](https://arxiv.org/html/2506.21878v1#A4.E57)) and ([58](https://arxiv.org/html/2506.21878v1#A4.E58)), we conclude that

$$ $\displaystyle|IV|=o_{p}(r\kappa_{k}^{2}+r^{1/2}\kappa_{k}).$ (59) $$

Step 3.5. Order of magnitude of $V$.
We now analyze the final term

$$ $\displaystyle V=$ $\displaystyle\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\big{\|}\mathbf{A}(t)-\mathbf{P}( \eta_{k})\big{\|}_{\mathrm{F}}^{2}-\big{\|}\mathbf{A}(t)-\mathbf{P}(\eta_{k+1} )\big{\|}_{\mathrm{F}}^{2}$ $\displaystyle=$ $\displaystyle\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\big{\|}\mathbf{P}(\eta_{k})- \mathbf{P}(\eta_{k+1})\big{\|}_{\mathrm{F}}^{2}-2\sum_{t=\eta_{k}+1}^{\eta_{k} +r}\big{\langle}\mathbf{A}(t)-\mathbf{P}(\eta_{k+1}),\mathbf{P}(\eta_{k})- \mathbf{P}(\eta_{k+1})\big{\rangle}$ $\displaystyle=$ $\displaystyle r\kappa_{k}^{2}-2V.1$ (60) $$

Using Lemma [5](https://arxiv.org/html/2506.21878v1#Thmtheorem5), we obtain that

$$ $\displaystyle|V.1|=O_{p}\Big{(}r^{1/2}\big{\|}\mathbf{P}(\eta_{k})-\mathbf{P}( \eta_{k+1})\big{\|}_{\mathrm{F}}\Big{)}=O_{p}(r^{1/2}\kappa_{k}).$ (61) $$

Step 3.6:
Combining ([38](https://arxiv.org/html/2506.21878v1#A4.E38)), ([44](https://arxiv.org/html/2506.21878v1#A4.E44)), ([50](https://arxiv.org/html/2506.21878v1#A4.E50)), ([55](https://arxiv.org/html/2506.21878v1#A4.E55)), ([59](https://arxiv.org/html/2506.21878v1#A4.E59)), ([D](https://arxiv.org/html/2506.21878v1#A4.Ex139)) and ([61](https://arxiv.org/html/2506.21878v1#A4.E61)) we have for all $r\kappa_{k}^{2}\geq 1$ that

$$ $r\kappa_{k}^{2}=O_{p}(1).$ $$

Step 4. Limiting Distributions. For any $t\in(\tilde{s}_{k},\tilde{e}_{k})$, define

$$ $\widetilde{\mathcal{Q}}_{k}(t)=\sum_{u=\tilde{s}_{k}+1}^{t}\|\mathbf{A}(u)- \mathbf{P}(\eta_{k})\|_{\mathrm{F}}^{2}+\sum_{u=t+1}^{\tilde{e}_{k}}\|\mathbf{ A}(u)-\mathbf{P}(\eta_{k+1})\|_{\mathrm{F}}^{2}.$ $$

Note that the term $V$ defined in ([38](https://arxiv.org/html/2506.21878v1#A4.E38)) satisfies

$$ $V=\widetilde{\mathcal{Q}}_{k}(\eta_{k}+r)-\widetilde{\mathcal{Q}}_{k}(\eta_{k}),$ $$

and hence by ([38](https://arxiv.org/html/2506.21878v1#A4.E38)), ([44](https://arxiv.org/html/2506.21878v1#A4.E44)), ([50](https://arxiv.org/html/2506.21878v1#A4.E50)), ([55](https://arxiv.org/html/2506.21878v1#A4.E55)), ([59](https://arxiv.org/html/2506.21878v1#A4.E59)) and $r\kappa_{k}^{2}=O_{p}(1)$, we have that

$$ $\big{|}\mathcal{Q}_{k}(\eta_{k}+r)-\mathcal{Q}_{k}(\eta_{k})-\big{\{} \widetilde{\mathcal{Q}}_{k}(\eta_{k}+r)-\widetilde{\mathcal{Q}}_{k}(\eta_{k}) \big{\}}\big{|}\leq|I|+|II|+|III|+|IV|\overset{p}{\rightarrow}0.$ $$

Therefore, by Slutsky’s theorem, it suffices to derive the limiting distributions of $\widetilde{\mathcal{Q}}_{k}(\eta_{k}+r)-\widetilde{\mathcal{Q}}_{k}(\eta_{k})$ as $T\to\infty$. We consider the two scenarios for $\kappa_{k}$.

Non-vanishing scenario. Suppose $\kappa_{k}\to\rho_{k}$, as $T\to\infty$, with $\rho_{k}>0$ being an absolute constant. For $r<0$, we have that

$$ $\displaystyle\widetilde{\mathcal{Q}}_{k}(\eta_{k}+r)-\widetilde{\mathcal{Q}}_{ k}(\eta_{k})=$ $\displaystyle\sum_{t=\eta_{k}+r+1}^{\eta_{k}}\Big{\{}\big{\|}\mathbf{A}(t)- \mathbf{P}(\eta_{k+1})\big{\|}_{\mathrm{F}}^{2}-\big{\|}\mathbf{A}(t)-\mathbf{ P}(\eta_{k})\big{\|}_{\mathrm{F}}^{2}\Big{\}}$ $\displaystyle=$ $\displaystyle\sum_{t=\eta_{k}+r+1}^{\eta_{k}}\big{\|}\mathbf{P}(\eta_{k})- \mathbf{P}(\eta_{k+1})\big{\|}_{\mathrm{F}}^{2}$ $\displaystyle\hskip 14.22636pt-2\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\big{\langle} \mathbf{A}(t)-\mathbf{P}(\eta_{k}),\mathbf{P}(\eta_{k+1})-\mathbf{P}(\eta_{k}) \big{\rangle}$ $\displaystyle\xrightarrow{\mathcal{D}}$ $\displaystyle\,-r\rho_{k}^{2}-2\rho_{k}\sum_{t=r+1}^{0}\langle\boldsymbol{ \mathbf{\Psi}}_{k},\mathbf{E}_{k}(t)\rangle,$ (62) $$

with $\boldsymbol{\mathbf{\Psi}}_{k}$ defined in Model [1](https://arxiv.org/html/2506.21878v1#Thmmodel1), and for any $k\in[K+1]$ and $t\in\mathbb{Z}$,
$\mathbf{E}_{k}(t)=\mathbf{A}_{k}(t)-\mathbf{P}(\eta_{k})$ with $\{\mathbf{A}_{k}(t)\}_{t\in\mathbb{Z}}\overset{\mathrm{i.i.d.}}{\sim}\mathrm{
MRDPG}(\{X_{i}\}_{i=1}^{n},\{W_{(l)}(\eta_{k})\}_{l\in[L]})$.

For $r>0$, we have that when $T\to\infty$,

$$ $\displaystyle\widetilde{\mathcal{Q}}_{k}(\eta_{k}+r)-\widetilde{\mathcal{Q}}_{ k}(\eta_{k})=$ $\displaystyle\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\Big{\{}\big{\|}\mathbf{A}(t)- \mathbf{P}(\eta_{k})\big{\|}_{\mathrm{F}}^{2}-\big{\|}\mathbf{A}(t)-\mathbf{P} (\eta_{k+1})\big{\|}_{\mathrm{F}}^{2}\Big{\}}$ $\displaystyle=$ $\displaystyle\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\big{\|}\mathbf{P}(\eta_{k})- \mathbf{P}(\eta_{k+1})\big{\|}_{\mathrm{F}}^{2}$ $\displaystyle\hskip 14.22636pt+2\sum_{t=\eta_{k}+1}^{\eta_{k}+r}\big{\langle} \mathbf{A}(t)-\mathbf{P}(\eta_{k+1}),\mathbf{P}(\eta_{k+1})-\mathbf{P}(\eta_{k })\big{\rangle}$ $\displaystyle\xrightarrow{\mathcal{D}}$ $\displaystyle\,r\rho_{k}^{2}+2\rho_{k}\sum_{t=1}^{r}\langle\boldsymbol{\mathbf {\Psi}}_{k},\mathbf{E}_{k+1}(t)\rangle.$ (63) $$

By Slutsky’s theorem and the argmin continuous mapping theorem (see e.g. Theorem 3.2.2 in Wellner et al., 2013), we obtain

$$ $\widehat{\eta}_{k}-\eta_{k}\xrightarrow{\mathcal{D}}\operatorname*{argmin} \mathcal{P}_{k}(r),$ $$

which completes the proof of part Theorem [3](https://arxiv.org/html/2506.21878v1#Thmtheorem3).

Vanishing scenario. Let $m=\kappa_{k}^{-2}$, noting that $m\to\infty$ as $T\to\infty$. For $r>0$, we have that

$$ $\displaystyle\widetilde{\mathcal{Q}}_{k}(\eta_{k}+rm)-\widetilde{\mathcal{Q}}_ {k}(\eta_{k})$ $\displaystyle=$ $\displaystyle\sum_{t=\eta_{k}+1}^{\eta_{k}+rm}\Big{\{}\big{\|}\mathbf{A}(t)- \mathbf{P}(\eta_{k})\big{\|}_{\mathrm{F}}^{2}-\big{\|}\mathbf{A}(t)-\mathbf{P} (\eta_{k+1})\big{\|}_{\mathrm{F}}^{2}\Big{\}}$ $\displaystyle=$ $\displaystyle\sum_{t=\eta_{k}+1}^{\eta_{k}+rm}\big{\|}\mathbf{P}(\eta_{k})- \mathbf{P}(\eta_{k+1})\big{\|}_{\mathrm{F}}^{2}+2\sum_{t=\eta_{k}+1}^{\eta_{k} +rm}\big{\langle}\mathbf{A}(t)-\mathbf{P}(\eta_{k+1}),\mathbf{P}(\eta_{k+1})- \mathbf{P}(\eta_{k})\big{\rangle}$ $\displaystyle=$ $\displaystyle r+\frac{2}{\sqrt{m}}\sum_{t=\eta_{k}+1}^{\eta_{k}+rm}\big{ \langle}\mathbf{A}(t)-\mathbf{P}(\eta_{k+1}),\boldsymbol{\mathbf{\Psi}}_{k} \big{\rangle}.$ $$

By the functional central limit theorem, we have that when $T\to\infty$,

$$ $\frac{1}{\sqrt{m}}\sum_{t=\eta_{k}+1}^{\eta_{k}+rm}\big{\langle}\mathbf{A}(t)- \mathbf{P}(\eta_{k+1}),\boldsymbol{\mathbf{\Psi}}_{k}\big{\rangle}\xrightarrow {\mathcal{D}}\sigma_{k,k+1}\mathbb{B}_{1}(r),$ $$

where $\mathbb{B}_{1}(r)$ is a standard Brownian motion and for any $k\in[K]$ and $k^{\prime}\in\{k,k+1\}$,
$\sigma_{k,k^{\prime}}^{2}=\operatorname{Var}\big{(}\langle\boldsymbol{\mathbf{
\Psi}}_{k},\mathbf{E}_{k^{\prime}}(1)\rangle\big{)}$.
Consequently, as $T\to\infty$

$$ $\widetilde{\mathcal{Q}}_{k}(\eta_{k}+rm)-\widetilde{\mathcal{Q}}_{k}(\eta_{k}) \xrightarrow{\mathcal{D}}r+2\sigma_{k,k+1}\mathbb{B}_{1}(r).$ $$

Similarly, for $r<0$, we have that when $T\to\infty$

$$ $\widetilde{\mathcal{Q}}_{k}(\eta_{k}+rm)-\widetilde{\mathcal{Q}}_{k}(\eta_{k}) \overset{D}{\to}-r+2\sigma_{k,k}\mathbb{B}_{2}(-r),$ $$

where $\mathbb{B}_{2}(r)$ is a standard Brownian motion.
Applying Slutsky’s theorem and the argmin continuous mapping theorem (see e.g. Theorem 3.2.2 in Wellner et al., 2013), we conclude that

$$ $\kappa_{k}^{2}(\widehat{\eta}_{k}-\eta_{k})\xrightarrow{\mathcal{D}} \operatorname*{argmin}\mathcal{P}_{k}^{\prime}(r),$ $$

which completes the proof of Theorem [2](https://arxiv.org/html/2506.21878v1#Thmtheorem2).

∎

###### Proof.

## Appendix E Additional details and results in Section 4

All experiments were run on a CPU with 16GB RAM. For each synthetic scenario with node size $n=100$, number of layers $L=4$ and time span $T=200$, the compute time is about $10$ hours to localize the change points and to construct the confidence intervals over $100$ Monte Carlo trials. For each real data experiment, the computation time is approximately $15$ minutes to perform change point localization and confidence interval construction.

### E.1 Additional results in Section 4.1

Table [3](https://arxiv.org/html/2506.21878v1#A5.T3) presents the results for Scenarios 3 and 4.

**Table 3: Means of evaluation metrics for networks simulated from Scenarios 3 and 4.**
|  |  | Scenario 3 | Scenario 4 |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $n$ | Method | $|\widehat{K}-K|\downarrow$ | $d(\widehat{\mathcal{C}}|\mathcal{C})\downarrow$ | $d(\mathcal{C}|\widehat{\mathcal{C}})\downarrow$ | $C(\mathcal{G},\mathcal{G^{\prime}})\uparrow$ | $|\widehat{K}-K|\downarrow$ | $d(\widehat{\mathcal{C}}|\mathcal{C})\downarrow$ | $d(\mathcal{C}|\widehat{\mathcal{C}})\downarrow$ | $C(\mathcal{G},\mathcal{G^{\prime}})\uparrow$ |
| $50$ | CPDmrdpg | $0.19$ | $9.64$ | $0.14$ | $95.11\%$ | $0.00$ | $0.02$ | $0.02$ | $99.98\%$ |
| gSeg (nets.) | $0.98$ | Inf | Inf | $68.93\%$ | $5.00$ | Inf | Inf | $0.00\%$ |  |
| kerSeg (nets.) | $0.16$ | $0.18$ | $2.06$ | $98.90\%$ | $0.36$ | $0.14$ | $2.65$ | $98.56\%$ |  |
| gSeg (frob.) | $0.92$ | Inf | Inf | $66.78\%$ | $1.53$ | Inf | Inf | $74.92\%$ |  |
| kerSeg (frob.) | $0.82$ | $48.52$ | $5.11$ | $73.55\%$ | $0.40$ | $0.05$ | $3.71$ | $98.12\%$ |  |
| $100$ | CPDmrdpg | $0.00$ | $0.02$ | $0.02$ | $99.98\%$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |
| gSeg (nets.) | $0.69$ | Inf | Inf | $80.10\%$ | $4.98$ | Inf | Inf | $0.77\%$ |  |
| kerSeg (nets.) | $0.17$ | $0.00$ | $3.26$ | $99.16\%$ | $0.34$ | $0.08$ | $2.93$ | $98.47\%$ |  |
| gSeg (frob.) | $0.79$ | Inf | Inf | $72.11\%$ | $1.86$ | Inf | Inf | $68.57\%$ |  |
| kerSeg (frob.) | $0.79$ | $48.82$ | $4.75$ | $73.80\%$ | $0.42$ | $0.06$ | $2.93$ | $98.63\%$ |  |

We then present a sensitivity analysis of the threshold constant $c_{\tau,1}$ used in Algorithm [1](https://arxiv.org/html/2506.21878v1#alg1). Specifically, based on Theorem [1](https://arxiv.org/html/2506.21878v1#Thmtheorem1), we set the threshold value $\tau=c_{\tau,1}n\sqrt{L}\log^{3/2}(T)$. The constant $c_{\tau,1}$ was chosen empirically assessing the false positive rate in an MSBM model without any change points. The construction follows the details provided in Section [4.1](https://arxiv.org/html/2506.21878v1#S4.SS1), with four evenly-sized communities. We found that $c_{\tau,1}=0.1$ detects a change point approximately $1\%$ of the time in this scenario, demonstrating ideal false positive control. Smaller values of $c_{\tau,1}$, such as $c_{\tau,1}\in\{0.05,0.08\}$, led to more frequent false detections, while larger values $c_{\tau,1}\in\{0.12,0.15,0.20\}$ failed to detect change points in this scenario.

Tables [4](https://arxiv.org/html/2506.21878v1#A5.T4)–[7](https://arxiv.org/html/2506.21878v1#A5.T7) report the results of the sensitivity analysis for Scenarios 1 through 4, varying $c_{\tau,1}\in\{0.05,0.08,0.10,0.12,0.15,0.20,0.25\}$. These results demonstrate that our proposed method is relatively robust against the choices of $c_{\tau,1}$.

**Table 4: Means of evaluation metrics for dynamic networks simulated from Scenario 1, varying $c_{\tau,1}$.**
| $n$ | $c_{\tau,1}$ | $|\widehat{K}-K|\downarrow$ | $d(\widehat{\mathcal{C}}|\mathcal{C})\downarrow$ | $d(\mathcal{C}|\widehat{\mathcal{C}})\downarrow$ | $C(\mathcal{G},\mathcal{G^{\prime}})\uparrow$ |
| --- | --- | --- | --- | --- | --- |
| $50$ | $0.25$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |
| $0.20$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |  |
| $0.15$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |  |
| $0.12$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |  |
| $0.10$ | $0.01$ | $0.00$ | $0.42$ | $99.86\%$ |  |
| $0.08$ | $0.25$ | $0.00$ | $6.68$ | $97.80\%$ |  |
| $0.05$ | $5.18$ | $0.00$ | $52.86$ | $67.50\%$ |  |
| $100$ | $0.25$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |
| $0.20$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |  |
| $0.15$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |  |
| $0.12$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |  |
| $0.10$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |  |
| $0.08$ | $0.15$ | $0.00$ | $4.98$ | $98.54\%$ |  |
| $0.05$ | $5.02$ | $0.00$ | $53.84$ | $67.56\%$ |  |

**Table 5: Means of evaluation metrics for dynamic networks simulated from Scenario 2, varying $c_{\tau,1}$.**
| $n$ | $c_{\tau,1}$ | $|\widehat{K}-K|\downarrow$ | $d(\widehat{\mathcal{C}}|\mathcal{C})\downarrow$ | $d(\mathcal{C}|\widehat{\mathcal{C}})\downarrow$ | $C(\mathcal{G},\mathcal{G^{\prime}})\uparrow$ |
| --- | --- | --- | --- | --- | --- |
| $50$ | $0.25$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |
| $0.20$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |  |
| $0.15$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |  |
| $0.12$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |  |
| $0.10$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |  |
| $0.08$ | $0.02$ | $0.00$ | $0.64$ | $99.68\%$ |  |
| $0.05$ | $3.79$ | $0.00$ | $28.46$ | $75.43\%$ |  |
| $100$ | $0.25$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |
| $0.20$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |  |
| $0.15$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |  |
| $0.12$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |  |
| $0.10$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |  |
| $0.08$ | $0.05$ | $0.00$ | $1.14$ | $99.38\%$ |  |
| $0.05$ | $3.53$ | $0.00$ | $28.60$ | $76.50\%$ |  |

**Table 6: Means of evaluation metrics for dynamic networks simulated from Scenario 3, varying $c_{\tau,1}$.**
| $n$ | $c_{\tau,1}$ | $|\widehat{K}-K|\downarrow$ | $d(\widehat{\mathcal{C}}|\mathcal{C})\downarrow$ | $d(\mathcal{C}|\widehat{\mathcal{C}})\downarrow$ | $C(\mathcal{G},\mathcal{G^{\prime}})\uparrow$ |
| --- | --- | --- | --- | --- | --- |
| $50$ | $0.25$ | $1.00$ | $50.00$ | $0.00$ | $75.00\%$ |
| $0.20$ | $0.96$ | $48.00$ | $0.00$ | $76.00\%$ |  |
| $0.15$ | $0.64$ | $32.00$ | $0.00$ | $84.00\%$ |  |
| $0.12$ | $0.39$ | $19.58$ | $0.08$ | $90.17\%$ |  |
| $0.10$ | $0.19$ | $9.64$ | $0.14$ | $95.11\%$ |  |
| $0.08$ | $0.09$ | $4.30$ | $0.50$ | $97.61\%$ |  |
| $0.05$ | $4.27$ | $0.36$ | $32.54$ | $71.55\%$ |  |
| $100$ | $0.25$ | $0.43$ | $21.50$ | $0.00$ | $89.25\%$ |
| $0.20$ | $0.15$ | $7.52$ | $0.02$ | $96.23\%$ |  |
| $0.15$ | $0.00$ | $0.02$ | $0.02$ | $99.98\%$ |  |
| $0.12$ | $0.00$ | $0.02$ | $0.02$ | $99.98\%$ |  |
| $0.10$ | $0.00$ | $0.02$ | $0.02$ | $99.98\%$ |  |
| $0.08$ | $0.06$ | $0.02$ | $1.08$ | $99.57\%$ |  |
| $0.05$ | $4.04$ | $0.02$ | $32.76$ | $73.62\%$ |  |

**Table 7: Means of evaluation metrics for dynamic networks simulated from Scenario 4, varying $c_{\tau,1}$.**
| $n$ | $c_{\tau,1}$ | $|\widehat{K}-K|\downarrow$ | $d(\widehat{\mathcal{C}}|\mathcal{C})\downarrow$ | $d(\mathcal{C}|\widehat{\mathcal{C}})\downarrow$ | $C(\mathcal{G},\mathcal{G^{\prime}})\uparrow$ |
| --- | --- | --- | --- | --- | --- |
| $50$ | $0.25$ | $2.67$ | $83.20$ | $0.00$ | $62.47\%$ |
| $0.20$ | $1.19$ | $28.40$ | $0.00$ | $85.63\%$ |  |
| $0.15$ | $0.13$ | $2.60$ | $0.00$ | $98.67\%$ |  |
| $0.12$ | $0.01$ | $0.22$ | $0.02$ | $99.88\%$ |  |
| $0.10$ | $0.00$ | $0.02$ | $0.02$ | $99.98\%$ |  |
| $0.08$ | $0.00$ | $0.02$ | $0.02$ | $99.98\%$ |  |
| $0.05$ | $0.75$ | $0.02$ | $11.94$ | $93.36\%$ |  |
| $100$ | $0.25$ | $0.01$ | $0.20$ | $0.00$ | $99.90\%$ |
| $0.20$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |  |
| $0.15$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |  |
| $0.12$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |  |
| $0.10$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |  |
| $0.08$ | $0.00$ | $0.00$ | $0.00$ | $100\%$ |  |
| $0.05$ | $0.89$ | $0.00$ | $12.46$ | $92.55\%$ |  |

### E.2 Additional details and results in Section 4.2

This section provides a detailed analysis of the U.S. air transportation network data, evaluates the performance of competing methods (introduced in Section [4.1](https://arxiv.org/html/2506.21878v1#S4.SS1)) on both real datasets and presents the constructed confidence intervals using the procedure in Section [3.1](https://arxiv.org/html/2506.21878v1#S3.SS1).

The U.S. air transportation network data consist of monthly data from January 2015 to June 2022 ($T=90$) and are available from Bureau of Transportation Statistics (2022). Each node corresponds to an airport and each layer represents a commercial airline. A directed edge in a given layer indicates a direct flight operated by a specific commercial airline between two airports. We choose the $L=4$ airlines with the highest flight volumes and the $n=50$ airports with the most departures and arrivals. Our method identifies change points in December 2015, June 2017, February 2019, February 2020 and February 2021, all corresponding to major abruptions in the U.S. aviation industry.

The change point in December 2015 coincides with increased regulatory scrutiny over airline consolidation, following concerns raised by the American Antitrust Institute about reduced market competition after a series of mergers. The June 2017 change point aligns with the proposal of the Aviation Innovation, Reform and Reauthorization Act, which advocated for privatizing air traffic control and influenced route planning among carriers. Moreover, the February 2019 change point follows the U.S. government shutdown (December 2018 - January 2019), which caused Transportation Security Administration staffing shortages and significant operational disruptions, prompting stabilization efforts in the months that followed. Lastly, the most significant structural disruptions emerged in February 2020 and February 2021, aligning with the initial shock and continued fallout of the COVID-19 pandemic, which triggered widespread flight cancellations, demand collapse and structural reconfiguration in the aviation industry.

Performance of competitors.
Table [8](https://arxiv.org/html/2506.21878v1#A5.T8) summarizes the change points detected by the proposed and competing methods
for the worldwide agricultural trade network data. Notably, the gSeg method fails to detect any change points after 2010, regardless of input type. Meanwhile, the kerSeg method detects change points in 1990 and 1992, which are temporally too close. In contrast, our proposed method (CPDmrdpg) identifies four major change points that align well with known geopolitical and policy-related events.

Table [9](https://arxiv.org/html/2506.21878v1#A5.T9) presents the results for the U.S. air transportation network data. Although the kerSeg method using networks as input demonstrates a good performance in the simulation study, it detects an excessive number of change points in this real data experiment, making the results unreliable and raising concerns about false positives. Similarly, the kerSeg method that uses layer-wise Frobenius norms as input has detected change points that are too close, yielding clusters of change points that could potentially be grouped together. On the contrary, the gSeg method that uses the Frobenius norms as input detects too few change points, while the gSeg method using networks as input has detected too many change points. The proposed CPDmrdpg method (Algorithm [1](https://arxiv.org/html/2506.21878v1#alg1)) yields five change points that align well with known disruptions and policy changes in the aviation sector.

While the competitor methods do detect important and relevant change points in both two real datasets, they tend to either over- or under-segment the time span. These patterns suggest that the change points identified by the competing methods may be less realistic or informative compared to those identified by the proposed method.

**Table 8: Detected change points for the worldwide agricultural trade network data.**
| Method | Detected change points |
| --- | --- |
| CPDmrdpg | 1991, 1999, 2005, 2013 |
| gSeg (nets.) | 1993, 2002, 2010 |
| kerSeg (nets.) | 1990, 1992, 1999, 2005, 2012 |
| gSeg (frob.) | 1993, 2002, 2009 |
| kerSeg (frob.) | 1990, 1992, 1997, 2003, 2012 |

**Table 9: Detected change points for the U.S. air transportation network data.**
| Method | Detected change points |
| --- | --- |
| CPDmrdpg | 2015-12, 2017-06, 2019-02, 2020-02, 2021-02 |
| gSeg (nets.) | 2015-11, 2016-10, 2017-09, 2018-09, 2019-09, 2020-10, 2021-08 |
| kerSeg (nets.) | 2015-11, 2016-03, 2016-10, 2017-05, 2017-09, 2018-05, 2018-10 |
|  | 2019-03, 2019-09, 2020-03, 2020-10, 2021-03, 2021-09 |
| gSeg (frob.) | 2015-11, 2020-01, 2021-03 |
| kerSeg (frob.) | 2015-11, 2017-10, 2020-01, 2021-03, 2021-05, 2021-09, 2022-01 |

Performance of constructed confidence intervals.
Tables [10](https://arxiv.org/html/2506.21878v1#A5.T10) and [11](https://arxiv.org/html/2506.21878v1#A5.T11) report the detected change point from Algorithm [1](https://arxiv.org/html/2506.21878v1#alg1) and the $95\%$ confidence intervals constructed via the procedure from Section [3.1](https://arxiv.org/html/2506.21878v1#S3.SS1), for the agricultural trade and air transportation networks, respectively.

**Table 10: Detected change point from Algorithm [1](https://arxiv.org/html/2506.21878v1#alg1) and $95\%$ confidence intervals via Section [3.1](https://arxiv.org/html/2506.21878v1#S3.SS1) for the worldwide agricultural trade network data.**
| Detected change points | Time point | Confidence interval |
| --- | --- | --- |
| 1991 | $6$ | $(5.97,6.03)$ |
| 1999 | $14$ | $(13.98,14.02)$ |
| 2005 | $20$ | $(17.97,18.05)$ |
| 2013 | $28$ | $(25.99,26.06)$ |

**Table 11: Detected change point from Algorithm [1](https://arxiv.org/html/2506.21878v1#alg1) and $95\%$ confidence intervals via Section [3.1](https://arxiv.org/html/2506.21878v1#S3.SS1) for the U.S. air transportation network data.**
| Detected change points | Time point | Confidence interval |
| --- | --- | --- |
| 2015-12 | $12$ | $(11.55,12.41)$ |
| 2017-06 | $30$ | $(28.79,30.98)$ |
| 2019-02 | $50$ | $(49.67,53.22)$ |
| 2020-02 | $62$ | $(59.66,60.36)$ |
| 2021-02 | $74$ | $(73.58,74.27)$ |

## Peer Reviews

confidence: 2
contribution: 3
initial_rating: 6
presentation: 3
questions: Please refer to the question in the Weakness section.
soundness: 3
strengths: 1. The theoretical contributions are novel, deriving limiting distributions for change point estimators in network data. The localization rates match minimax optimality for single-layer networks while extending to multilayer settings.

2. The paper provides consistency results, limiting distributions, a data-driven confidence interval procedure, and extensive robustness checks.
summary: This paper proposes a two-stage algorithm for offline change point localization in dynamic multilayer random dot product graphs. Stage I uses seeded binary segmentation with CUSUM statistics; Stage II refines estimates via tensor heteroskedastic PCA. The authors establish consistency for change point detection/localization, derive limiting distributions under vanishing and non-vanishing jump regimes, and provide a data-driven confidence interval procedure. Experiments on synthetic and real networks demonstrate strong performance.
weaknesses: While online change point detection in multilayer networks has been recently studied (Wang et al., 2025), the paper does not clearly articulate what specific technical challenges arise in the offline setting or why such an extension from online methods is non-trivial.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 2
questions: - Middle of page 3: “We allow all model parameters … to diverge with T”. Please clarify this, how do they relate to T exactly?
soundness: 3
strengths: The paper is original in that they consider a model that has not been studied before. This model is well motivated from the literature and appears in practice. To the best of their (and my) knowledge, this is the first result of its kind in the context of dynamic network data. 

The paper is well written, they state the scalings very clearly which is important in the area of change points. Their theory is sound, and they use standard techniques/sub-methods from the change point literature such as CUSUM statistics. They also acknowledge potential improvements in the conclusion section which are sound. 

The contribution is significant to some extent, as the experiments indicate an experimental improvement over previous models that seems to come from the fact that they are looking at a more specific problem than competing methods (some competing methods are nonparametric). However, I must mention that the techniques are not in themselves novel (at least at first glance), but it would be good if the authors could comment on this.
summary: The paper proposes and analyzes a new moethod for offline change point localization and inference in dynamic multilayer networks (following the D-MRDPGs model). They prove consistency of their estimator, and derive limiting distributions of the refined estimators in general regimes. They provide methods for computing confidence intervals, and demonstrate the favourable performance of their estimator on both synthetic and real data.
weaknesses: From my understanding, competing methods only work change point detection in single layer dynamic networks, or they are inherently nonparametric methods. Can these be adapted to the multilayer setting? Can you comment on how fair this comparison is? It is not clear to me at this point. 

Also, implications from the experimental results could be made more clear, for example, by including figures in the main body of the text. You state that “across all scenarios, our method achieves best performance” but that does not seem to be true, as if I dig into Table 3 in Appendix F.1, it seems that in the n=50 case you do not clearly outperform kerSeg (nets.). Could you comment on this, and at least make it clear in the paper that you do not outperform in all cases, and maybe some intuition as to why? Showing some of these results graphically also would help for clarity. 

The comparison with other methods on real world data should be made more clear (in Appendix F.2 for example), and brought into the main paper. This is very important, and can help better clarify the benefit of your method, which is currently not so clear.

---

confidence: 4
contribution: 2
initial_rating: 4
presentation: 3
questions: 1. Comparison with $\ell_0$ approaches. Please include a discussion comparing against spacing-free $\ell_0$ methods. Proper references should be cited and compared; I list some, but I think more can likely be found.

2. I think including a vectorized version of multiplayer networks and then using binary segmentation is worth being included as a baseline to see how much the improvement in stage 2 is. 

3. To demonstrate robustness beyond the current (near-)even spacing and global-weight-jump setting, please add stress tests that randomly change times i.i.d. from a distribution. 


The paper has very strong assumptions, and the model is not very applicable given those assumptions. However, the paper is technically solid, and the TH-PCA refinement plus distributional results are valuable. I would not oppose acceptance if proper questions are addressed or at least discussed.
soundness: 2
strengths: The paper is clearly written; the model, assumptions, and the two-stage algorithm are easy to follow.

Clean multilayer insight: the Stage-2 projection/matched filter makes precise how low rank and layer pooling provide $\sqrt{L}$ SNR gains.

Technical development is careful (ensuring consistency and limiting distribution for the refined estimator); a standardized pipeline (screen then refine) is appealing.
summary: This paper develops a two-stage CUSUM pipeline for change-point detection in dynamic multilayer random dot product graphs. Stage-1 performs seeded/binary segmentation with a tensor CUSUM; Stage-2 applies a low-rank TH-PCA/matched-filter refinement that exploits shared latent positions across layers, yielding a $\sqrt{L}$ SNR gain and $O(\kappa^{-2})$ localization rate under large minimal spacing and global layer-weight changes. Proofs establish consistency and limit laws; simulations and small real data are provided.
weaknesses: 1. $\textbf{Scope (extension vs.~novelty):}$ The contribution appears closely aligned with prior multilayer RDPG pipelines [1] that (i) \emph{screen} for candidate changes via CUSUM/binary segmentation–type statistics and then (ii) refine via a low-rank projection step (and the target is similar). In the current draft, the primary difference seems to be a regime shift (online $\rightarrow$ offline) while keeping the same modeling assumptions (multilayer RDPG/GRDPG) and the same screen–then–refine structure and target. As written, this reads more like an extension to the offline setting than a fully new methodological framework.


2. $\textbf{Narrow regime.}$ The methodology focuses on time-invariant latent positions with only changes in global layer weights. This means that the communities of nodes remain consistent over time; in this scenario, the change points appear to be narrowly defined, as they cannot disrupt the community structure of the networks. Many practical settings involve node-localized changes alongside shifts in community structure, where the latent positions $X(t)$ should also be able to vary (See [2] for an example). The current scope does not account for these regimes. 


3. $\textbf{Spacing dependence.}$ The paper assumes a strong gap at the order $T$--with a gap at the order of $T$ and the number of change points is $O(1)$, spacing-free $\ell_0$ approach can already achieve $\textbf{exact selection}$ and the same $O(\kappa^{-2})$ localization approach without a minimal spacing condition. In contrast, using binary segmentation requires setting up the minimal spacing window, which is not necessary and could be misleading. For example. The paper effectively must generate evenly spaced change points. If, instead, change times are drawn i.i.d.\ from a distribution—e.g.,
$
\eta_1,\ldots,\eta_K \sim \mathrm{Unif} \{ 2,\ldots,T \},
$
the minimal spacing
$
\Delta_{\min}:=\min_{1\le k\le K+1}\bigl(\eta_k-\eta_{k-1}\bigr)
$
typically collapses. This destroys the large-spacing assumption required by a CUSUM screen and makes any fixed spacing window ill-posed. By contrast, spacing-free $\(\ell_0\)$ approaches remain optimal/consistent in this setting and do not require enforcing even spacing.


For the $\ell_0$ approach to estimation of dynamic (possible multilayer) networks, I refer to [2], but I believe more literature can be found.

For comparison between CUSUM and $\ell_0$ approach, I refer to [3]. See Section 4, 'However, compared to the $l_0$-penalization methods, WBS is computationally more expensive and involves more tuning parameters.''.

4. Simulations use (nearly) evenly spaced changes, reinforcing the strong-spacing assumption; real-data experiments also use very small $n,T,L$, This leaves the reader with the impression that the method hinges on large spacing and may not scale or adapt to realistic, short-gap or node–localized scenarios.

[1] Wang, Fan, et al. "Multilayer random dot product graphs: Estimation and online change point detection." Journal of the Royal Statistical Society Series B: Statistical Methodology (2025): qkaf051.

[2] Zhao, P., Bhattacharya, A., Pati, D., & Mallick, B. K. (2022). Factorized fusion shrinkage for dynamic relational data. arXiv preprint arXiv:2210.00091.

[3] ''Wang, Daren, Yi Yu, and Alessandro Rinaldo. "Univariate mean change point detection: Penalization, CUSUM and optimality." Electronic Journal of Statistics 14 (2020): 1917-1961.


Decision: Accept (Poster)