Fast Gaussian Process Approximations for Autocorrelated Data

AI summary

Overview

Research area: Machine learning / Gaussian process regression, specifically scalable approximation methods for data with temporal autocorrelation.

Technical level: Advanced. The paper assumes familiarity with Gaussian process regression, covariance kernels, Vecchia approximations, maximin ordering, and partial autocorrelation functions.

Scope: The paper adapts three existing fast Gaussian process approximations (local approximate GP, twin GP, and scaled Vecchia) to work with the data-thinning scheme used in the tempGP model so that autocorrelated data can be modeled quickly without temporal overfitting.

What This Paper Is About

Gaussian process regression is a popular nonlinear modeling tool, but inverting its covariance matrix costs O(n³), which is prohibitive for large datasets; fast approximations such as laGP, twinGP, and scaled Vecchia solve this by using nearest-neighbor subsets. Those approximations assume independently and identically distributed (i.i.d.) noise, so applying them directly to autocorrelated data reinforces autocorrelation and produces a failure mode called temporal overfitting, which harms prediction on new data. The goal of this paper is to modify each of the three approximations so that they operate on "thinned" (de-correlated) data blocks in the manner of the tempGP model, thereby recovering speed without sacrificing accuracy.

Key Contributions

1. Thinned twinGP. A scheme that trains a twinGP on each of the T thinned data blocks, averages the block-wise predictions as the point estimate, and combines the block predictive standard deviations into a single predictive standard deviation for the averaged prediction (Equation 7). The authors also report that the three internal size parameters of twinGP (global-set, local-set, and hold-out validation-set sizes) must be taken from a twinGP fit on the original dataset rather than re-selected automatically on the smaller thinned blocks; smaller automatically selected parameters led to poor performance.

2. Thinned laGP. Because laGP tunes its own hyperparameters per test point and is slow, the authors use a single-block approach: find the training point with smallest Euclidean distance to the test point, locate the thinned block containing it, and train/predict with only that block. This is faster than training laGP on every block, but since T ranges from 6 to 30 in their applications, only 3–15% of the data is used for training and prediction.

3. Thinned scaled Vecchia (SV). A modified conditioning-set construction: after the tempGP-style thinning, maximin ordering (using the algorithm of Schäfer et al. 2021) is applied within each block, and each point's m nearest neighbors are chosen only from points that precede it in that block's ordering. For prediction, a separate conditioning set of size m_p is built by maximin ordering over the full (unthinned) training and test data, with test locations predicted in order and appended to the ordered training set as predictions are made. Generally m_p > m.

4. Comparative evaluation on four datasets. The paper compares the three thinned methods against tempGP and against two unthinned SV variants — SV(x), which ignores time, and SV(x, t), which cascades x and t — on the robot arm simulation, the satellite drag coefficient data, and DSWE Datasets 5 and 6, using RMSE and NLPD.

Main Findings

• Thinned SV has the best overall accuracy on autocorrelated data. The introduction states that "thinned SV has the overall best model performance on autocorrelated data."
• Thinned twinGP is the fastest approximation. Among the thinned methods, twinGP is reported as fastest.
• Thinned laGP does not appear competitive relative to the other two thinned alternatives, which the authors attribute in part to the single-block design that uses only 3–15% of the data.
• Speedup motivation is substantial. Fitting tempGP on 50,000 data points takes more than 3 hours on high-performance computers, whereas a fast GP method may take just a few minutes on the same size.
• Including time as an input does not reliably fix the problem. SV(x, t) helps under some circumstances but not always; on the robot arm and DSWE Dataset 6 test sets, the time indices are far from those in training, so SV(x, t) is expected to extrapolate poorly in the time dimension.
• The g(·) term is small in variance. When the predictive variance for the thinned methods is approximated by adding var(f(·)) and var(g(·)), the magnitude of var(g(·)) is about 1% of var(f(·)); ignoring the covariance between f and g does not fundamentally change the overall picture. This approximation applies only to the satellite drag coefficient data and DSWE Dataset 5.
• Autocorrelation strength is measured by T. For the robot arm scenarios, T falls in the range of 6 to 26 for M = 13 and M = 24; across the paper's applications T ranges from 6 to 30.
• A residual caveat on metric coverage. tempGP is included only in the RMSE comparison and not in NLPD, because it does not produce predictive variances as model output.

Methodology in Plain English

The tempGP model splits the response into two additive parts: a time-invariant function f(x) of the inputs, and a time-indexed function g(t) that absorbs the autocorrelation, so that the leftover noise can be treated as i.i.d. To estimate f, tempGP "thins" the time-ordered data into T interleaved blocks constructed so that points within a block are essentially uncorrelated — the thinning number T is chosen as the smallest lag h at which the partial autocorrelation function of every covariate and the response drops below 2/√n, which is roughly the 95% confidence threshold under no autocorrelation. The paper then asks how to plug each of the three fast GP approximations into this thinned block structure.

For twinGP, the answer is to fit one model per block and average. For laGP, which is too slow to fit T times, the answer is to pick the single block nearest in input space to the test point. For scaled Vecchia, the answer is to build the conditioning sets entirely within blocks, so that the neighbors chosen are spatially similar but temporally far apart. After f is estimated by whichever fast method, the residuals y − f̂(x) are used to estimate g(t) using the standard approach from the original tempGP paper, with a Matérn covariance of smoothness ν = 1.5 and only the temporal neighborhood [t − T, t + T]. Because T is small — usually in the tens and definitely under hundreds — estimating g does not need a fast approximation. The final prediction is f̂(x*) + ĝ(t*).

Experiments use four datasets: the robot arm simulation (with Latin hypercube sampling plus two autocorrelation-injected scenarios, using autoregressive structure on joint angles and segment lengths), the satellite drag coefficient data, and DSWE Datasets 5 and 6. Randomized cross-validation is used for the satellite drag and DSWE Dataset 5; the robot arm uses separate simulated training and test sets; DSWE Dataset 6 trains on one year of historical data and tests on a different year. Performance is measured by RMSE (point accuracy) and NLPD (posterior distribution alignment); smaller is better. Parameter settings were m = 30 and m_p = 140 for the SV models (matching Katzfuss et al. 2022), and laGP's neighborhood size was set to 30. All experiments ran on five Intel CPU cores on the RHL7 Phoenix cluster at Georgia Tech's PACE facility.

Why This Matters

Impact on research. The work bridges two largely separate lines of methodological development — scalable GP approximation and autocorrelation-aware modeling — by showing that the two can be combined through a modest restructuring of neighborhood selection. It also supplies concrete evidence that naively adding a time index to a scalable GP's inputs (SV(x, t)) is not a reliable substitute for proper thinning.

Real-world applications (drawn from the paper's own motivating examples):

• Wind power production forecasting — the authors use this as their lead example, since wind speed is autocorrelated through air-mass inertia and drives power output, while unmodeled factors leave residual autocorrelation in the response.
• Wind energy analytics — two of the four datasets (DSWE Dataset 5 and Dataset 6) are the companion datasets for the book Data Science for Wind Energy.
• Robotic arm trajectory modeling — joint angles and segment lengths cannot change instantaneously, so smooth temporal transitions are physically required.
• Aerospace / satellite drag coefficient estimation — the satellite drag coefficient dataset is one of the four benchmarks.

More broadly, the paper targets data generated by physical and engineering systems with automated collection (inline sensors producing evenly spaced sequences), where the i.i.d. assumption typically fails.

Industry relevance. Operators of large sensor-instrumented systems — wind farms, manufacturing equipment, robotics, vehicles — face exactly this trade-off: models that are accurate enough to run on tens of thousands of points within minutes, while still respecting the temporal structure that makes naive validation optimistic. Moving a 50,000-point fit from more than 3 hours to minutes is operationally significant for iterative model development and deployment.

Future Directions

1. Extending the variance approximation. The additive var(f) + var(g) simplification was only needed for the two datasets where g is non-zero; a principled treatment of the covariance between f and g remains open.
2. Improving thinned laGP. The single-block design discards 85–97% of the data when T is large; a cheaper way to aggregate across blocks without retraining laGP T times would address its uncompetitive accuracy.
3. Better use of time information within the fast approximations. SV(x, t) failed in cases of extreme temporal extrapolation, raising the question of how to incorporate time-aware structure without extrapolating along the time axis.
4. Stress-testing across autocorrelation regimes. The paper varies autocorrelation strength through M = 13 and M = 24 in simulation; broader coverage of autocorrelation levels, block sizes, and dataset scales would clarify when each thinned method is preferable.

Target Audience

Researchers and practitioners working on scalable Gaussian process methods, spatio-temporal statistics, and surrogate modeling for expensive simulations will get the most value, along with applied statisticians and machine learning engineers handling large autocorrelated sensor or time-series datasets where i.i.d. assumptions are untenable. Readers who already know the tempGP model, Vecchia approximations, or the twinGP/laGP literature will find the methodological modifications immediately actionable; readers without GP background should expect a steep entry.

Authors’ abstract

This paper is concerned with the problem of how to speed up computation for Gaussian process models trained on autocorrelated data. The Gaussian process model is a powerful tool commonly used in nonlinear regression applications. Standard regression modeling assumes random samples and an independently, identically distributed noise. Various fast approximations that speed up Gaussian process regression work under this standard setting. But for autocorrelated data, failing to account for autocorrelation leads to a phenomenon known as temporal overfitting that deteriorates model performance on new test instances. To handle autocorrelated data, existing fast Gaussian process approximations have to be modified; one such approach is to segment the originally correlated data points into blocks in which the blocked data are de-correlated. This work explains how to make some of the existing Gaussian process approximations work with blocked data. Numerical experiments across diverse application datasets demonstrate that the proposed approaches can remarkably accelerate computation for Gaussian process regression on autocorrelated data without compromising model prediction performance.