{"id":"7ffaf4ce-7e30-48d3-a3cd-c4c3db7543c3","arxiv_id":"2507.06839","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":4.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":6,"one_line_summary":"The thesis shows that iterative linear solvers plus pathwise conditioning scale Gaussian processes to millions of data points, introducing SGD-based, dual-descent, warm-started, and latent-Kronecker methods for inference and hyperparameter learning.","lead":"This PhD thesis combines iterative solvers with pathwise conditioning to make Gaussian process models work on datasets with millions of points. It introduces stochastic gradient and dual descent methods, warm-start tricks, and latent Kronecker structure, reporting large speedups and competitive accuracy in regression, Bayesian optimization, and drug-binding prediction.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Central 'exact GP' claim rests on a fixed m=2000 random-feature prior with no sensitivity analysis; if this approximation is inaccurate, posterior samples and learned hyperparameters are biased.","rationale":"The reader's weakest assumption already identifies the random-feature approximation as the most fragile link, and my read agrees: the central claim that iterative solvers perform inference in the exact GP model depends on prior samples being faithful draws from GP(0, K), yet the dissertation fixes m=2000 random features without a sensitivity study. This concern is load-bearing because it affects posterior samples, predictive uncertainties, and the marginal-likelihood gradients used for hyperparameter optimisation. The missing proof for warm-starting bias (Proposition 5.2, deferred to an external appendix) is a genuine completeness issue, but it is less central than the random-feature approximation because Table 5.1 already gives strong empirical evidence that warm starting does not degrade predictive performance, whereas no analogous evidence is provided that m=2000 is adequate across all kernels and dataset scales. The inducing-point substitution f_X^[Z] -> f_X in Section 3.2.3 is also asserted rather than demonstrated, compounding the concern. A concrete sensitivity check on the number of random features would settle whether the approximation error is material. Since the reader already conditioned acceptance on reproducibility and the missing proof, my stress test does not alter the verdict; it strengthens the rationale for requiring the random-feature sensitivity analysis to be included.","tokens_in":56712,"tokens_out":5901,"duration_ms":64981,"concrete_test":"Run a sensitivity analysis on a representative large dataset (e.g., HOUSEELECTRIC or one of the five-million-point datasets) using the pathwise estimator with m = 500, 2000, 8000, and 32000 random features, keeping all other hyperparameters and solver settings fixed. Compare test RMSE, test NLL, and the final marginal-likelihood value. If the m=2000 results differ from the m=32000 results by more than a small margin (e.g., >5% NLL or >2% RMSE), the fixed random-feature budget is load-bearing and the exact-GP claim needs qualification. Additionally, on a small dataset, compare the inducing-point sampling objective with exact f_X^[Z] (drawn via Cholesky on the inducing points) against the f_X approximation to quantify the error in Section 3.2.3.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The dissertation positions iterative methods as performing 'approximate inference in the exact Gaussian process model up to a specified numerical tolerance' (Section 2.2.4), and the abstract claims scalable inference 'in the exact GP model.' However, every posterior sample and the pathwise marginal-likelihood gradient estimator use prior samples drawn from GP(0, Phi Phi^T) with m=2000 random Fourier features (Sections 2.2.2, 3.2.2, 5.2). Unless Phi Phi^T equals K exactly, the covariance of a pathwise sample is not the exact GP posterior covariance, so the 'exact model' claim is false whenever the random-feature approximation is imperfect. The fixed budget m=2000 is used for all datasets and kernels with no sensitivity analysis; for slow-spectral-decay kernels or high-dimensional inputs, 2000 features can leave substantial approximation error. This biases posterior samples, predictive uncertainties, and the gradient estimates used for hyperparameter learning in Chapter 5. In the inducing-point variant (Section 3.2.3), the sampling objective further replaces f_X^[Z] with f_X, justified only by a non-quantified assertion that the error is small 'when m is large and inducing points are close.' If these approximations do not hold at the scales tested, the state-of-the-art empirical claims may still stand, but the conceptual guarantee of inference in the exact GP model is unsupported.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The dissertation develops scalable Gaussian process inference by combining iterative linear-system solvers with pathwise conditioning. Chapter 3 introduces stochastic gradient descent for posterior means and samples, including an inducing-point extension and a spectral analysis of its implicit bias. Chapter 4 proposes stochastic dual descent, based on a dual objective with more favorable conditioning, random-coordinate gradient estimators, momentum, and geometric averaging. Chapter 5 contributes a pathwise estimator of marginal-likelihood gradients, warm starting of inner-loop solvers, and an analysis of early stopping, reporting speed-ups of up to 72×. Chapter 6 introduces latent Kronecker structure for kernel matrices and demonstrates scalability to datasets with up to five million examples. The main theoretical derivations in Chapters 3–5 are internally consistent, while the empirical claims are scoped to large-scale or ill-conditioned regression, Bayesian optimization, and molecular binding-affinity prediction.","tokens_in":57144,"tokens_out":5996,"duration_ms":68494,"significance":"If the technical caveats discussed below are resolved, this is a valuable consolidation and extension of iterative methods for Gaussian processes. The dual-objective equivalence in Proposition 4.1, the initial-distance comparison in Eqs. (5.12)–(5.13), the variance bound in Proposition 5.1, and the spectral characterization of Chapter 3 are clean, well-explained contributions. The empirical program is broad and honestly compared against conjugate gradients, SVGP, and graph neural network baselines. The central limitation is that the claim of inference 'in the exact Gaussian process model' is weakened by the fixed random-feature approximation used for prior samples; this is a model approximation rather than a solver tolerance, and it is not subjected to sensitivity analysis.","major_comments":[{"comment":"The claim that iterative methods perform 'approximate inference in the exact Gaussian process model up to a specified numerical tolerance' conflates solver tolerance with model approximation. In Eq. (3.4) and Eqs. (5.7)–(5.8), every pathwise posterior sample and every probe vector is constructed from a prior sample f_X = Φ_X w with Φ_X Φ_X^T ≈ K_XX using m = 2000 random Fourier features. For finite m, the law of the posterior sample is not the exact GP posterior; this is an approximation of the model itself, not merely of the linear solve. Figure 5.5 confirms that the pathwise marginal-likelihood trajectory deviates from exact optimization and attributes the deviation to random features. Since m = 2000 is fixed for all datasets and kernels with no sensitivity analysis, the 'exact GP' guarantee is unsupported. Please provide either (i) a quantitative bound on ∥K - ΦΦ^T∥ or on the induced posterior perturbation in terms of m and the kernel spectrum, (ii) an m-ablation for representative datasets, including a slow-spectral-decay kernel and higher-dimensional inputs, showing that posterior samples and hyperparameter gradients are insensitive, or (iii) a revised claim that the method performs inference in the random-feature GP model rather than the exact GP model.","section":"Section 2.2.4 / 3.2.2 / 5.2.4"},{"comment":"The inducing-point sampling objective replaces f_X^[Z] ∼ N(0, K_XZ K_ZZ^{-1} K_ZX) with f_X ∼ N(0, K_XX), immediately after Eq. (3.24). This changes the sampling target, and the only justification is the unquantified statement that the approximation error is small 'when m is large and the inducing points are close enough to the data.' The sublinear-cost inducing-point results in Figure 3.2 rely on this replacement, but no bound or controlled experiment isolates its effect. Please provide a quantitative error analysis, or a comparison for moderate n between sampling with exact f_X^[Z] and with f_X, to show the replacement is benign in the regimes used.","section":"Section 3.2.3"},{"comment":"Proposition 5.2, which is used to justify that warm starting introduces only negligible bias, is stated informally as 'Under reasonable assumptions' and its formal proof is deferred to Lin et al. (2024c, Appendix A), which is not included in the dissertation. Warm starting fixes the probe vectors across outer-loop steps, so gradient estimates become coupled and the objective being optimized is no longer exactly the marginal likelihood; this is precisely the step that supports the 72× speed-up claim in Table 5.1. The dissertation should either state the formal theorem with explicit assumptions and a proof, or provide a stronger empirical verification, such as multiple seeds with different fixed probe vectors and a comparison of the final marginal likelihood values against exact optimization. The trajectory histograms in Figures 5.5 and 5.8 are suggestive but do not by themselves establish convergence to the exact optimum.","section":"Section 5.3.2"}],"minor_comments":[{"comment":"In the initialization line, 'α0 = 0' appears twice; presumably the third assignment should be to the geometric average \\bar{α}_0. Please correct this.","section":"Algorithm 4.1"},{"comment":"The abstract claims 'state-of-the-art performance' without qualification, while Chapter 3 itself states the result holds 'on sufficiently large-scale or ill-conditioned regression tasks.' The abstract should carry the same qualification.","section":"Abstract / Section 1.1"},{"comment":"The variance comparison in Proposition 5.1 assumes Gaussian probe vectors for both the standard and pathwise estimators. Since Rademacher or other probe distributions are common in Hutchinson estimators, the text should state explicitly that the comparison is restricted to Gaussian probes.","section":"Section 5.2.3"},{"comment":"In the SVGP negative log-likelihood row for ELEVATORS, the entry '0.43, ± 0.00' contains a stray comma between the value and the standard error.","section":"Table 3.1"}],"recommendation":"major_revision","confidential_remarks":"The dissertation is a compilation of the author's own published papers. For a journal submission, the marginal novelty over the underlying papers is a matter for the editor to weigh. I also note that the 'state-of-the-art' empirical claims are supported by the author's own implementations compared against baselines chosen in the original papers; this is not a correctness issue, but it should be considered when assessing the strength of the empirical contributions."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Short version: this is a competent dissertation-length synthesis of the author's own peer-reviewed papers on iterative methods and pathwise conditioning for scalable GPs. If you already know the NeurIPS/ICLR/ICML papers, there is no new result here; the contribution is the unified presentation, the connections drawn between stochastic dual descent, pathwise estimators, warm starting, and Kronecker structure, and the clean exposition of the derivations. That is a real contribution for someone entering the area, and the empirical comparisons against CG, SVGP, and GNN baselines on datasets up to five million points are useful and honestly scoped.\n\nThe soft spots are mostly about framing. The persistent claim that iterative methods give approximate inference in the exact GP model is too strong. The linear system solves are iterative, but every prior sample, posterior sample, and pathwise gradient uses a random-feature approximation with m=2000 features throughout the whole thesis. That approximation is inherited by the posterior and by the hyperparameter gradients in Chapter 5, and no sensitivity analysis is given for m across different datasets or kernels. For slowly decaying spectra or higher-dimensional inputs, 2000 features can leave real error. The inducing-point variant in Section 3.2.3 also replaces f_X^{[Z]} with f_X based on a non-quantified closeness assertion. So the central 'exact model' language is unsupported; what is actually delivered is approximate inference in an approximate prior model, with the approximation controlled partly by m and partly by solve tolerance. That does not kill the empirical claims--the comparisons to SVGP and CG are apples-to-apples in using the same prior approximations--but it should be stated plainly.\n\nA smaller practical issue: the warm-start bias result in Chapter 5 is deferred to an external appendix, and the thesis ships no code or data. For a dissertation this is normal; for a standalone arXiv paper it reduces reproducibility. I also note some of the 'state-of-the-art' claims lean on the author's own implementations against baselines chosen in the same papers, which is not a flaw per se, but readers should check the underlying papers for independent comparisons.\n\nWho is this for? A graduate student or researcher wanting a single, readable account of where iterative GP methods and pathwise conditioning currently stand. I would cite the underlying papers rather than the thesis, but I would recommend the thesis to students. Deserves a serious referee: yes, though the referee should ask for the m-sensitivity analysis and for the missing proof to be included before the exact-GP claim is advertised.","headline":"A well-written thesis that usefully unifies the author's published iterative-GP work, but the 'exact GP' framing is overstated because it rests on a fixed 2000-random-feature prior approximation.","tokens_in":57555,"tokens_out":1591,"would_cite":false,"duration_ms":22843,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"A dissertation argues that rewriting Gaussian process computations as linear systems and pathwise sampling scales them to millions of data points with order-of-magnitude speed-ups.","keywords":["Gaussian processes","scalable inference","iterative linear system solvers","pathwise conditioning","stochastic dual descent","marginal likelihood optimisation","Kronecker structure","random Fourier features"],"falsifier":"Run the full pipeline (pathwise sampling, stochastic dual descent, warm-started pathwise marginal likelihood gradients) on a regression problem small enough for exact Cholesky inference (say $n=10^4$) using both the standard 2000 random Fourier features and exact prior samples or a much larger feature budget, while solving all linear systems to machine precision; if predictive log-likelihood, posterior sample moments, or the converged marginal likelihood differ systematically between the two, the random-feature approximation, not solver non-convergence, is the limiting assumption.","tokens_in":56547,"feed_emoji":"⚡","tokens_out":9249,"duration_ms":97366,"temperature":0.7,"pith_summary":"This dissertation tries to establish that Gaussian process inference can be made scalable to large modern datasets without moving to an approximate model. Its strategy is to write every expensive GP computation, posterior means, posterior samples, marginal-likelihood gradients, as the solution of a positive-definite linear system and to solve those systems with iterative solvers whose main operation is matrix-vector multiplication. Pathwise conditioning then turns the same solver output into posterior samples that can be evaluated at arbitrary locations, and the chapters build stochastic gradient and dual-descent solvers, warm-started and pathwise gradient estimators, and latent Kronecker structure on top of this core. If the claims are right, a practitioner gets GP inference in the exact model up to a numerical tolerance, with linear or sublinear scaling in time and memory, demonstrated on regression, Bayesian optimisation, and molecular binding affinity tasks, with speed-ups up to $72\\times$ and datasets up to five million examples.","feed_headline":"Gaussian process inference scales to five million points","feed_subtitle":"Rewriting GP computations as linear systems and pathwise sampling cuts memory and runtimes up to 72x.","key_machinery":"The workhorse identity is $v=(K_{XX}+\\sigma^2I)^{-1}b \\iff v=\\arg\\min_u \\tfrac12 u^\\top (K_{XX}+\\sigma^2I)u - u^\\top b$, which turns GP computations into convex quadratic optimisation; for any iterative solver the gradient is the residual $Au-b$, so matrix-vector multiplication with the kernel matrix is the dominant cost. Pathwise conditioning, $f(\\cdot)|y = f(\\cdot)+K(\\cdot)X(K_{XX}+\\sigma^2I)^{-1}(y-(f_X+\\varepsilon))$, converts a single linear solve into a posterior sample and decouples the expensive data-dependent solve from the number of prediction locations. Random Fourier features provide approximate prior samples $f_X$; the dual objective $\\tfrac12\\|\\alpha\\|^2_{K_{XX}+\\sigma^2I}-\\alpha^\\top b$ with random-coordinate gradients gives stochastic dual descent better curvature than the primal objective; and latent Kronecker structure expresses the observed covariance as the projection of a Kronecker product so that fast matrix multiplication survives missing values and irregular inputs. These pieces are the machinery that carries every chapter's method.","core_discovery":"The central claim is that the apparent cubic barrier of Gaussian processes is not inherent: the same kernel matrix inverse appears in prediction, sampling, and hyperparameter learning, and each appearance can be re-expressed as the solution of $(K_{XX}+\\sigma^2I)v=b$, obtained by iterative optimisation of the equivalent convex quadratic. The dissertation argues that stochastic optimisation is a legitimate solver for these systems: stochastic gradient descent and its stochastic dual descent variant have favourable geometry and implicit-bias properties, so they produce accurate predictions even before full convergence. It further shows that a pathwise gradient estimator for the marginal likelihood both accelerates solver convergence and supplies posterior samples at no extra cost, that warm starting solvers across hyperparameter steps yields large speed-ups with negligible bias, and that projecting a latent Kronecker product keeps fast matrix multiplication available for non-grid data. The author would state the contribution as: iterative methods plus pathwise conditioning make Gaussian process inference practical at scales where exact Cholesky-based inference is impossible, while preserving the exact GP model up to a user-chosen numerical tolerance.","pith_inferences":["Beyond the paper: because the solver and sampler are decoupled, replacing random Fourier features with a higher-fidelity approximate prior sampler should remove the main approximation bias without changing the iterative or pathwise machinery; this is directly testable on the dissertation's own benchmarks.","Beyond the paper: the observation that low relative residual norms do not reliably predict predictive performance suggests that convergence criteria for GP linear systems should be tied to downstream quantities such as predictive log-likelihood or acquisition value rather than residual alone.","Beyond the paper: warm-started, pathwise marginal-likelihood gradients and latent Kronecker structure are complementary, so the large-scale hyperparameter learning of Chapter 5 should extend directly to the five-million-example non-grid setting of Chapter 6.","Beyond the paper: the dual-objective insight applies to kernel ridge regression generally, which may shift large-scale kernel methods toward first-order stochastic solvers on ill-conditioned problems where conjugate gradients struggle."],"forward_implications":["Posterior means and posterior samples can be computed with asymptotically linear (or with inducing points, sublinear) time and memory in dataset size, replacing cubic Cholesky-based inference.","Iterative and pathwise inference targets the exact Gaussian process model up to a numerical tolerance rather than an approximate sparse or variational model, preserving calibrated uncertainty for decision-making.","Because pathwise conditioning requires one linear solve per posterior sample independent of evaluation locations, acquisition-function optimisation in Bayesian optimisation can reuse samples across many candidate points.","Hyperparameter learning can be accelerated by a pathwise marginal-likelihood gradient estimator and warm-started solvers, yielding up to $72\\times$ speed-ups when solving to tolerance and up to $7\\times$ lower residual norms on a fixed budget.","Latent Kronecker structure extends scalable GP inference to non-grid data, demonstrated on real-world datasets with up to five million examples in robotics, automated machine learning, and climate modelling."],"supporting_citations":[{"why":"Supplies the pathwise conditioning formula used throughout the dissertation to generate posterior samples from prior samples and linear solves.","marker":"Wilson et al. (2020, 2021)"},{"why":"Supplies random Fourier features used to approximate prior function samples and to estimate kernel-norm regularisers.","marker":"Rahimi and Recht (2008)"},{"why":"Established the conjugate-gradient framework for scalable GP marginal likelihood and posterior inference that the dissertation extends and benchmarks against.","marker":"Wang et al. (2019)"},{"why":"Provides the sparse variational GP whose inducing-point posterior the stochastic-gradient extension is designed to approximate.","marker":"Titsias (2009a)"},{"why":"Provides stochastic variational inference with inducing points, a principal comparison baseline in the empirical chapters.","marker":"Hensman et al. (2013)"},{"why":"Contributes stochastic dual coordinate ascent, the algorithmic inspiration for the stochastic dual descent objective and updates.","marker":"Shalev-Shwartz and Zhang (2013)"},{"why":"Contributes the low-variance stochastic gradient estimator and variance reduction used for posterior sampling objectives.","marker":"Antorán et al. (2023)"},{"why":"Provides the DOCKSTRING molecular binding affinity benchmark and preprocessing used to compare GP regression with graph neural networks.","marker":"García-Ortegón et al. (2022)"},{"why":"Establishes preconditioned conjugate-gradient inference for GPs and is a baseline for linear-system solver performance.","marker":"Gardner et al. (2018a)"}],"fun_headline_variants":["GP inference as iterative linear solves: 5M points","Rewrite GP ops as solves, scale past five million","Pathwise conditioning gets GPs to five million","No Cholesky needed: GPs scale via iterative solves","Solvers and pathwise sampling push GPs to 5M"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise is that the approximate prior samples used for pathwise conditioning and marginal-likelihood gradients, drawn from a random-feature model with a fixed number of features rather than from the exact GP, are faithful enough; if they are not, every posterior sample and hyperparameter gradient inherits the bias.","fun_headline_variants_meta":{"raw":{"variants":["GP inference as iterative linear solves: 5M points","Rewrite GP ops as solves, scale past five million","Pathwise conditioning gets GPs to five million","No Cholesky needed: GPs scale via iterative solves","Solvers and pathwise sampling push GPs to 5M"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000823,"raw_usage":{"total_tokens":3561,"prompt_tokens":869,"completion_tokens":2692,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":485,"completion_tokens_details":{"reasoning_tokens":2611}},"tokens_in":485,"tokens_out":2692,"duration_ms":23066,"temperature":1.0,"reasoning_tokens":2611,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-06T18:53:03.083715+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run the full pipeline (pathwise sampling, stochastic dual descent, warm-started pathwise marginal likelihood gradients) on a regression problem small enough for exact Cholesky inference (say $n=10^4$) using both the standard 2000 random Fourier features and exact prior samples or a much larger feature budget, while solving all linear systems to machine precision; if predictive log-likelihood, posterior sample moments, or the converged marginal likelihood differ systematically between the two, the random-feature approximation, not solver non-convergence, is the limiting assumption.","supporting_citations":[],"review_version":1}