Pith. sign in

REVIEW 4 major objections 5 minor 34 references

Adaptive Nystr\"om for Gaussian Process Regression

T0 review · 4 major / 5 minor · reviewed 2026-07-31 · grok-4.5

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read Greedy trace-residual landmarks, updated with the kernel itself, let Gaussian process regression match exact accuracy at linear cost in sample size.

desk verdict Solid engineering paper: greedy residual Nyström interleaved with hyperparameter fits beats random on small smooth benchmarks, but the distinctive claim lacks the obvious static-greedy control. read the letter →

arxiv 2607.27427 v1 pith:C33QXHYJ submitted 2026-07-29 stat.ME stat.ML

Lulu Kang This is my paper · ORCID
classification stat.MEstat.ML MSC 62M2062G0865F3068T05
keywords GaussianprocessregressionNyströmapproximationlandmarkselectiontraceresidualscalableGPRcomputerexperimentsWoodburyidentityhyperparameteroptimization
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

Gaussian process regression gives reliable predictions and uncertainty, but the usual matrix work grows as the cube of the number of data points, so it stalls on large computer experiments. A standard fix is the Nyström low-rank shortcut that rebuilds the kernel from a small set of landmark points; the catch is that random landmarks often miss the structure that matters. This paper’s claim is that you can choose those landmarks greedily by repeatedly picking the training point with the largest leftover kernel residual, and that you should re-estimate the kernel hyperparameters after each addition so the residual is measured under the current covariance. On five standard benchmark simulators the adaptive loop reaches accuracy close to a full Gaussian process while using only a small fraction of the points as landmarks, and it is more stable than random selection. The practical payoff is a linear-in-n procedure that still behaves like exact GP inference for the kinds of smooth design problems common in computer experiments.

What carries the argument

The adaptive Nyström loop: at each step add the training point that maximizes the residual variance r(x_i) = K(x_i,x_i) − k_Z(x_i)^T W^{-1} k_Z(x_i) (equivalently the pivoted-Cholesky pivot), then re-minimize the Woodbury-reduced Nyström profiled likelihood until both the hyperparameter change and the mean residual fall below fixed tolerances.

What would settle it

On a larger or less smooth computer-experiment suite, run the same adaptive loop against exact GP and random Nyström with matched landmark count: if adaptive S-RMSPE stays far above exact GP, or fails to beat random selection with lower variance across random starts, the central claim fails.

Watch

Extended reading notes

Core claim

Interleaving greedy maximization of the kernel trace residual with hyperparameter re-optimization produces a Nyström landmark set that yields predictive accuracy comparable to exact Gaussian process regression and clearly better and more stable than uniform random landmarks, while the dominant cost stays linear in the training size.

Load-bearing premise

The method assumes that stopping when hyperparameter change and average residual are small, with landmarks forced to be a subset of the observed design, is enough to keep the low-rank surrogate faithful for near-exact predictions in general—not only on the smooth Latin-hypercube benchmarks tested.

Editorial extensions

If this is right

  • For smooth computer experiments, near-exact GP predictions can be obtained with landmark counts on the order of a few percent of n once the residual criterion and interleaved updates are used.
  • Uniform random landmark selection is an unreliable baseline for Nyström GPR accuracy and stability under the same computational budget.
  • Restricting landmarks to observed design points still leaves enough flexibility for the residual-greedy rule to cover the kernel feature space on the reported problems.
  • The same residual-driven expansion can be stopped by simple relative tolerances, giving a practical accuracy–cost dial without continuous inducing-point optimization.
  • When n is only a few hundred, sequential re-optimization can erase the asymptotic gain; the method’s advantage appears as n grows.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • The equivalence to pivoted Cholesky under a moving kernel suggests the method could inherit pivoted-Cholesky error bounds if those bounds were re-proved with slowly changing length-scales.
  • Because landmarks stay discrete and inside the training set, the approach may be easier to warm-start or parallelize than variational inducing-point methods that move continuous coordinates.
  • A natural stress test left open by the paper is multimodal or discontinuous response surfaces, where greedy residual peaks might cluster and leave large unexplored regions.
  • Automated schedules for the residual and hyperparameter tolerances could remove the remaining hand-tuned knobs the authors flag for future work.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

4 major / 5 minor

Summary. The paper proposes an adaptive Nyström approximation for Gaussian process regression in which landmark points are chosen greedily to maximize the kernel trace residual, and landmark expansion is interleaved with re-optimization of the profiled Nyström likelihood (via Woodbury and the matrix-determinant lemma). Landmarks are restricted to the observed design, and stopping is controlled by relative change in θ and mean residual. On five VLSE benchmark functions (n≤1000, Matérn-1.5, constant and polynomial means), the method is compared to uniform random Nyström of matched final m and to exact GP, and is reported to reach near-exact predictive accuracy with greater stability than random selection while scaling as O(T E n m²).

Significance. If the interleaving procedure reliably yields near-exact GP accuracy at linear cost in n, the work would be a useful, implementable contribution for large-scale computer experiments, where practitioners often prefer discrete landmarks drawn from the design and derivative-free optimizers. The algorithmic reductions in §3.1 are standard and correctly stated, the selection rule is cleanly tied to residual variance / pivoted Cholesky (§3.5), and the empirical protocol (unified custom R code, matched-m random baseline, S-RMSPE on held-out LHDs, 10 seeds) is transparent. The distinctive claim, however, is the value of interleaving relative to static greedy selection; that claim is only partially supported by the present evidence, so the incremental significance over known pivoted-Cholesky / residual-greedy Nyström is currently modest.

major comments (4)
  1. [§3.5, Algorithm 1; §4] §3.5 identifies interleaving of landmark growth with hyperparameter updates (Algorithm 1) as the primary novelty relative to static pivoted Cholesky on K. Section 4 compares only to uniform random Nyström of matched final m and to full GP. There is no static control that (i) builds the same greedy/pivoted-Cholesky landmark set under fixed or once-estimated θ and (ii) optimizes the profiled likelihood (2) a single time. Without that ablation, gains on the five smooth LHD problems cannot be attributed to interleaving rather than greedy coverage alone, so the central methodological claim is under-supported.
  2. [Abstract; §3.4; Table 2; §5] The abstract and §3.4–§5 claim linear scaling and a framework for large-scale experiments where exact O(n³) inference is infeasible. All reported runs use n∈{500,1000}. In Table 2 (OTL, polynomial mean) adaptive Nyström is slower than full GP (30.66s vs 23.28s), which the text correctly attributes to sequential overhead at modest n. No experiment with n large enough that full GP is impractical is provided, so the scaling claim remains an asymptotic argument rather than a demonstrated advantage.
  3. [Abstract; Table 3; §4.2] Table 3 (Wing Weight, polynomial mean) shows random Nyström with lower S-RMSPE than adaptive (0.0112 vs 0.0169). The text notes possible high-d clustering of greedy landmarks but does not quantify when residual-greedy selection fails or how defaults δθ=10^{-3}, δr=10^{-4} and m_max interact with that failure mode. Given that the abstract asserts the method “significantly outperforms random landmark selection in both accuracy and stability,” this counterexample needs either a clearer scope restriction or additional diagnostics (e.g., residual maps, length-scale trajectories) before the claim is left unqualified.
  4. [§3.3; §4 experimental design] §3.3 defaults (δθ=10^{-3}, δr=10^{-4}, landmarks ⊂ X, modest m0/m_max) are presented as consistently balancing accuracy and cost, yet the weakest modeling assumption is that the resulting low-rank surrogate remains faithful when the residual surface is multimodal or length-scales shift sharply. The five smooth VLSE functions with maximin LHDs do not stress this regime. At minimum, a sensitivity study over (δθ, δr) and a brief comparison to free inducing-point / variational sparse GP baselines (already cited in §3.5) would better bound the claim.
minor comments (5)
  1. [§4.1; Figure 1] Figure 1 caption and §4.1 refer to “Eq. (17)” for the residual, but the residual selection criterion is equation (3) in §3.2; numbering appears inconsistent.
  2. [§4] Abstract and introduction say “five benchmark functions,” but §4.1–4.2 detail Borehole plus four others (Piston, OTL, Wing Weight, Steel Column)—fine, yet Borehole results are only in Figure 1 and not tabulated like Tables 1–4; a short table for Borehole would aid comparison.
  3. [§2] Typo in §2: “conditional variane” → “variance.”
  4. [References] References include future-dated items (e.g., R Core Team 2026, lhs 2026, arXiv date 2026); verify bibliographic years for consistency with the venue.
  5. [§3.3] Predictive variance formula after Algorithm 1 drops the usual mean-correction term present in the exact GP variance in §2; state explicitly whether this is intentional under the Nyström approximation or an omission.

Circularity Check

0 steps flagged · score 0.0 of 10

No circularity: methods paper with held-out predictive benchmarks; selection criterion and claims are not self-defining.

full rationale

The paper proposes a greedy trace-residual landmark rule interleaved with hyperparameter optimization (Algorithm 1, Eqs. 2–3) and evaluates it via S-RMSPE on independent maximin-LHD test sets against full GP and random Nyström baselines (§4, Tables 1–4, Figure 1). Predictive metrics are not algebraic restatements of the training residual or profiled Nyström likelihood. Woodbury/determinant reductions and the residual formula are standard identities, not self-definitions of the claimed accuracy. §3.5 explicitly notes equivalence of the selection rule to pivoted Cholesky and positions interleaving as the novelty rather than renaming a known result as a new theorem. Citations justifying the trace criterion (e.g., KL bounds) are external. There is no fitted constant presented as a first-principles prediction, no load-bearing self-citation uniqueness claim, and no derivation that reduces by construction to its inputs. Ordinary hyperparameter fitting on training data with held-out comparison is not circularity under these criteria.

Assumptions & free parameters 5 free parameters · 5 assumptions · 0 invented entities

The work rests on classical GP and kernel-approximation machinery plus a few implementation choices. No new physical entities. Load-bearing modeling choices are stationary anisotropic kernels, landmarks subset of training data, and user tolerances that stop landmark growth. Empirical claims further depend on LHD synthetic benchmarks and BOBYQA optimization behavior.

free parameters (5)
  • initial landmark count m0 = 20–40 depending on experiment
    Hand-chosen starting size for the greedy set (e.g., 20–40); affects path of adaptive selection and runtime.
  • maximum landmarks m_max = 100–200
    Hard cap on landmark growth; chosen per experiment (100–200).
  • hyperparameter change tolerance δθ = 1e-3
    Stopping threshold on relative change in θ; default fixed by authors, directly controls early exit.
  • mean trace residual tolerance δr = 1e-4
    Stopping threshold on average residual variance; primary accuracy–cost knob.
  • kernel hyperparameters θ, noise ratio η, variance τ², mean β = data-dependent MLE / profile estimates
    Standard GP parameters optimized from the Nyström profiled likelihood each outer iteration; not physical constants.
assumptions (5)
  • domain assumption Responses follow a Gaussian process with stationary positive-definite kernel (Gaussian or Matérn) plus optional nugget η.
    Section 2 model (1) and covariance definitions; standard GP regression assumption.
  • standard math Woodbury identity and matrix-determinant lemma give exact likelihood reductions for the rank-m Nyström-plus-nugget covariance.
    Section 3.1; classical linear algebra used to obtain O(nm²) likelihood.
  • domain assumption Trace residual tr(K−K̃) upper-bounds (via cited sparse variational theory) the KL gap relevant to approximation quality, justifying greedy residual maximization.
    Section 3.2 citing Titsias / Wild et al.; selection criterion inherits that bound’s relevance.
  • ad hoc to paper Inducing/landmark points may be restricted to a subset of the observed training inputs without needing continuous variational optimization over free inducing locations.
    Section 3.5 design choice contrasting with variational sparse GPs; simplifies optimization but can limit approximation power.
  • domain assumption BOBYQA derivative-free optimization with ≤2000 evaluations reliably finds adequate hyperparameters for the profiled Nyström likelihood.
    Section 4 experimental protocol; unproven but common in kriging practice.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Adaptive Nystr\"om for Gaussian Process Regression." pith.science (2026). https://pith.science/paper/C33QXHYJ

@misc{pith2026260727427,
  author       = {Pith},
  title        = {Pith review of: Adaptive Nystr\"om for Gaussian Process Regression},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/C33QXHYJ}},
  note         = {Machine review of arXiv:2607.27427}
}
abstract

Gaussian Process Regression (GPR) is a robust framework for uncertainty quantification, yet its $O(n^3)$ complexity limits its scalability. Low-rank Nystr\"om approximations can reduce this burden to $O(nm^2)$, but their accuracy depends heavily on the selection of landmark points. We propose an adaptive Nystr\"om approach that greedily selects landmarks to minimize the trace residual of the kernel approximation error. Unlike static approximations, our method interleaves landmark expansion with hyperparameter optimization, allowing the selection process to adapt as the covariance structure is refined. Numerical experiments on five benchmark functions demonstrate that this method significantly outperforms random landmark selection in both accuracy and stability. It achieves predictive performance comparable to exact GP inference while maintaining linear scaling with respect to the sample size, providing a principled and efficient framework for large-scale computer experiments.

Figures

Figures reproduced from arXiv: 2607.27427 by the authors.

Figure 1
Figure 1. Borehole Example: S-RMSPE of the three methods with respect to increasing number of landmark [PITH_FULL_IMAGE:figures/full_fig_p008_1.png] view at source ↗

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

34 extracted references · 1 linked inside Pith

  1. [1]

    Kernel independent component analysis.Journal of Machine Learning Research, 3(Jul):1–48, 2002

    Francis R Bach and Michael I Jordan. Kernel independent component analysis.Journal of Machine Learning Research, 3(Jul):1–48, 2002

    2002

  2. [2]

    R package version 1.2.1

    Rob Carnell.lhs: Latin Hypercube Samples, 2026. R package version 1.2.1

    2026

  3. [3]

    Deep mixed effect model using gaussian processes: A personalized and reliable prediction for healthcare

    Ingyo Chung, Saehoon Kim, Juho Lee, Kwang Joon Kim, Sung Ju Hwang, and Eunho Yang. Deep mixed effect model using gaussian processes: A personalized and reliable prediction for healthcare. InProceedings of the AAAI conference on artificial intelligence, volume 34, pages 3649–3657, New York Hilton Midtown, New York, NY , USA, February 7–12 2020

    2020

  4. [4]

    John Wiley and Sons, New York, 2015

    Noel Cressie.Statistics for Spatial Data. John Wiley and Sons, New York, 2015

    2015

  5. [5]

    Distributed gaussian processes

    Marc Deisenroth and Jun Wei Ng. Distributed gaussian processes. In Francis Bach and David Blei, editors,Proceedings of the 32nd International Conference on Machine Learning, volume 37 ofPro- ceedings of Machine Learning Research, pages 1481–1490, Lille, France, 07–09 Jul 2015. PMLR

    2015

  6. [6]

    Gaussian processes for data-efficient learning in robotics and control.IEEE Transactions on Pattern Analysis and Machine Intelligence, 37(2):408–423, 2015

    Marc Peter Deisenroth, Dieter Fox, and Carl Edward Rasmussen. Gaussian processes for data-efficient learning in robotics and control.IEEE Transactions on Pattern Analysis and Machine Intelligence, 37(2):408–423, 2015

    2015

  7. [7]

    Forecasting global climate drivers using gaussian processes and convolutional autoencoders.Engineering Applications of Artificial Intelligence, 128:107536, 2024

    James Donnelly, Alireza Daneshkhah, and Soroush Abolfathi. Forecasting global climate drivers using gaussian processes and convolutional autoencoders.Engineering Applications of Artificial Intelligence, 128:107536, 2024

    2024

  8. [8]

    Petros Drineas and Michael W. Mahoney. On the nystrom method for approximating a gram matrix for improved kernel-based learning.Journal of Machine Learning Research, 6(72):2153–2175, 2005

    2005

Show all 34 references
  1. [9]

    Spectral grouping using the nys- trom method.IEEE Transactions on Pattern Analysis and Machine Intelligence, 26(2):214–225, 2004

    Charless Fowlkes, Serge Belongie, Fan Chung, and Jitendra Malik. Spectral grouping using the nys- trom method.IEEE Transactions on Pattern Analysis and Machine Intelligence, 26(2):214–225, 2004

    2004

  2. [10]

    Harville.Matrix Algebra From a Statistician’s Perspective

    D.A. Harville.Matrix Algebra From a Statistician’s Perspective. Springer, 2008

    2008

  3. [11]

    Higham.Accuracy and Stability of Numerical Algorithms

    Nicholas J. Higham.Accuracy and Stability of Numerical Algorithms. Society for Industrial and Applied Mathematics, second edition, 2002

    2002

  4. [12]

    Johnson.The NLopt Nonlinear-Optimization Package, 2008

    Steven G. Johnson.The NLopt Nonlinear-Optimization Package, 2008

    2008

  5. [13]

    Parallel inference for massive distributed spatial data using low-rank models.Statistics and Computing, 27(2):363–375, Mar 2017

    Matthias Katzfuss and Dorit Hammerling. Parallel inference for massive distributed spatial data using low-rank models.Statistics and Computing, 27(2):363–375, Mar 2017. 11

    2017

  6. [14]

    Ensemble nystrom method

    Sanjiv Kumar, Mehryar Mohri, and Ameet Talwalkar. Ensemble nystrom method. In Y . Bengio, D. Schuurmans, J. Lafferty, C. Williams, and A. Culotta, editors,Advances in Neural Information Processing Systems, volume 22. Curran Associates, Inc., 2009

    2009

  7. [15]

    Sampling methods for the nyström method.The Journal of Machine Learning Research, 13(1):981–1006, 2012

    Sanjiv Kumar, Mehryar Mohri, and Ameet Talwalkar. Sampling methods for the nyström method.The Journal of Machine Learning Research, 13(1):981–1006, 2012

    2012

  8. [16]

    Making large-scale nyström approximation possible

    Mu Li, James T Kwok, and Bao-Liang Lu. Making large-scale nyström approximation possible. In Proceedings of the 27th International Conference on International Conference on Machine Learning, pages 631–638, 2010

    2010

  9. [17]

    Chien-Yu Peng and C. F. Jeff Wu. On the choice of nugget in kriging modeling for deterministic computer experiments.Journal of Computational and Graphical Statistics, 23(1):151–168, 2014

    2014

  10. [18]

    The bobyqa algorithm for bound constrained optimization without derivatives

    Michael JD Powell. The bobyqa algorithm for bound constrained optimization without derivatives. Cambridge NA Report NA2009/06, University of Cambridge, Cambridge, 26(26-46):1, 2009

    2009

  11. [19]

    R Foundation for Statistical Computing, Vienna, Austria, 2026

    R Core Team.R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, 2026

    2026

  12. [20]

    Carl Edward Rasmussen and Christopher K. I. Williams.Gaussian Processes for Machine Learning. The MIT Press, 11 2005

    2005

  13. [21]

    DiceKriging, DiceOptim: Two R packages for the analysis of computer experiments by kriging-based metamodeling and optimization.Journal of Statistical Software, 51(1):1–55, 2012

    Olivier Roustant, David Ginsbourger, and Yves Deville. DiceKriging, DiceOptim: Two R packages for the analysis of computer experiments by kriging-based metamodeling and optimization.Journal of Statistical Software, 51(1):1–55, 2012

    2012

  14. [22]

    Springer, 2003

    Thomas J Santner, Brian J Williams, William I Notz, and Brain J Williams.The Design and Analysis of Computer Experiments, volume 1. Springer, 2003

    2003

  15. [23]

    Sparse variational inference for generalized gp mod- els

    Rishit Sheth, Yuyang Wang, and Roni Khardon. Sparse variational inference for generalized gp mod- els. In Francis Bach and David Blei, editors,Proceedings of the 32nd International Conference on Machine Learning, volume 37 ofProceedings of Machine Learning Research, pages 1302...

    2015

  16. [24]

    Global versus local methods in nonlinear dimensionality reduc- tion

    Vin Silva and Joshua Tenenbaum. Global versus local methods in nonlinear dimensionality reduc- tion. In S. Becker, S. Thrun, and K. Obermayer, editors,Advances in Neural Information Processing Systems, volume 15. MIT Press, 2002

    2002

  17. [25]

    Smola and Bernhard Schökopf

    Alex J. Smola and Bernhard Schökopf. Sparse greedy matrix approximation for machine learning. InProceedings of the Seventeenth International Conference on Machine Learning, ICML ’00, page 911–918, San Francisco, CA, USA, 2000. Morgan Kaufmann Publishers Inc

    2000

  18. [26]

    Local and global sparse gaussian process approximations

    Edward Snelson and Zoubin Ghahramani. Local and global sparse gaussian process approximations. In Marina Meila and Xiaotong Shen, editors,Proceedings of the Eleventh International Conference on Artificial Intelligence and Statistics, volume 2 ofProceedings of Machine Learning ...

    2007

  19. [27]

    Genton and

    Yan Song, Wenlin Dai, and Marc G. Genton and. Large-scale low-rank gaussian process prediction with support points.Journal of the American Statistical Association, 0(0):1–12, 2024

    2024

  20. [28]

    A review of nyström methods for large-scale machine learning

    Shiliang Sun, Jing Zhao, and Jiang Zhu. A review of nyström methods for large-scale machine learning. Information Fusion, 26:36–48, 2015. 12

    2015

  21. [29]

    Surjanovic and D

    S. Surjanovic and D. Bingham. Virtual library of simulation experiments: Test functions and datasets. Retrieved June 8, 2026, fromhttp://www.sfu.ca/~ssurjano, 2013

    2026

  22. [30]

    Variational learning of inducing variables in sparse gaussian processes

    Michalis Titsias. Variational learning of inducing variables in sparse gaussian processes. In David van Dyk and Max Welling, editors,Proceedings of the Twelfth International Conference on Artificial Intelligence and Statistics, volume 5 ofProceedings of Machine Learning Resear...

    2009

  23. [31]

    Wild, Motonobu Kanagawa, and Dino Sejdinovic

    Veit D. Wild, Motonobu Kanagawa, and Dino Sejdinovic. Connections and Equivalences between the Nyström Method and Sparse Variational Gaussian Processes.ArXiv e-prints:2106.01121, 2021

    2021 arXiv

  24. [32]

    Using the nyström method to speed up kernel ma- chines

    Christopher KI Williams and Matthias Seeger. Using the nyström method to speed up kernel ma- chines. InProceedings of the 14th International Conference on Neural Information Processing Sys- tems, NIPS’00, page 661–667, Cambridge, MA, USA, 2000. MIT Press

    2000

  25. [33]

    Woodbury.Inverting Modified Matrices

    M.A. Woodbury.Inverting Modified Matrices. Memorandum Report / Statistical Research Group, Princeton. Department of Statistics, Princeton University, 1950

    1950

  26. [34]

    Tsang, and James T

    Kai Zhang, Ivor W. Tsang, and James T. Kwok. Improved nyström low-rank approximation and error analysis. InProceedings of the 25th International Conference on Machine Learning, ICML ’08, page 1232–1239, New York, NY , USA, 2008. Association for Computing Machinery. 13

    2008

Pith tools

Reviewed July 31, 2026 · model on record in the stance chip above.