A synthetic data approach for FDR control in change-point detection

Change-point analysis, a process of detecting structural changes in a data sequence, has been an active area of research and attracted increasing attention with the growing availability of temporal data. It has applications across a wide range of fields, including but not limited to environmental sciences, econometrics, biology, geosciences, and linguistics. In this context, the accurate and efficient detection of multiple change-points (MCP) is undoubtedly one of the most crucial issues. For example, Liu et al.^[1] introduced a general framework for high-dimensional change-point detection by constructing a U-statistic-based cumulative sum matrix $ \mathcal{C} $ and aggregating it based on the adjusted $ L_{p} $-norm, while Liu et al.^[2] focused on high-dimensional linear models, providing asymptotic validity and an extension to multiple change points via binary segmentation. For more comprehensive reviews of various existing approaches to MCP inference, see Aue & Horváth^[3], and Niu et al.^[4].

However, obtaining consistent estimators for the number and locations of MCP typically requires stringent conditions on the magnitude of changes, as extensively documented in prior studies^[5–13]. Unfortunately, such requirements are often unrealistic, as small change magnitudes tend to cause underfitting. Consequently, an overfitted selection set is often obtained via some conservative algorithms. Furthermore, the empirical performance of certain detection methods is intrinsically linked to the choice of tuning parameters, which requires access to unavailable population-level information. These two issues, referred to as the unreliability of assumptions and the unreliability of algorithms, can introduce false discoveries, potentially leading to the reproducibility crisis if an excessive number of false detections occur.

To tackle this problem, a natural solution is to detect the active set while controlling the false discovery rate (FDR) at a pre-specified level. A widely adopted strategy is to treat change point detection as a multiple hypothesis testing problem, utilizing classical p-value-based methods ^[14–18] to control the FDR. Notable works include Hao et al.^[19], Li et al.^[20], and Cheng et al.^[21]. These methods work properly for the univariate mean change problem, but extending them to a multi-dimensional setting is challenging, since the model's complexity renders the derivation of p-values intractable. Leveraging the knockoff framework^[22–26], a related study is Liu et al.^[27], which proposed a generalized knockoff procedure (GKnockoff) to control the FDR for detecting structural changes in the coefficients of a linear regression model. More recently, Du et al.^[28], and Dai et al.^[29] proposed a data-splitting-based FDR control framework, which outperforms knockoff methods in power under moderate to strong dependence and is more robust than asymptotic p-value-based methods. Chen et al.^[30] adopted this data-splitting philosophy and proposed a data-driven selection procedure for MCP detection while controlling the error rate. This approach is quite general and can handle complex MCP scenarios. However, on the one hand, the symmetry property of the proposed statistic depends on the sample size. If the sample size is too small, the comparison statistic may become asymmetric, distorting the FDR control. On the other hand, data-splitting inevitably reduces the power of change-point detection as only half of the information is utilized.

The limitations of data-splitting motivate us to develop a novel framework for controlling the error rates in high-dimensional MCP detection. This framework, termed the Synthetic Data Filter (SD filter), is designed to enhance the accuracy of multiple change-point identification. Basically, the SD filter procedure consists of three steps. First, the dataset is divided according to the temporal order's parity, as inspired by Chen et al.^[30]. The change-point detection is performed on the odd-part and can be carried out using the adaptive $ \ell_q $ aggregated CUSUM-type statistic introduced by Liu et al.^[1]. Next, information from the odd part is leveraged to generate a synthetic dataset. In the final step, the synthetic dataset is merged with the reserved even-part dataset to construct symmetric statistics and control the false discovery rate (FDR). Figure 1 presents a flow chart summarizing the procedure of the proposed SD filter. The SD filter exhibits the capability to asymptotically control the FDR at the desired level while achieving a power that approaches one under mild conditions. Moreover, this work demonstrates its competence in addressing high-dimensional MCP problems and represents, to the best of current knowledge, the first attempt at controlling the error rate in this context.

Figure 1. 

Flow chart of SD filter.

The synthetic data framework clearly departs from the aforementioned error rate control procedures. Unlike conventional methods that quantify the distribution of the test statistic, the Gaussian multiplier bootstrap is employed ^[31,32] to construct a mirror statistic for FDR control. This framework is highly general, easily extendable to various model settings, and applicable to any scenario requiring data splitting for simultaneous hypothesis formulation and testing. By integrating the synthetic random sample with the second dataset in the testing phase, the original sample size is maintained, thereby boosting statistical power.

Suppose that a sequence of independent data have been observed, $ \mathcal{Z} = \{ \mathbf{z}_1, \ldots, \mathbf{z}_{2n}\} $ collecting from

where, K_n is the number of change-points which could diverge with sample size $ n $ and $ \tau_{k}'s $ are the change locations with the convention that $ \tau_0 = 0 $ and $ \tau_{K_n+1} =2n $. $ F(\cdot | \boldsymbol{\beta}_k) $ represents the model structure of segment k, where $ \boldsymbol{\beta}_k $ is a d-dimensional parameter vector of interest, satisfying $ \boldsymbol{\beta}_{k} \ne \boldsymbol{\beta}_{k+1} $. The setting of MCP in Eq. (1) is quite general and encompasses many classical models, such as the multivariate mean change model and the regression model with structural breaks ^[33].

The objective of this study is to detect the active set $ \mathcal{S} = \{\tau_k, k = 1, \ldots, K_n\} $ while controlling the FDR. However, the definition of FDR in this setting is context-specific, as the optimal rate of change-point estimation is characterized by O_p(1)^[34]. Therefore, it is essential to first establish the corresponding concepts in this context. Given a candidate change-point set $ \hat{\mathcal{S}} = \{ \hat{\tau}_k, k =1,\ldots, \hat{K}_n\} $, the definition of false discovery in MCP detection is as follows:

Definition 2.1 (False discovery). The candidate change-point $ \hat{\tau}_{k} \in \hat{\mathcal{S}} $ is a False Discovery if there is no true change-point falls in

where, $ \hat{\tau}_0 = 0 $ and $ \hat{\tau}_{ \hat{K}_n} = 2n $.

Further, let the set $ \mathcal{I}_0 $ encompass all the false discoveries. Then, the set $ \mathcal{I}_1 = \hat{\mathcal{S}}\cap \mathcal{I}_0^c $ contains all the true discoveries in the selection set $ \hat{\mathcal{S}} $. Throughout this paper, $ \mathcal{I}_1 $ and $ \mathcal{I}_0 $ are referred to as the informative and uninformative sets, respectively. This definition is well-defined, as each candidate change point is unambiguously classified as either a true or false discovery, with no overlap between the two categories. Based on this, the corresponding FDR is defined as follows:

where $ \mathcal{T} $ represents a subset of $ \hat{\mathcal{S}} $ yielded by the selection procedure, and |A| denotes the cardinality of set A. The FDR is the expected value of the False Discovery Proportion (FDP), which represents the ratio of false discoveries to the total number of discoveries.

Synthetic data generating process

This subsection first introduces a synthetic data generating procedure. Following the order-preserving splitting procedure in Zou et al.^[33], the data $ \mathcal{Z} $ is partitioned into odd and even parts

On the subset $ \mathcal{Z}^{O} $, a candidate set of change points is estimated $ \hat{\mathcal{S}} = \{ \hat{\tau}_1, \ldots, \hat{\tau}_{ \hat{K}_n}\} $ using a suitable detection algorithm, such as the aggregated CUSUM method. This training phase allows for the possibility that $ \hat{\mathcal{S}} $ may overestimate the true set of change-points $ \mathcal{S} $. Based on the candidate set $ \hat{\mathcal{S}} $, sets G_k are defined according to Definition 2.2. Then the odd sample $ \mathcal{Z}^O $ and the even sample $ \mathcal{Z}^E $ are partitioned into segments $ \mathcal{Z}_{G_k}^O := \{ \mathbf{z}_{2i-1}: i \in G_k\} $ and $ \mathcal{Z}_{G_k}^E := \{ \mathbf{z}_{2i}: i \in G_k\} $, respectively. In the validation stage, $ \mathcal{Z}^{O} $ is treated as given to avoid dealing with intractable post-selection inference.

Let $ \ell(\boldsymbol{\beta}; \mathbf{z}_i) $ be a suitable loss function evaluated at data point $ \mathbf{z}_i $, with its derivative denoted as $ \mathbf{s}_{ \boldsymbol{\beta}}(\mathbf{z}_i) = {\partial \ell(\boldsymbol{\beta}; \mathbf{z}_i)}/{\partial \boldsymbol{\beta}} $. Ideally, for a given d-dimensional reference vector $ \boldsymbol{\gamma} $, $ \text*{E}\{ \mathbf{s}_{\gamma}(\mathbf{z}_i)\} \ne \text*{E}\{ \mathbf{s}_{ \boldsymbol{\gamma}}(\mathbf{z}_{i^\prime})\} $ is expected when there is a change between $ i $ and $ i^\prime $, since the score function remains constant in regions without change. This motivates the decomposition of the score as

where, $ \boldsymbol{\mu}_i = \text*{E}[\mathbf{s}_{ \boldsymbol{\gamma}}(\mathbf{z}_i)] $ is the expected score at $ \mathbf{z}_i $, and $ \boldsymbol{\zeta}_i = \mathbf{s}_{ \boldsymbol{\gamma}}(\mathbf{z}_i) - \text*{E}[\mathbf{s}_{ \boldsymbol{\gamma}}(\mathbf{z}_i)] $ represents the residual. It is further assumed that $ \text*{Cov}(\boldsymbol{\zeta}_i) = \boldsymbol{\Sigma}^{(k)} $ for all $ i $ such that $ \tau_{k}+1 \leqslant i \leqslant \tau_{k+1} $. This score-type framework was first introduced by Zou et al.^[33] for selecting the number of change-points. The procedure is often invariant to the choice of $ \boldsymbol{\gamma} $, which can therefore be set as $ \boldsymbol{\gamma} := \arg\min_{ \boldsymbol{\beta}} \sum_{ \mathbf{z}_i \in \mathcal{Z}^{O}} \ell(\boldsymbol{\beta}; \mathbf{z}_i) $ when no prior information is available. Given a specific $ \boldsymbol{\gamma} $, the score function $ \mathbf{s}_{ \boldsymbol{\gamma}}(\mathbf{z}_i) $ is denoted as $ \mathbf{s}_i $ for simplicity.

To monitor the change magnitude in the dataset$ \mathcal{Z}_{G_k}^E $, the CUSUM statistic is employed for score $ \mathbf{s}_i^E $ in $ \mathcal{Z}_{G_k}^E $, which is defined as,

where, $ n_k = |G_k| $ and $ s $ varies in $ G_k $ with the exception of points that are too close to both ends. Substitute $ \mathbf{s}^E_i $ by $ \boldsymbol{\mu}_i^E + \boldsymbol{\zeta}_i^E $, it becomes

where, $ \Delta_k(s) = \sqrt{{s(n_k-s)}/{n_k}}\left({s^{-1}} \sum_{i \leqslant s, i \in G_k} \boldsymbol{\mu}_{i}^E-(n_k-s)^{-1} \sum_{i \gt s, i \in G_k} \boldsymbol{\mu}_{i}^E\right). $ If $ \hat{\tau}_k $, located at the s-th position of G_k, is a false discovery, then $ \Delta_k(s) = 0 $. Hence the CUSUM statistic $ \mathbf{c}_k $ is just a combination of the errors,

Observing that for a certain false discovery $ \hat{\tau}_k $ lies in the s-th position of G_k, $ \mathbf{c}_{k}(s) $ consists of two parts: $ s^{-1} \sum_{i \leqslant s, i \in G_k} \boldsymbol{\zeta}^E_{i} $ and $ {(n_k - s)^{-1}}\sum_{i \gt s, i \in G_k} \boldsymbol{\zeta}^E_{i} $. According to the central limit theorem, it can be inferred that the two scaled summations $ s^{-1/2} \sum_{i \leqslant s, i \in G_k} \boldsymbol{\zeta}^E_{i} $ and $ (n_k - s)^{-1/2} \sum_{i \gt s, i \in G_k} \boldsymbol{\zeta}^E_{i} $ converge to a d-dimensional normal distribution with mean zero and covariance matrix $ \boldsymbol{\Sigma}^{(k)} $, provided that both s and $ n_k-s $ are sufficiently large.

Based on this intuition, i.i.d random variables $ \xi_1, \ldots, \xi_{n_k} $ are introduced, each of which follows a standard normal distribution $ N(0,1) $. These variables are also independent of the original dataset $ \mathcal{Z} $ and define two partial sums based on $ \mathcal{Z}^O $ as $ \bar{ \mathbf{s}}_k^{O,-}(s) = s^{-1}\sum_{i \leqslant s, i \in G_k} \mathbf{s}^O_i $ and $ \bar{ \mathbf{s}}_k^{O,+}(s) = (n_k-s)^{-1} \sum_{i \gt s, i \in G_k} \mathbf{s}^O_{i} $. Then a synthetic dataset $ \tilde{ \mathcal{Z}}_{G_k} $ is generated, consisting of the following two parts.

Definition 2.2. Define the synthetic data based on the training sample as :

By using the synthetic dataset $ \tilde{ \mathcal{Z}}_{G_k} $, a synthetic CUSUM can further be constructed

Given the dataset $ \mathcal{Z}_{G_k}^O $, two summations in the $ \tilde{ \mathbf{c}}_{k}(s) $ are also normally distributed, that is

where,

are two plausible estimates of $ \boldsymbol{\Sigma}^{(k)} $. Hence, the synthetic CUSUM statistic $ \tilde{\mathbf{c}}_k(s) $ mimics the distributional behavior of the original CUSUM statistic $ \mathbf{c}_k(s) $ when there are no change points, and it is independent of $ \mathbf{c}_k(s) $ given the odd samples $ \mathcal{Z}^O $.

To conclude, for each candidate change point $ \hat{\tau}_k $, if it is a false discovery, the distributions of $ \mathbf{c}_k(s) $ and $ \tilde{ \mathbf{c}}_{k}(s) $ are highly similar. In contrast, if $ \hat{\tau}_k $ corresponds to a true change point, the two distributions differ significantly. This distinction is visually demonstrated in Fig. 2, where the distributions align closely in the absence of a change point but diverge otherwise. The synthetic data thus serves as a diagnostic tool for detecting false discoveries, thereby enhancing selection power in the validation step.

Figure 2. 

Densities of the synthetic-data-based statistic and the original statistic under stationary (left) and non-stationary (right) settings. Detailed settings are described in the simulation study.

Moreover, it is worth noting that solely comparing the difference at a fixed point $ \hat{\tau}_k $ may not yield the optimal result because $ \hat{\tau}_k $ maximizes the difference in the odd part dataset $ \mathcal{Z}^O $ instead of the even part dataset $ \mathcal{Z}^E $. When the candidate change-points in $ \mathcal{Z}^E $ are validated, it is more reasonable to iterate through all possible locations in the interval G_k and select the point with the largest absolute difference. Furthermore, the CUSUM statistics are aggregated via a $ \ell_{q} $ norm, where q takes value in {1, 2, ..., $ \infty\ $} to adapt the change structure under the alternatives. Hence, the comparison statistics for the original dataset $ \mathcal{Z}^E_{G_k} $ and the synthetic dataset $ \tilde{ \mathcal{Z}}_{G_k} $are defined as:

where, for a vector $ \mathbf{x} \in \mathbb{R}^d $, $ \|\mathbf{x}\|_q := \big(\sum_{j=1}^d |x_j|^q\big)^{1/q} $ for $ 1 \leqslant q \lt \infty $, $ \|\mathbf{x}\|_{\infty} := \max_{1 \leqslant j \leqslant d}|x_j| $ and $ G_k^* $ is a truncated version of $ G_k $, obtained by removing the first $ \underline{s} $ and last $ \underline{s} $ indices. This truncation is necessary since $ \widehat{\boldsymbol{\Sigma}}^{(k)-}(s) $ and $ \widehat{\boldsymbol{\Sigma}}^{(k)+}(s) $ may not be reliable estimates of $ \boldsymbol{\Sigma}^{(k)} $ when sample sizes within the left or right segments are too small. Therefore, if $ \hat{\tau}_k $ is a false discovery, the distribution of $ \tilde{T}_k^q $ is expected to closely approximate that of $ T_k^q $, given the similarity between the synthetic and validation data.

However, if the candidate change-point $ \hat{\tau}_k $ belongs to the informative set, i.e., $ \hat{\tau}_k \in \mathcal{I}_1 $, then $ \Delta_k(s) $ in Eq. (6) should be the leading term in the $ \mathbf{c}_{k}(s) $. Denote

then in this case,

where $ \mathbf{r}_{k}(s), s \in G_k^* $ are all zero mean vectors. Furthermore, the CUSUM statistics defined on the synthetic dataset $ \tilde{ \mathcal{Z}}_{G_k} $ also have zero mean, as they are transformations of zero-mean random vectors. Hence, the original comparison statistic $ T_k^q $ should be sufficiently larger than the synthetic comparison statistic $ \tilde{T}_k^q $. That is,

provided that the change magnitudes within $ G_k $ are sufficiently large and well-separated.

FDR control via SD filter

Since the synthetic data in Eq. (8) mimic the distributional behavior under the null, the distributions of $ T_k^{q} $ and $ \tilde{T}_k^{q} $ are close when $ \hat{\tau}_k \in \mathcal{I}_0 $. Therefore, a comparison statistic is defined for each candidate change-point $ \hat{\tau}_k $ as follows:

discussed previously.

Moreover, the odd dataset $ \mathcal{Z}^O $, used to detect the candidate change-points, provides valuable side information. For each $ \hat{\tau}_k $, the CUSUM statistic can also be computed using $ \mathcal{Z}_{G_k}^O $, denoted by $ T_k^q(\mathcal{Z}_{G_k}^O) $. If the candidate set $ \hat{\mathcal{S}} $ is a plausible estimate of the true set $ \mathcal{S} $, it is typically observed that $ T_k^q(\mathcal{Z}_{G_k}^O) \gt T_{k^\prime}^q(\mathcal{Z}_{G_{k^\prime}}^O) $ for $ \hat{\tau}_k \in \mathcal{I}_1 $ and $ \hat{\tau}_{k^\prime} \in \mathcal{I}_0 $. By incorporating this information, a new statistic $ W_k^{side,q} $ is introduced. This blends the original comparison statistic with the side information:

After incorporating the odd dataset $ \mathcal{Z}^O $, the statistic still exhibits symmetry around zero in the case of a false discovery. The advantage of $ W_k^{side,q} $ lies in its ability to enhance the separation between the comparison statistics corresponding to informative and uninformative sets. Specifically, when $ \hat{\tau}_k \in \mathcal{I}_1 $ and $ \hat{\tau}_{k^\prime} \in \mathcal{I}_0 $, the ratio

since the second factor is typically greater than or equal to one.

Building on the definitions of $ W_k^q $ and $ W_k^{side,q} $, the candidate set is refined by selecting positions with large values of the comparison statistics. Specifically, it is defined as:

where, $ \mathcal{T}(t) $ refered to as the selection set. Given such a set, the number of false discoveries can be estimated by exploiting the symmetric of $ W_k^q $ or $ W_k^{side,q} $ around zero when $ \hat{\tau}_k \in \mathcal{I}_0 $. This leads to the following relationship:

where, $ W_k^{side,q} $ can be replaced by $ W_k^q $. Based on this property, the FDR can be approximated by

To control the FDR at a target level $ \alpha $, the knockoff+ procedure ^[22] is followed and a data-dependent threshold $ T(\alpha) $ is computed as follows:

where, $ \mathcal{W} = \{ W_1, \ldots W_{ \hat{K}_n}\} \backslash \{0\} $ and the extra term 1 in the numerator makes the choice of $ T(\alpha) $ slightly conservative.

In summary, compared with pure data-splitting methods such as MOPS and M-MOPS, the core methodological distinction of this approach lies in how the test statistics are constructed. Instead of performing inference on only half of the data, the method used generates a synthetic dataset via a Gaussian multiplier bootstrap, which integrates information from both the odd part and the reserved even part. This design increases the effective sample size and thus improves statistical power. Moreover, although the construction of $ W_{k}^{q} $ in Eq. (13) differs from that in MOPS and M-MOPS, it serves a conceptually similar role as the comparison statistic used in data-splitting–based procedures. As shown in Eq. (13), the formulation effectively amplifies the contrast between informative and uninformative subsets, which in turn facilitates better control of the false discovery rate.

The error rate control results rely on the symmetry property of the comparison statistics $ W_k^q $ and $ W_k^{side,q} $ when $ \hat{\tau}_k $ is a false discovery. To lay the groundwork for FDR control, it is essential to systematically examine this symmetry property. Before presenting the main theorem, some regular conditions are first imposed.

Condition 3.1 (Moments and tails). Let $ \underline{b}, \overline{b} $, and any vector $ \boldsymbol{\vartheta} \in \mathbb{S}^d $, and for all $ \ell = 1, 2 $, $ i = 1, \ldots, 2n $ and $ j = 1, \ldots, d $. (1)$ \text*{E}[(\boldsymbol{\vartheta}^\top \boldsymbol{\zeta}_{i})^2] \geqslant \underline{b} $; (2)$ \text*{E}[(\boldsymbol{\vartheta}^\top \boldsymbol{\zeta}_{i})^{\ell+2}] \leqslant \overline{b}^{\ell} $; (3) $ \| \boldsymbol{\vartheta}^\top \boldsymbol{\zeta}_{i}\|_{\psi_1} \leqslant \overline{b} $, where for $ \beta \in [1, \infty) $, $ \| \cdot \|_{\psi_\beta} $ represents an Orlicz norm.

Condition 3.2 (Detection ability). Assume $ \hat{K}_n \geqslant K_n $. There exist $ \hat{\tau}_{j_1} \lt \ldots \lt \hat{\tau}_{j_{K_n}} $ belonging to $ \mathcal{T} $ such that $ \max_{1 \lt k \lt K_n} | \hat{\tau}_{j_k} - \tau_k| \leqslant \delta_n $ holds with probability approaching one as $ n \to \infty $, where $ \delta_n $ is some positive sequence.

Condition 3.3 (Minimum distance). Assume that $ \mathcal{T} \subseteq \mathcal{T}(\omega_n) = \{ \mathcal{T}: \min_{j}(\tau_{j+1} - \tau_j) \geqslant \underline{\lambda}_n\} $, where $ \underline{\lambda}_n $ is a positive sequence such that $ \underline{\lambda}_n \geqslant n^{\eta} $ for some $ 0<\eta \lt 1 $ and $ \underline{\lambda}_n \geqslant 2\delta_n $.

A sufficient condition for Condition 3.1 (1) is that the minimum eigenvalue of $ \boldsymbol{\Sigma}^{(k)} $ is uniformly bounded below for all $ k = 1, \ldots, K_n $. This is a common assumption when the dimension $ d $ is fixed. Condition 3.1 (2) is a mild moment condition of the linear transformation of $ \boldsymbol{\zeta}_i $. Condition 3.1 (3) assumes a sub-exponential tail for the transformation of the residuals. Note that $ \| \boldsymbol{\vartheta}^\top \boldsymbol{\zeta}_i\|_{\psi_1} \leqslant \sum_{j =1}^d \|\vartheta_j \zeta_{ij}\|_{\psi_1} \leqslant \sum_{j =1}^d \| \zeta_{ij} \|_{\psi_1} $. Therefore, if $ \sum_{j =1}^d \| \zeta_{ij} \|_{\psi_1} \leqslant \bar{b} $, then Condition 3.1 (3) is satisfied. Conditions similar to Condition 3.1 can be found in Liu et al.^[1] and Yu & Chen^[13] for change-point inference. Condition 3.2 requires that the set of candidate change-points $ \hat{\mathcal{S}} $ is sufficiently accurate. Condition 3.3 imposes sparsity on the true change-points. Both Condition 3.2 and Condition 3.3 are also considered in Chen et al.^[30].

Let $ \overline{n} = \max_{k = 1,\ldots, K_n} n_k $ and $ \underline{n} = \min_{k = 1,\ldots, K_n} n_k $ denote the maximum and minimum sample sizes of the intervals G_k, respectively. The following lemma is now established.

Lemma 3.1. Assume Condition 3.1, 3.2 and 3.3 holds, then

where, $ \rho(T_1, T_2) = \sup_{t \in (0, \infty)} | \text{Pr}(T_1 \leqslant t) - \text{Pr}(T_2 \leqslant t) | $ denotes the Kolmogorov distance between $ T_1 $ and $ T_2 $. Here, $ \underline{s} $ is defined as $ \min_{k = 1,\ldots, K_n} \tau_{k+1} - \tau_{k} $, and $ \kappa $, $ C $, and $ c $ are positive constants.

Lemma 3.1 demonstrates that the distribution of $ T_k^q $ approximates that of its original CUSUM counterpart under mild conditions. Based on Lemma 3.1, the main results on FDR control are therefore presented.

Theorem 3.1. Under Condition 3.1, 3.2 and 3.3, and $ \log^7(\hat{K}_n \bar{n})/\underline{s} \to 0 $, the SD filter satisfies

for any $ \alpha \in (0,1) $.

Theorem 3.1 establishes the asymptotic FDR control property of the SD filter. The condition $ \log^7(\hat{K}_n \bar{n} d^2)/\underline{s} \to 0 $ implies that the bound parameter $ \underline{s} $ should not approach the endpoints too closely. This requirement is reasonable, since accurate covariance estimation requires sufficiently large samples on both sides of s. Unlike Theorem 1 in Chen et al.^[30], Theorem 3.1 does not require $ \delta_n/\underline{\lambda_n} \to 0 $, indicating that the SD filter does not rely on highly accurate candidate localization.

The FDR control for the general $ \ell_q $-norm cannot be extended to high-dimensional MCP settings due to the constraints of high-dimensional central limit theorem for simple and sparse convex sets, as discussed in Chernozhukov et al.^[32]. To address the high-dimensional MCP challenges, a practical approach is to set $ q = \infty $. In this case, the high-dimensional central limit theorem for hyperrectangles can be employed to justify Lemma 3.1 even in high-dimensional MCP scenarios.

Theorem 3.2 (FDR control for high-dimensional MCP). When $ d \to \infty $, under Condition 3.1, 3.2 and 3.3, and $ \log^7(\hat{K}_n \bar{n} d^2)/\underline{s} \to 0 $. Let $ q = \infty $, the SD filter satisfies

for any $ \alpha \in (0,1) $.

Power analysis

Next, the power of the SD filter is analyzed under the following signal condition.

Condition 3.4 (Minimum signal). Let $ \boldsymbol{\delta}_{k} = \boldsymbol{\mu}_{k+1} - \boldsymbol{\mu}_k $ be the change magnitude, which satisfies

where, $ t_k = \tau_k/n_k $, $ \alpha_n $ is a sequence that converges to infinity with a slow rate and $ C $ is a positive constant.

Condition 3.4 imposes a minimum signal separation between any two true change-points, ensuring their asymptotic identifiability. Similar conditions can be found in Harchaoui & Lévy-Leduc^[5], Fryzlewicz^[6], and Yu & Chen^[13].

Theorem 3.3. Under Condtion 3.1, 3.2, 3.3 and 3.4 and $ \log^7(\hat{K}_n \bar{n})/\underline{s} \to 0 $. The power of SD filter satisfies

Theorem 3.3 states that the power of SD filter approaches 1 asymptotically. Furthermore, the selection consistency property can be established.

Corollary 3.1. Under Conditions in Theorem 3.3, there is

Compared with the condition for selection consistency in Chen et al.^[30], which requires $ \min_{k \in \mathcal{I}_1}\| \boldsymbol{\delta}^{(k)}\|_2 \gg \sqrt{\log n/\underline{\lambda}_n} $, Condition 3.4 is of approximately the same order, noting that $ \underline{\lambda}_n \approx \underline{n} $.

In this section, aseries of change-point detection experiments is conducted to evaluate the empirical performance of the SD filter. Before presenting the results, the competing methods, the Mirror with Order-Preserved Splitting (MOPS) method and its variant, the Modified-MOPS (M-MOPS) are briefly summarized, both introduced in Chen et al.^[30].

The M-MOPS method controls the FDR via a mirror statistic

where, $ \overline{\mathbf{S}}_k^{O,-} $, $ \overline{\mathbf{S}}_k^{O,+} $, $ \overline{\mathbf{S}}_k^{E,-} $, $ \overline{\mathbf{S}}_k^{E,+} $ are previously defined, and $ \Omega_n $ is a positive matrix. Since the performance of M-MOPS is not sensitive to $ \Omega_n $, the study sets $ \Omega_n = \mathbf{I}_d $, the d-dimensional identity matrix, in the simulations. Compared with M-MOPS, the MOPS statistic differs only in that it uses a larger sample to compute the sample means of the score functions:

where, $ \tilde{\mathbf{S}}_k^{O,-} $ and $ \tilde{\mathbf{S}}_k^{O,+} $ denote the sample means of $ \{ \mathbf{s}_i^O, \hat{\tau}_{k-1} \lt i \leqslant \hat{\tau}_k\} $ and $ \{ \mathbf{s}_i^O, \hat{\tau}_{k} \lt i \leqslant \hat{\tau}_{k+1}\} $, respectively; $ \tilde{\mathbf{S}}_k^{E,-} $ and $ \tilde{\mathbf{S}}_k^{E,+} $ are defined analogously using the even subsample.

Then the computational complexity of the methods is compared. Treating basic arithmetic operations as O(1), computing $ W_{k}^{\mathrm{M}-\mathrm{MOPS}} $ requires O(dn) opreations, where d is the data dimension. Similarly, the quantities $ \mathbf{c}_{k}(s) $ and $ \tilde{\mathbf{c}}_{k}(s) $ also involve O(dn) operations, implying that both $ T_{k}^{q} $ and $ \tilde{T}_{k}^{q} $, and therefore $ W_{k}^{side,q} $, have the same complexity. This method generates only a single set of multiplier bootstrap samples and therefore adds no extra computational overhead. In practice, one may run $ B $ bootstrap repetitions, yielding a cost of O(Bdn); however, the procedure can be parallelized with ease, so the runtime can approach that of a single iteration. Overall, the SD filter does not incur high computational cost.

Beyond computational considerations, an important issue is statistical reliability. In particular, MOPS may fail to control the FDR when the discrepancy between the candidate and true change-point sets is large. To assess the performance of all methods, 200 simulation replications were conducted and each method was evaluated using the empirical FDR and power:

where, $ \mathcal{T}_i $ denotes the estimated selection set in the $ i $th replication.

The detailed pseudocode of this approach is as the Algorithm 1.

Simulation for multiple mean changes model

Consider a sequence of $ d $-dimensional mean vectors $ \boldsymbol{\mu}_{i}, i = 1,\ldots, n $, where the means are piecewise constant with change-points at positions $ \{\tau_k, k = 1, \ldots, K\} $. Specifically,

The sequence is initialized with $ \boldsymbol{\mu}_{\tau_1} = (A/2) \mathbf{1}_d $, where $ \mathbf{1}_d $ denotes the $ d $-dimensional vector of ones. To define $ \boldsymbol{\mu}_{\tau_2} $, randomly selecting $ r $ positions in $ \boldsymbol{\mu}_{\tau_1} $ and flip the sign of those entries. Subsequent vectors $ \boldsymbol{\mu}_{\tau_k} $ are constructed recursively from their predecessors using the same procedure. Under this setup, the study has

The data points are generated as $ \mathbf{z}_i = \boldsymbol{\mu}_i + \boldsymbol{\epsilon}_i, i = 1, \ldots, n $, and the signal strength is controlled by the value of $ A $. Throughout this subsection, the dimensionality is set to $ d = 50 $, the sample size to $ n=4,000 $. The true change-points are set to $ \tau_k=200\ k,\ k=1,\ldots,19 $, in Example 1, and to $ \tau_k=400\ k,\ k=1,\ldots,9 $, in Example 2. To mitigate the undesired bias from detection algorithms, a sequence of candidate change-points are manually constructed as $ \mathcal{T} = \{150 k+ (-1)^{\rm{B_k}} {\rm{P_k}} \mid k = 1, \ldots 26\} $, where $ \text{B}_{k} $, $ \text{P}_{k} $ are independently drawn from $ \text{Bernoulli}\ (1/2) $ and $ \text{Poisson }(5) $, respectively. Theoretically, $ q = \infty $ is required in high-dimensional settings, while any $ q \geqslant 1 $ is admissible in low-dimensional cases. Since the optimal choice of $ q $ depends on unknown signal conditions, $ q = \infty $ is adopted as a practical default in both simulations and empirical analysis. For the truncation size, $ s \in [10, 30] $ is chosen to ensure that both $ s $ and $ n_k - s $ are sufficiently large, thereby mitigating boundary effects and improving the stability of the algorithm. In addition, the study sets $ r = 1 $ and FDR level $ \alpha = 0.15 $.

Example 1: Normal distribution

Consider the error term $ \boldsymbol{\epsilon}_i $ to be drawn from multivariate normal distribution with covariance matrix $ \boldsymbol{\Sigma} = \{\rho^{|i-j|}\}_{(i,j)} $. The study chose the change magnitude $ A $ and the correlation coefficient $ \rho $ as follows:

● Fix $ \rho $ = 0, and let A vary in {1.5, 1.7, 1.9, 2.1, 2.3, 2.5}.

● Fix A = 1.5, and let $ \rho $ vary in {0, 0.2, 0.4, 0.6, 0.8}.

Example 2: Beyond normal distribution

Consider the error term $ \boldsymbol{\epsilon}_i $ to be drawn from either a multivariate t distribution or a multivariate chi-square distribution, each having a covariance matrix $ \boldsymbol{\Sigma} = \mathbf{I}_d $. We choose the change magnitude A and the degrees of freedom df as follows:

● Fix A = 2, and let df vary in {8, 9, 10, 11, 12}.

● Fix A = 3, and let df vary in {3, 4, 5, 6, 7}.

The simulation results summarized in Figs 3 and 4 demonstrate the superior performance of the SD filter across various signal strengths and dependence structures. The results indicate that both the SD filter and M-MOPS exhibit the capability to maintain control over the false discovery rate (FDR) at the predetermined level. In contrast, MOPS struggles to maintain FDR control due to the deliberate reduction in the quality of the candidate change-point set. When the candidate change positions are significantly distant from the actual change-points, MOPS fails to perform effectively. In terms of empirical power, the SD filter consistently outperforms the M-MOPS method across all settings.

Figure 3. 

FDR and power trends with respect to A and $ \rho $ for SD filter, MOPS, and M-MOPS under the mean change model, where n = 4,000, d = 50, $ \alpha $ = 0.15.

Figure 4. 

FDR and power trends with respect to df of the t-distribution and chi-square distribution for SD filter, MOPS, and M-MOPS under the mean change model, where n = 4,000, d = 50, $ \alpha $ = 0.15.

Structural breaks in linear regression model

Consider a linear regression model with structural breaks, defined as $ \mathbf{y}_i = \mathbf{x}_i^\top \boldsymbol{\beta}_{\tau_k} +\epsilon_i, \,\, \text{for}\,\, \tau_{k-1} \leqslant i \leqslant \hat{\tau}_k $. Initially, let $ \boldsymbol{\beta}_{\tau_1} = (A/2) \mathbf{1}_d $. Then, define $ \boldsymbol{\beta}_{\tau_2} $ by randomly selecting $ s $ positions in $ \boldsymbol{\beta}_{\tau_1} $ and flipping the signs of those entries. Each subsequent vector $ \boldsymbol{\beta}_{\tau_k} $ is generated from $ \boldsymbol{\beta}_{\tau_{k-1}} $ using the same procedure. Under this setup, the study has

The covariates $ \mathbf{x}_i $ are drawn from a multivariate normal distribution $ N(\mathbf{0}_{d}, \boldsymbol{\Sigma}) $, where $ \mathbf{0}_{d} $ is the d-dimensional vector of zeros and the covariance matrix $ \boldsymbol{\Sigma} = \{\rho^{|i-j|}\}_{(i,j)} $. The error terms $ \epsilon_i $ are independently drawn from $ N(0,1) $. In this model, the sample size is set as n = 8,000 and the number of covariates d = 10. The true change-point set is defined as $ \mathcal{S} = \{1000 k, k = 1, \ldots, 7\} $, and the candidate change-point set is given by $ \hat{\mathcal{S}} = \{ 450 k + (-1)^{\text{B}}_{\rm k} \text{P}_{\rm k} \mid k = 1, \ldots, 16\} $, where $ \text{B}_{k} \sim \text{Bernoulli}(1/2) $ and $ \text{P}_k\sim\text{Poisson }(5) $. Except that the FDR level is set to $ \alpha = 0.2 $, all other parameters remain unchanged. The change magnitude $ A $ and correlation coefficient $ \rho $ were chosen as follows:

● Fix $ \rho$ = 0 , and let A vary in {0.20, 0.22, 0.24, 0.26, 0.28, 0.30}.

● Fix A = 0.25, and let $ \rho $ vary in {0, 0.1, 0.2, 0.3, 0.4, 0.5}.

The simulation results presented in Figs 5 and 6 display the estimated FDR and empirical power results for the linear regression model with structural breaks under various scenarios. As before, it was observed that both the SD filter and M-MOPS successfully control the FDR at the pre-specified level. Again, MOPS fails to do so due to the deliberately reduced quality of the candidate change-point set. In terms of empirical power, the SD filter still consistently outperforms M-MOPS in this setting.

Figure 5. 

FDR and power trends with respect to A for SD filter, MOPS, and M-MOPS under the structural breaks linear regression model, where n = 8,000, d = 10, $ \alpha $ = 0.2.

Figure 6. 

FDR and power trends with respect to $ \rho $ for SD filter, MOPS and M-MOPS under the structural breaks linear regression model, where n = 8,000, d = 10, $ \alpha $ = 0.2.

In this section, the proposed SD filter is applied to analyze the bladder tumor micro-array dataset sourced from Bleakley & Vert^[35], which is conveniently available in the ecp R package. The dataset consists of log–intensity-ratio measurements for 2,215 genetic loci obtained from 43 individuals diagnosed with bladder tumors. The primary objective is to identify change-points within the genetic loci, enabling the study to pinpoint potential influential genes related to bladder tumors. This dataset has been widely used as a benchmark in several prior studies on change-point detection^[35–37], making it a representative and well-established dataset for evaluating the empirical performance of new methods. The analysis was conducted on the full dataset. However, for visualization purposes and to provide a clearer and more interpretable presentation of the results, the findings for the first ten individuals were reported, specifically individuals 3, 4, 5, 6, 7, 8, 9, 10, 14, and 15.

Firstly, the inspect method^[9] is applied to narrow the scope and obtain a candidate set. To ensure optimal performance, a minimum difference of $ 50 $ was established between two change-points. The set of identified change-points is as follows:

Subsequently, the SD filter, MOPS and M-MOPS are applied to further refine the results, controlling the FDR at a level of 0.1. The final sets of detected change-points from MOPS and M-MOPS are identical to $ \hat{\mathcal{S}} $, whereas the set of change-points detected by the SD filter is as follows:

This result of SD filter excludes the position 1,831. Figure 7 visually demonstrates the change-points identified through the SD filter. In each plot, individual data points represent log-intensity ratios on a specific genetic locus, with each plot corresponding to a different test subject. Vertical lines are used to indicate the locations of detected change-points. The change-points identified by the SD filter are shown as dashed lines. Notably, the only solid line—positioned at 1,831—does not correspond to any apparent change across the ten individuals, highlighting the greater precision and accuracy of the SD filter in identifying true change-points.

Figure 7. 

Detected change-points on bladder tumor micro-array dataset (first ten persons are presented).

To overcome the limitations of existing FDR control methods for multiple change-point detection-particularly the drawbacks associated with data-splitting approaches, the study proposes a synthetic data filter (SD filter) for change-point detection and FDR control. After identifying potential change-points, Gaussian multiplier bootstrap is applied to generate synthetic data based on information from the detection dataset. This synthetic data is then used to construct a mirror statistic that enables control of the FDR, offering the flexibility to leverage information from the entire dataset and improve statistical power under a variety of alternatives and dimensions. The symmetry property of the mirror statistic is then established andits ability to rigorously control FDR asymptotically is proven. The detection power is also demonstrated under mild conditions. Simulation studies empirically verified the outstanding performance of the SD filter in terms of FDR control and power. The study also applies the proposed method to analyze a micro array dataset that describes the change loci of bladder tumor patients. As mentioned above, the framework of the SD filter is quite general, and it would be interesting to extend it to a broader range of cases where data splitting is required for formulating and testing hypotheses within the same dataset.