RuLSIF

Change-point detection by relative density-ratio estimation.

Description

The series is cut into overlapping subsequences of length $k$ (subsequence_length); a multivariate subsequence concatenates all channels. For every candidate split, the $n$ subsequences before it (distribution $P$) are compared with the $n$ subsequences after it (distribution $P'$, n_subsequences) by the symmetrised $\alpha$-relative Pearson divergence

\[\mathrm{PE}_\alpha(P \| P') + \mathrm{PE}_\alpha(P' \| P), \qquad \mathrm{PE}_\alpha(P \| P') = \frac{1}{2} \int p'_\alpha(Y) \left(r_\alpha(Y) - 1\right)^2 \mathrm{d}Y,\]

with $p'_\alpha = \alpha p + (1 - \alpha) p'$ and the relative density ratio $r_\alpha = p / p'_\alpha$. Each divergence is estimated without estimating densities: RuLSIF fits $r_\alpha$ by regularised least squares with a Gaussian kernel model centred on the samples of $P$, which has a closed-form solution. As in the authors’ code, the subsequences of the two sets are first standardised coordinate by coordinate, so that channels on different scales weigh the same; the kernel width (a multiple of the median distance between subsequences) and the regularisation parameter are chosen by 5-fold cross-validation for every estimate. The relative ratio is bounded by $1 / \alpha$, which makes the estimate more stable than that of the plain ratio ($\alpha = 0$, uLSIF).

Change points are the peaks of this score at least min_distance points apart: the n_cps highest ones in guided mode, those above threshold otherwise. The paper only evaluates the score over all thresholds (ROC curves); the default threshold (2.25, for $\alpha = 0.1$) is calibrated by tsseg: it maximises the mean F1 score (margin 1 %) over the paper’s four synthetic datasets. The score saturates near its bound $(1 - \alpha) / \alpha$ on smooth or drifting series (motion capture, random walk), where no threshold separates the changes: use n_cps there. The score is kept in the scores_ attribute after predict.

Type: change point detection
Supervision: unsupervised (threshold) or semi-supervised (n_cps)
Scope: univariate and multivariate
Complexity: \(O(T \cdot F S L \cdot n^3)\) for a series of length
\(T\), \(F\) folds, \(S\) kernel widths and \(L\)
regularisation parameters; about 40 s for 5,000 points with the defaults
Requires: nothing beyond the core dependencies

Parameters

Name	Type	Default	Description
`subsequence_length`	int	`10`	Subsequence length $k$ (paper value).
`n_subsequences`	int	`50`	Subsequences per compared set, $n$ (paper value). The series needs at least $2n + k - 1$ points.
`alpha`	float	`0.1`	Mixture weight of the relative density ratio, in $[0, 1)$ (paper value).
`sigma_factors`	tuple of float	`(0.6, 0.8, 1.0, 1.2, 1.4)`	Candidate kernel widths, as multiples of $\sqrt{\operatorname{median}(d^2) / 2}$ over the pairwise distances $d$ between subsequences: the paper’s factors, the authors’ code’s base width (the paper’s median distance over $\sqrt 2$).
`lambdas`	tuple of float	`(1e-3, 1e-2, 1e-1, 1, 10)`	Candidate regularisation parameters (paper values).
`n_folds`	int	`5`	Cross-validation folds (paper value).
`n_cps`	int / None	`None`	Number of change points; `None` thresholds the peaks.
`threshold`	float	`2.25`	Minimum peak score without `n_cps`, calibrated for $\alpha = 0.1$ (tsseg). The score lies in $[0, (1 - \alpha) / \alpha]$: scale the threshold with this bound for another $\alpha$.
`min_distance`	int	`20`	Minimum distance between change points; the higher of two closer peaks is kept (the paper drops alarms closer than 20 points to the previous one).
`axis`	int	`0`	Time axis.

Usage

from tsseg.algorithms import RuLSIFDetector

# peaks above the threshold
cps = RuLSIFDetector().fit_predict(X)

# known number of change points
detector = RuLSIFDetector(n_cps=3)
cps = detector.fit_predict(X)
score = detector.scores_

Implementation: NumPy port of the authors’ MATLAB code (anewgithubname/change_detection, MIT licence), checked against a line-by-line transcription of it; it reproduces the code’s demo on well-log data (reproduction test).

Reference: Liu, Yamada, Collier & Sugiyama (2013), Change-point detection in time-series data by relative density-ratio estimation, Neural Networks 43; Yamada, Suzuki, Kanamori, Hachiya & Sugiyama (2013), Relative density-ratio estimation for robust distribution comparison, Neural Computation 25(5).

API reference

RuLSIF change-point detector (Liu, Yamada, Collier and Sugiyama, 2013).

class tsseg.algorithms.rulsif.detector.RuLSIFDetector(*, subsequence_length=10, n_subsequences=50, alpha=0.1, sigma_factors=(0.6, 0.8, 1.0, 1.2, 1.4), lambdas=(0.001, 0.01, 0.1, 1.0, 10.0), n_folds=5, n_cps=None, threshold=2.25, min_distance=20, axis=0)[source]

Bases: BaseSegmenter

Change-point detection by relative density-ratio estimation (RuLSIF).

Each split of the series is scored by how much the distribution of the n_subsequences subsequences before it differs from that of the n_subsequences subsequences after it. The difference is the symmetrised $\alpha$-relative Pearson divergence $\mathrm{PE}_\alpha(P \| P') + \mathrm{PE}_\alpha(P' \| P)$, each term estimated directly, without density estimation, by RuLSIF (a regularised least-squares fit of the relative density ratio with a Gaussian kernel model). Change points are peaks of this score.

Parameters:

subsequence_length (int) – Length k of the subsequences that form the samples; a multivariate subsequence concatenates all channels. Paper value.
n_subsequences (int) – Number n of subsequences in each of the two compared sets. A score needs 2 * n_subsequences + subsequence_length - 1 points, which is also the minimum series length. Paper value.
alpha (float) – Mixture weight $\alpha \in [0, 1)$ of the relative density ratio $p / (\alpha p + (1 - \alpha) p')$; 0 gives the plain density ratio (uLSIF). Paper value.
sigma_factors (tuple[float, ...]) – Candidate Gaussian kernel widths, as multiples of the median scale $\sqrt{\operatorname{median}(d^2) / 2}$ of the pairwise distances $d$ between the subsequences of the two sets (after standardisation, see Notes). The paper’s factors; the base width is the authors’ code’s (the paper’s median distance divided by $\sqrt 2$).
lambdas (tuple[float, ...]) – Candidate regularisation parameters. Paper values.
n_folds (int) – Number of cross-validation folds used to choose the kernel width and the regularisation parameter, separately for every divergence estimate. Paper value.
n_cps (int | None) – Number of change points to return (guided mode): the n_cps highest peaks of the score at least min_distance apart, or fewer if the score has fewer peaks. None keeps the peaks above threshold.
threshold (float) – Minimum score of a peak in unsupervised mode. Not from the paper, which evaluates the score over all thresholds (ROC curves). Calibrated by tsseg for the default parameters (alpha=0.1): it maximises the mean F1 score (margin 1 %) over the paper’s four synthetic datasets (Section 4.1, three draws each), 0.700 against 0.678 for 3.0; on change-free series of 2,000 points it gives 2 to 4 false alarms (white noise, noisy sine, correlated Gaussian, 5-channel noise), 18 on an AR(2) process. The score is bounded by (1 - alpha) / alpha (9 for the default alpha; its estimate can fall slightly below 0), so for another alpha scale the threshold with this bound; a threshold at or above it finds nothing and raises a warning.
min_distance (int) – Minimum number of points between two change points: of two peaks closer than that, the higher one is kept. The paper, when it evaluates the score, drops an alarm closer than 20 points to the previous one (Section 4.1).
axis (int) – Time axis.

scores_

Change-point score of the last series passed to predict, nan where it is undefined (about n_subsequences + subsequence_length / 2 points at each end).

Type:: np.ndarray of shape (n_timepoints,)

Notes

Port of the authors’ MATLAB code (change_detection.m, RelULSIF.m, demo.m; https://github.com/anewgithubname/change_detection, MIT licence), checked against a line-by-line NumPy transcription of it (equal divergence estimates to rounding). Where the code is more specific than the paper, the code is followed:

Before each estimate, every coordinate of the 2 * n_subsequences subsequences is divided by its standard deviation over them (not centred), so channels on different scales weigh the same; a nearly constant channel is brought to unit variance too, which amplifies its noise.
The base kernel width is $\sqrt{\operatorname{median}(d^2) / 2}$ over the non-zero pairwise distances (comp_med.m), i.e. the paper’s “median distance” divided by $\sqrt 2$.
The centres of the kernel model are the samples of the numerator set. Cross-validation draws, for every (kernel width, regularisation) candidate, a random partition of each set into n_folds folds; the code resets its generator before every estimate, so the partitions are the same for all estimates. MATLAB’s generator cannot be reproduced: NumPy’s is used with a fixed seed, which keeps the detector deterministic. The code subsamples 100 centres when n_subsequences > 100; tsseg keeps them all.
The symmetric score is the sum of the code’s score on the series and on the reversed series (demo.m), i.e. both directions of the divergence on the same pair of sets.

tsseg’s own choice: the score of the sets starting at t and t + n_subsequences is placed at t + n_subsequences + (subsequence_length - 1) // 2, where the subsequences that straddle the boundary between the two sets are split evenly between them (demo.m plots it at the end of the second set, 2 * n_subsequences + subsequence_length - 2 points after t, which suits online monitoring).

The cost is $O(T \cdot F \cdot S \cdot L \cdot n^3)$ for a series of length $T$, $F$ folds, $S$ kernel widths and $L$ regularisation parameters: about 8 ms per point with the defaults, i.e. 40 s for 5,000 points.

References

Examples

>>> import numpy as np
>>> from tsseg.algorithms import RuLSIFDetector
>>> rng = np.random.default_rng(0)
>>> X = np.concatenate([rng.normal(0, 1, 300), rng.normal(0, 4, 300)])
>>> RuLSIFDetector(n_cps=1).fit_predict(X)
array([300])

set_fit_request(*, axis: bool | None | str = '$UNCHANGED$') → RuLSIFDetector

Configure whether metadata should be requested to be passed to the fit method.

Note that this method is only relevant when this estimator is used as a sub-estimator within a meta-estimator and metadata routing is enabled with enable_metadata_routing=True (see sklearn.set_config()). Please check the User Guide on how the routing mechanism works.

The options for each parameter are:

True: metadata is requested, and passed to fit if provided. The request is ignored if metadata is not provided.
False: metadata is not requested and the meta-estimator will not pass it to fit.
None: metadata is not requested, and the meta-estimator will raise an error if the user provides it.
str: metadata should be passed to the meta-estimator with this given alias instead of the original name.

The default (sklearn.utils.metadata_routing.UNCHANGED) retains the existing request. This allows you to change the request for some parameters and not others.

Added in version 1.3.

Parameters:: axis (str, True, False, or None, default=sklearn.utils.metadata_routing.UNCHANGED) – Metadata routing for axis parameter in fit.
Returns:: self – The updated object.
Return type:: object

set_predict_request(*, axis: bool | None | str = '$UNCHANGED$') → RuLSIFDetector

Configure whether metadata should be requested to be passed to the predict method.

Note that this method is only relevant when this estimator is used as a sub-estimator within a meta-estimator and metadata routing is enabled with enable_metadata_routing=True (see sklearn.set_config()). Please check the User Guide on how the routing mechanism works.

The options for each parameter are:

True: metadata is requested, and passed to predict if provided. The request is ignored if metadata is not provided.
False: metadata is not requested and the meta-estimator will not pass it to predict.
None: metadata is not requested, and the meta-estimator will raise an error if the user provides it.
str: metadata should be passed to the meta-estimator with this given alias instead of the original name.

The default (sklearn.utils.metadata_routing.UNCHANGED) retains the existing request. This allows you to change the request for some parameters and not others.

Added in version 1.3.

Parameters:: axis (str, True, False, or None, default=sklearn.utils.metadata_routing.UNCHANGED) – Metadata routing for axis parameter in predict.
Returns:: self – The updated object.
Return type:: object