RFF-GP-HSMM: Scalable Time-Series Segmentation

• The paper introduces RFF-GP-HSMM, which replaces computationally heavy kernel inversion in GP-HSMM with a Random Fourier Feature approximation for scalable segmentation.
• It transforms GP regression into linear regression on M-dimensional features, significantly reducing training cost while preserving segmentation quality.
• Empirical results on human motion data demonstrate up to 278× speedup with nearly identical segmentation accuracy compared to traditional GP-HSMM.

RFF-GP-HSMM is a scalable, unsupervised time-series segmentation method that replaces the Gaussian-process emissions in a Gaussian Process Hidden Semi-Markov Model with a Random Fourier Feature approximation. Its purpose is the same as GP-HSMM: to segment a continuous sequence into interpretable, recurring patterns with explicit durations and without labels. The method addresses the central computational limitation of GP-HSMM, namely the need to invert an $N \times N$ kernel matrix during training for each state’s Gaussian process, which yields an $O(N^3)$ bottleneck and $O(N^2)$ memory use as the number of data points grows. By approximating a stationary kernel with a finite-dimensional feature map, RFF-GP-HSMM reduces GP regression to linear regression or Bayesian linear regression on $M$-dimensional features with $M \ll N$, thereby preserving expressive power while eliminating the need for kernel-matrix inversion (Saito et al., 14 Jul 2025).

1. Definition and objective

RFF-GP-HSMM was introduced in "Scalable Unsupervised Segmentation via Random Fourier Feature-based Gaussian Process" (Saito et al., 14 Jul 2025). It is formulated as a fast unsupervised time-series segmentation method designed for settings in which standard GP-HSMM becomes computationally prohibitive. The target problem is segmentation of sequential observations into latent states that recur over time, with each segment having an explicit duration distribution rather than the geometric-duration assumption of a standard HMM.

The motivating contrast is between expressiveness and scalability. GP-HSMM models state-conditioned emissions with Gaussian processes, which are flexible enough to capture structured temporal patterns, but this flexibility requires inversion of a kernel matrix whose size grows with the number of data points assigned to a state. RFF-GP-HSMM retains GP-like emission structure through a stationary-kernel approximation while moving inference into a finite-dimensional linear model. This suggests a design goal of maintaining segmentation quality while substantially reducing both training and inference cost.

The paper’s empirical setting is human motion segmentation. In that application, the input variable is time $t$, the observations are 8-D vectors corresponding to hands and feet coordinates, and conditional independence across output dimensions is assumed given the state and input (Saito et al., 14 Jul 2025).

2. Probabilistic structure

The underlying sequence model is an HSMM with explicit segment durations. Let $s_1, s_2, \ldots$ denote the sequence of states, let $d_1, d_2, \ldots$ denote segment durations, and let $y_{1:T}$ denote the observations. The generative structure is

$$p(s_1) \prod_j p(d_j \mid s_j) p(s_{j+1} \mid s_j) \prod_j p(y_{\tau_j} : \tau_j + d_j - 1 \mid s_j),$$

where $O(N^3)$0 is the start time of segment $O(N^3)$1.

State transitions follow a categorical distribution with a Dirichlet prior over rows. The transition probability used is

$O(N^3)$2

where $O(N^3)$3 is the number of states, $O(N^3)$4 counts transitions from $O(N^3)$5 to $O(N^3)$6, and $O(N^3)$7 is the symmetric Dirichlet hyperparameter. Segment durations are explicitly modeled by a Poisson distribution,

$O(N^3)$8

with inference restricted to $O(N^3)$9.

In the original GP-HSMM emission model, each state $O(N^2)$0 has a Gaussian-process regression function

$O(N^2)$1

and observations are generated as

$O(N^2)$2

Given state-specific training inputs $O(N^2)$3 and outputs $O(N^2)$4, the GP predictive distribution for a new input $O(N^2)$5 is

$O(N^2)$6

where $O(N^2)$7 has entries $O(N^2)$8, $O(N^2)$9, and $M$0. The dominant cost is computation of $M$1, which is $M$2 per state (Saito et al., 14 Jul 2025).

3. Random Fourier Feature approximation

The distinctive step in RFF-GP-HSMM is replacement of the GP emission by a Random Fourier Feature approximation of a stationary kernel. For any shift-invariant kernel $M$3 with Fourier transform $M$4, the method samples frequencies and phases as

$M$5

and defines the feature map

$M$6

so that

$M$7

For the RBF kernel used in the paper’s exposition,

$M$8

the corresponding spectral density is Gaussian,

$M$9

that is, $M \ll N$0. The paper derives the cosine-with-random-phase construction by writing the inverse Fourier integral, introducing a random phase $M \ll N$1, and using the identity $M \ll N$2, followed by Monte Carlo estimation of the integral (Saito et al., 14 Jul 2025).

Once $M \ll N$3 is introduced, the GP emission becomes a linear-Gaussian model. For each output dimension,

$M \ll N$4

with Gaussian prior

$M \ll N$5

typically with $M \ll N$6 and $M \ll N$7. For state $M \ll N$8, the posterior covariance and mean are

$M \ll N$9

$t$0

and the predictive distribution is

$t$1

with

$t$2

An equivalent ridge-regression view is

$t$3

where $t$4. The computational consequence is that the model no longer requires inversion of an $t$5 kernel matrix, replacing it with an $t$6 system where $t$7 (Saito et al., 14 Jul 2025).

4. Inference procedure and computational profile

The learning procedure uses blocked Gibbs sampling together with forward filtering-backward sampling over segments and state labels, and Bayesian linear regression updates for the emissions. Initialization consists of random segmentation and random segment-class assignment, followed by initialization of transition counts and class-specific linear-Gaussian emission parameters.

For each sequence $t$8, the forward recursion computes messages over segment ending time $t$9, duration $s_1, s_2, \ldots$0, and class $s_1, s_2, \ldots$1:

$s_1, s_2, \ldots$2

For $s_1, s_2, \ldots$3-dimensional observations, the emission likelihood is factorized as

$s_1, s_2, \ldots$4

Backward sampling then draws segment durations and state labels from the end of the sequence using the paper’s form

$s_1, s_2, \ldots$5

assigns the segment $s_1, s_2, \ldots$6, and repeats with $s_1, s_2, \ldots$7 until the start of the sequence is reached. After segments are reassigned, the method updates state-specific sufficient statistics $s_1, s_2, \ldots$8 and $s_1, s_2, \ldots$9, recomputes $d_1, d_2, \ldots$0 and $d_1, d_2, \ldots$1 or the corresponding ridge estimate, and updates transition probabilities from the Dirichlet-multinomial counts. In the reported experiments, the Poisson mean $d_1, d_2, \ldots$2 and the truncation bounds $d_1, d_2, \ldots$3 are fixed rather than learned. The blocked Gibbs outer loop repeats sequence removal, parameter update from the remaining data, and FFBS resampling for all sequences; the paper uses 5 iterations (Saito et al., 14 Jul 2025).

The complexity reduction is explicit. For GP-HSMM, training per state is $d_1, d_2, \ldots$4 with $d_1, d_2, \ldots$5 memory, and prediction is $d_1, d_2, \ldots$6 per test point. For RFF-GP-HSMM, sufficient-statistics accumulation is $d_1, d_2, \ldots$7 if $d_1, d_2, \ldots$8 is formed directly, or $d_1, d_2, \ldots$9 with streaming updates and small $y_{1:T}$0; posterior or MAP computation is $y_{1:T}$1; emission evaluation is $y_{1:T}$2 per time point, multiplied by the observation dimensionality $y_{1:T}$3. The HSMM forward-backward dynamic program costs $y_{1:T}$4 for transitions, and with RFF emissions the additional emission cost is $y_{1:T}$5. The summarized complexity is therefore

• GP-HSMM: dominated by $y_{1:T}$6 with $y_{1:T}$7 memory.
• RFF-GP-HSMM: $y_{1:T}$8 training, $y_{1:T}$9 inference, and $$p(s_1) \prod_j p(d_j \mid s_j) p(s_{j+1} \mid s_j) \prod_j p(y_{\tau_j} : \tau_j + d_j - 1 \mid s_j),$$0 memory.

A plausible implication is that the principal scalability gain comes not only from replacing cubic-statewise GP training but also from removing statewise predictive costs that would otherwise recur inside FFBS.

5. Hyperparameters and empirical evaluation

The reported experiments use an RBF kernel with $$p(s_1) \prod_j p(d_j \mid s_j) p(s_{j+1} \mid s_j) \prod_j p(y_{\tau_j} : \tau_j + d_j - 1 \mid s_j),$$1 and $$p(s_1) \prod_j p(d_j \mid s_j) p(s_{j+1} \mid s_j) \prod_j p(y_{\tau_j} : \tau_j + d_j - 1 \mid s_j),$$2, together with $$p(s_1) \prod_j p(d_j \mid s_j) p(s_{j+1} \mid s_j) \prod_j p(y_{\tau_j} : \tau_j + d_j - 1 \mid s_j),$$3 random features. Bayesian linear regression uses the prior $$p(s_1) \prod_j p(d_j \mid s_j) p(s_{j+1} \mid s_j) \prod_j p(y_{\tau_j} : \tau_j + d_j - 1 \mid s_j),$$4 and observation noise parameter $$p(s_1) \prod_j p(d_j \mid s_j) p(s_{j+1} \mid s_j) \prod_j p(y_{\tau_j} : \tau_j + d_j - 1 \mid s_j),$$5, although the numerical constants are not specified. The duration model is $$p(s_1) \prod_j p(d_j \mid s_j) p(s_{j+1} \mid s_j) \prod_j p(y_{\tau_j} : \tau_j + d_j - 1 \mid s_j),$$6 with $$p(s_1) \prod_j p(d_j \mid s_j) p(s_{j+1} \mid s_j) \prod_j p(y_{\tau_j} : \tau_j + d_j - 1 \mid s_j),$$7, truncated in FFBS with $$p(s_1) \prod_j p(d_j \mid s_j) p(s_{j+1} \mid s_j) \prod_j p(y_{\tau_j} : \tau_j + d_j - 1 \mid s_j),$$8 and $$p(s_1) \prod_j p(d_j \mid s_j) p(s_{j+1} \mid s_j) \prod_j p(y_{\tau_j} : \tau_j + d_j - 1 \mid s_j),$$9. The number of states $O(N^3)$00 is fixed but not explicitly reported. Initialization is random, and each experiment uses 5 blocked-Gibbs iterations. The implementation extends a GP-HSMM implementation with fast FFBS attributed to Sasaki et al., 2023; a code release is listed at https://github.com/naka-lab/RFF-GP-HSMM (Saito et al., 14 Jul 2025).

The dataset is CMU motion capture, Subject 14, using three sequences containing multiple exercises: jumping jack, jogging, squatting, knee raise, arm circle, twist, side reach, boxing, arm wave, side bend, and toe touch. The data are downsampled to 4 fps, represented as 8-D inputs consisting of 2D positions of left and right hands and feet, and min-max normalized to $O(N^3)$01.

For segmentation accuracy, the three sequences are duplicated 10 times to produce 30 sequences. Each method is run 10 times with different initializations, and the best-likelihood run is reported. The metric is normalized Hamming distance, where 0 is best. The reported values are $O(N^3)$02 for GP-HSMM and $O(N^3)$03 for RFF-GP-HSMM, which the paper characterizes as essentially unchanged segmentation accuracy under the RFF approximation (Saito et al., 14 Jul 2025).

For speed evaluation, the same three sequences totaling 490 frames are duplicated to create datasets with 3, 30, 60, and up to 240 sequences, corresponding to 490, 4,900, 9,800, 14,700, 19,600, 24,500, 29,400, 34,300, and 39,200 frames. Each method is run 5 times and average wall-clock time is reported on an Apple M1 Ultra with 64 GB RAM.

Frames	Time comparison	Speedup
9,800	GP-HSMM 708 s; RFF-GP-HSMM 21 s	34×
19,600	3,227 s vs 48 s	67×
34,300	21,294 s vs 97 s	220×
39,200	32,583 s vs 117 s	278×

At 39,200 frames, which the paper describes as approximately 2 h 43 min of data at 4 fps, the runtime drops from about 9 h to about 2 min. Across all evaluated sizes, RFF-GP-HSMM consistently outperforms GP-HSMM, and the speedup increases with dataset size. The paper also compares this result to prior GP-side optimizations, noting that approximately $O(N^3)$04 acceleration at approximately 8.5k frames had been reported in Sasaki et al., 2023, whereas RFF-GP-HSMM achieves approximately $O(N^3)$05 at 9.8k frames (Saito et al., 14 Jul 2025).

6. Relation to adjacent methods, limits, and usage considerations

The paper positions RFF-GP-HSMM among several families of alternatives. Inducing-point or sparse GP methods and Nyström methods reduce the effective rank of the kernel matrix and replace $O(N^3)$06 with $O(N^3)$07 for $O(N^3)$08, but they still retain GP-specific machinery, require per-state inducing sets, and often involve nontrivial optimization of inducing locations. By contrast, RFF uses explicit features and reduces GP regression to linear regression, yielding constant-time-in-$O(N^3)$09 matrix inversion at the level of an $O(N^3)$10 system and straightforward integration with linear-Gaussian updates and HSMM dynamic programs. The stated trade-off is approximation error controlled by $O(N^3)$11 and random sampling variance (Saito et al., 14 Jul 2025).

Relative to segmentation baselines, HMMs or HSMMs with Gaussian emissions are described as fast but limited in emission expressiveness for complex motion. GP-HSMM addresses this limitation with nonparametric emissions but incurs cubic cost in the number of state-assigned data points. RFF-GP-HSMM is therefore presented as restoring GP-like emission flexibility to the HSMM framework at near-HMM runtimes. Transformer-based and deep sequence models are described as capable of capturing long-range structure but as typically requiring labels, large data, and heavy compute, while offering less interpretability than HSMMs.

The paper does not report explicit ablations varying the number of random features, kernel families, or duration models. Instead, it fixes $O(N^3)$12, the RBF kernel, and Poisson durations with $O(N^3)$13, $O(N^3)$14, and $O(N^3)$15. It nevertheless provides guidance on expected behavior rather than reported results. Increasing $O(N^3)$16 is said to improve the kernel approximation and typically improve accuracy up to diminishing returns while increasing per-state $O(N^3)$17 to $O(N^3)$18 cost. Very small $O(N^3)$19 can underfit, while very large $O(N^3)$20 erodes the speed advantage. For periodic or quasi-periodic motion, the text suggests that periodic or mixed stationary kernels could help, provided their spectral densities are used to sample $O(N^3)$21. It also notes that more flexible duration families such as log-normal or negative binomial could better match real segment lengths but would complicate inference. Truncation bounds $O(N^3)$22 and $O(N^3)$23 strongly affect runtime because the dynamic program scales linearly in $O(N^3)$24.

The practical guidance given in the paper recommends starting with $O(N^3)$25 in $O(N^3)$26, increasing it if accuracy lags GP-HSMM or a nonparametric baseline, and considering $O(N^3)$27 to $O(N^3)$28 for high-dimensional inputs or very long sequences. With one-dimensional time inputs, smaller values around $O(N^3)$29 to $O(N^3)$30 are said often to suffice. For the number of states $O(N^3)$31, the recommendation is to set a plausible upper bound based on domain knowledge; overestimation is described as often safer because the Dirichlet prior regularizes unused transitions, although per-time-step computation grows linearly or quadratically in $O(N^3)$32 depending on whether emissions or transitions dominate. Tight duration bounds are recommended when approximate segment lengths are known. A plausible implication is that RFF-GP-HSMM is most suitable when a stationary kernel is appropriate, segment durations are semantically meaningful, and the application benefits from interpretable unsupervised segmentation together with strong computational scaling (Saito et al., 14 Jul 2025).