Online Change Point Detection for Multivariate Inhomogeneous Poisson Processes Time Series

Technical Summary: Online Change Point Detection for Multivariate Inhomogeneous Poisson Processes Time Series

Problem Statement
This work addresses the challenge of online change point detection for multivariate inhomogeneous Poisson point process (PPP) time series. While PPPs are widely used to model random events in space with location-dependent rates (e.g., earthquakes, wildfires, epidemics), existing literature largely focuses on offline estimation or parametric models. A significant gap remains in developing nonparametric, online methods capable of handling multivariate PPPs where consecutive observations exhibit temporal dependence. Real-world data often displays such dependence due to environmental persistence or latent dynamics, yet standard change detection methods typically assume independence across time or rely on strong parametric assumptions that do not fit general intensity function changes.

Methodology
The authors propose a computationally efficient, single-pass online detection algorithm that maps infinite-dimensional intensity functions to finite-dimensional low-rank matrices. The core methodology proceeds as follows:

1. Matrix Representation via Basis Expansion:
   The intensity function $\lambda^*$ is represented using an orthonormal basis expansion (specifically within a Reproducing Kernel Hilbert Space framework). The domain $\mathcal{X} \subset \mathbb{R}^d$ is partitioned into two disjoint coordinate subsets of dimensions $p$ and $q$ (where $p+q=d$). The intensity function is approximated by a matrix $M(\lambda^*) \in \mathbb{R}^{M^p \times M^q}$, where entries are inner products of the intensity with tensor-product basis functions. This representation preserves distances and allows for efficient computation of $L^2$ norms without Monte Carlo integration.

2. CUSUM Statistic Construction:
   For each observed PPP sample $X^{(i)}$, an empirical intensity matrix $\hat{M}^{(i)}$ is constructed by summing the basis function evaluations over the observed points. The algorithm maintains running sums of these matrices in two lists, $L$ and $R$, corresponding to pre-change and post-change windows. At each time step $j$, a CUSUM statistic matrix $D$ is computed as the difference between the average of matrices in the pre-change window and the post-change window.

3. Restricted SVD and Trimming:
   To handle high dimensionality and reduce variance, the algorithm applies a Restricted Singular Value Decomposition (RSVD) to the difference matrix $D$.

   • Adaptive Trimming: The basis size is adaptively trimmed based on the effective sample size of the current window. Smaller windows only support the reliable estimation of a smaller number of matrix coefficients.
   • Low-Rank Projection: The trimmed matrix is projected onto its top $r$ singular vectors. This leverages the assumption that the difference between pre- and post-change intensity functions ($\lambda^* - \lambda^*_a$) is approximately low-rank, a property supported by functional PCA theory for smooth functions.

4. Detection Rule:
   An alarm is raised if the Frobenius norm of the rank-$r$ projected matrix, $\|D^{[r]}\|_F$, exceeds a data-driven threshold that scales with the window size and smoothness parameters.

Key Theoretical Contributions
The paper establishes finite-sample guarantees for the proposed method under a geometrically $\beta$-mixing temporal dependence assumption.

• New Matrix Bernstein Inequality: A primary theoretical contribution is the derivation of a new Matrix Bernstein inequality tailored for geometrically $\beta$-mixing Poisson point process time series. This inequality controls the deviation of the sum of dependent random matrices, which is essential for handling the temporal dependence in the PPP samples.
• False Alarm Control: The authors prove that with a sufficiently large threshold constant, the overall probability of a false alarm over the entire time horizon is bounded by $\alpha$.
• Detection Delay: Under a single change point scenario, the algorithm is shown to detect a change within a delay $\Delta$ that depends explicitly on the magnitude of the change in the $L^2$-norm ($\kappa$) and the smoothness of the intensity functions. Notably, the detection delay bound scales as $O(\kappa^{-2-(p \vee q)/\gamma})$, which is strictly better in order than existing offline nonparametric bounds ($O(\kappa^{-2-d/\gamma})$) because the method exploits the coordinate partition ($p \vee q < d$).

Empirical Results
The method, termed the "Matrix detector," was evaluated against two benchmarks: an MMD-based detector (Li et al., 2015) and a Kernel Intensity Estimator (KIE) based detector (Madrid Padilla et al., 2023).

• Simulation Studies: Experiments in 3D and 4D settings with temporally dependent latent variables demonstrated that the Matrix detector achieves significantly lower average detection delays (ADD) compared to the benchmarks while maintaining comparable false alarm probabilities.
• Robustness: Additional 10D simulations showed the method remains effective even when the intensity difference is not exactly low-rank, outperforming a neural-network-based CUSUM detector. The method was also shown to be robust to the choice of rank parameter $r$ and the coordinate partition strategy.
• Real Data Application: The algorithm was applied to earthquake data from Oklahoma (2000–2018). It successfully detected a change point in June 2009, aligning with known increases in seismic activity linked to wastewater injection, whereas the KIE detector raised an alarm too early (2007) and the MMD detector failed to raise an alarm by the end of the study period.

Significance and Claims
The authors claim that this work provides the first general-purpose approach with provable guarantees for detecting nonparametric changes in multivariate PPP intensities under realistic temporal dependence. The significance of the work lies in:

1. Handling Temporal Dependence: Moving beyond the independence assumption common in point process change detection by utilizing $\beta$-mixing conditions.
2. Computational Efficiency: Offering a single-pass algorithm with constant computational cost per new observation, independent of the elapsed time series length, making it scalable for long-term monitoring.
3. Theoretical Rigor: Providing non-asymptotic bounds that jointly account for basis-truncation bias and stochastic variance in dependent settings, supported by a novel Matrix Bernstein inequality.

The paper concludes that the proposed method is both statistically robust and computationally efficient, offering a practical solution for applications such as earthquake seismology, climate monitoring, and epidemic surveillance where event data is multivariate, inhomogeneous, and temporally dependent.