{"id":"927969ed-523e-42db-9a06-de7194196077","arxiv_id":"2505.16481","paper_version":3,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":2,"one_line_summary":"Two neighbour-based approximations, HPA and SPA, make Gaussian Process variational autoencoders scalable while preserving local latent correlations.","lead":"This paper makes Gaussian Process variational autoencoders scalable by approximating the latent prior with two nearest-neighbour schemes, HPA and SPA. It reports strong speed and accuracy on video, physics, and geostatistics benchmarks, with one benchmark where the claimed dominance does not hold.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Nearest-neighbour truncation is only faithful for local kernels; periodic/global latent structure breaks the core 'preserves essential dependencies' claim, so the flexible-kernel claim needs a scope restriction.","rationale":"The reader's weakest assumption identifies exactly the load-bearing premise: local adjacency in the auxiliary input space must capture essentially all GP-prior dependence that matters. I agree with that diagnosis. The paper's strongest claim, that the method 'preserves essential latent dependencies' while 'allowing more flexible kernel choices', is universal in tone, and Table 15 explicitly marks arbitrary-kernel support for both HPA and SPA. Existing NNGP theory justifies neighbour truncation under ordering and conditional-independence assumptions that hold for kernels with decaying local dependence; they do not hold for periodic or other globally coupled kernels, where the most informative conditioning set is not the Euclidean nearest-neighbour set. A synthetic periodic-GP experiment would settle this cleanly because the failure would be a statistical misspecification of the prior, not a tuning or optimisation artefact. I deliberately focus on this rather than the HPA ELBO gap in Appendix B.1: that gap is real and should be fixed, but it affects the theoretical status of one variant, whereas the locality premise is shared by both variants and by the central claim itself. The reader's CONDITIONAL verdict remains appropriate: the method is scalable and empirically strong on local-adjacency data, but its central claim needs a stated restriction to kernels and datasets where dependence is dominantly local.","tokens_in":26985,"tokens_out":6260,"duration_ms":64590,"concrete_test":"Run a controlled synthetic experiment: draw latent trajectories from a periodic GP with period T such that T > H times the typical spacing (e.g., T=50, H=10, spacing 1) and generate observations through a fixed nonlinear decoder. Train GPVAE-SPA and GPVAE-HPA with Euclidean nearest temporal neighbours, a standard VAE, a full-batch GPVAE, and a modified SPA whose conditioning set includes the H period-lagged points. Compare test NLL and RMSE. If the Euclidean-H models are no better than the VAE while the period-lagged variant matches the full-batch GPVAE, the local-adjacency premise is falsified as a general claim; if the Euclidean-H models track the full-batch baseline, the concern does not land.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central premise (Section 3.2) is that conditioning on the H nearest neighbours of each location preserves the GP-prior dependence that matters. This is not a property of flexible kernels in general: it is a conditional-independence assumption. For a periodic kernel such as k(x,x') = sigma^2 exp(-2 sin^2(pi(x-x')/T)/ell^2) with period T larger than the H-neighbour window, the nearest neighbours are nearly uncorrelated while points one period apart are strongly correlated, so the H-neighbour prior (12) and the predictive rule in Section 3.3 discard the only strong dependencies. The reported benchmarks use RBF, Matern-3/2 and Cauchy kernels on temporal or spatial grids where correlation decays with Euclidean distance, so they never exercise this failure regime. Consequently the abstract's claim that the method 'preserves essential latent dependencies' and 'allows more flexible kernel choices' is not supported for global or periodic structure; the claim should be scoped to kernels whose dependence is dominantly local.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes two neighbour-driven approximations to the Gaussian-process prior in GPVAEs: Hierarchical Prior Approximation (HPA), Eq. (10), and Sparse Precision Approximation (SPA), Eq. (13). Both restrict prior dependencies to the H nearest neighbours of each auxiliary location, giving O(N_b H^3) training complexity and enabling mini-batch inference with arbitrary kernel functions. The authors evaluate the methods on latent trajectory reconstruction, rotated-MNIST imputation and generation, MuJoCo conditional generation, and two geostatistical datasets, reporting improvements over several GPVAE baselines in prediction quality and wall-clock time.","tokens_in":27144,"tokens_out":12820,"duration_ms":111766,"significance":"If the claims hold, the paper provides a practical alternative to inducing-point-based GPVAEs for structured latent modelling on data with local dependence, with an open-source implementation and a careful set of benchmarks. The SPA construction is a standard Vecchia-style approximation and is coherent as a scalable surrogate for the full GP prior. The reported experiments are extensive, use 10 random trials, and include both small and large datasets, which makes the empirical contribution valuable. The main caveats are that the HPA objective is not derived as a valid ELBO for the stated hierarchical model, and that the nearest-neighbour premise restricts the class of kernels for which the approximation is statistically faithful.","major_comments":[{"comment":"The displayed objective L_HPA is not presented as a valid ELBO for the hierarchical model in Eqs. (6)-(9). In the derivation, the full likelihood term E_{p(w)} E_{q(Z|w)} log p(Y|Z) is replaced by N/|I| times a sum over the mini-batch I, and the full KL term is replaced by an average over neighbourhood KL divergences. When the neighbourhoods n(i) overlap, the KL term counts overlapping latent variables multiple times, while the likelihood term ignores selected neighbours that are not in I. In addition, q(z_i|w) from Eq. (8) is replaced by q(z_i|y_i), which is not justified. Thus L_HPA is, for finite H, a block-wise surrogate objective rather than a lower bound on log p(Y). The statement that H=N recovers Eq. (4) only establishes agreement in the limit and does not validate the finite-H objective. The authors should either define a variational family and a proper distribution p(w) for which Eq. (10) is a genuine ELBO, or explicitly present HPA as an approximate training objective and analyse the bias introduced by the block approximation.","section":"Appendix B.1, Eq. (10)"},{"comment":"The claims that the method 'preserves essential latent dependencies' and supports 'arbitrary kernels' are scoped too broadly. The H-neighbour truncation in Eqs. (10) and (13) is a conditional-independence assumption: it is faithful only when kernel-induced dependence is dominated by short-range correlations. For a periodic kernel such as k(x,x') = sigma^2 exp(-2 sin^2(pi(x-x')/T)/ell^2) with period T larger than the H-neighbour window, the nearest neighbours are nearly uncorrelated while points one period apart are strongly correlated; the H-neighbour prior and the predictive rule in Section 3.3 would discard the dominant dependencies. Every experiment in the paper uses RBF, Matérn-3/2, or Cauchy kernels on temporal or spatial grids with monotone distance decay, so this failure regime is not exercised. The claims should be restricted to kernels with dominantly local correlation structure, or a long-range/periodic-kernel experiment should be added to demonstrate where the approximation breaks down.","section":"Abstract, Section 3.2, Section 3.3, Table 15"},{"comment":"The central empirical claim that the approach 'improves both predictive accuracy and training speed compared to existing GPVAE baselines' is contradicted by the Rotated-MNIST corrupted-frame results in Table 1. There, MGPV AE achieves NLL 0.090 +/- 0.001 versus 0.095 +/- 0.000 for HPA and 0.096 +/- 0.000 for SPA, RMSE 0.197 +/- 0.001 versus 0.199 and 0.202, and a faster training time of 110.8 s/epoch versus 137.9 and 121.2 s/epoch. The paper acknowledges this in Section 5.2, but the Abstract and Section 1 still make an unqualified superiority claim. The superiority claims should be qualified to the settings and metrics where they actually hold.","section":"Introduction and Table 1"}],"minor_comments":[{"comment":"The complexity statement 'Locating the H nearest neighbours of each point takes O(HN) in the worst case' is dimensionally incomplete: a per-query cost of O(HN) gives O(H N^2) for all N points, and the Faiss-based approximate search does not have this worst-case guarantee as stated.","section":"Section 3.4"},{"comment":"The sentence 'setting H = 0 will cause the model to degenerate into conventional VAEs' is only true when the marginal variances k(x_j, x_j) are 1; in general H=0 gives a VAE with a diagonal Gaussian prior with per-point variances. Please clarify.","section":"Section 3.2, SPA discussion"},{"comment":"The predictive posterior is derived using the SPA-style conditional p(z_* | Z_n(*)); for HPA, whose approximate posterior is not defined through a chain factorisation, the same formula is used without justification. Please state explicitly that this is a shared approximation for both variants.","section":"Section 3.3"},{"comment":"The notation for the variational distribution is inconsistent: Eq. (10) writes E_{q(z_i|y_i)}, while the derivation uses E_{q(z_i|w)}. Please align these notations and define the marginal q(z_i) used in the final objective.","section":"Appendix B.1 and Eq. (10)"}],"recommendation":"major_revision","confidential_remarks":"The stress-test concern about local-only fidelity lands: the nearest-neighbour premise is the load-bearing statistical assumption, and the paper's broad 'arbitrary kernel' claim needs to be scoped. The HPA derivation is the other substantive issue; if the authors reframe HPA as an approximate surrogate objective rather than a rigorous ELBO, the paper's empirical contribution remains viable. I would be willing to re-review a revised version."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"What you should know: this is a legitimate and reproducible extension of NNGP ideas to the latent space of GPVAEs. The SPA variant is essentially a clean Vecchia approximation of the GP prior with amortised mini-batch inference, and the experiments are thorough enough that the scalability story mostly checks out. The abstract overclaims, though, and the HPA derivation does not actually establish an ELBO.\n\nThe genuinely new bit is applying nearest-neighbour GP approximations—Vecchia and sparse precision—not to GP regression but to the latent prior of a VAE, with a standard encoder and mini-batch training. The code is open-sourced, the experimental detail is unusually complete, and the SPE10 result with 1.4e5 points is a real demonstration that this buys scalability without the inducing-point headaches. On the moving-ball and MuJoCo tasks the neighbour-driven models are competitive with or better than full-batch and inducing-point baselines. That part should be credited.\n\nNow the soft spots, in rough order of severity. First, the central claim that this 'preserves essential latent dependencies' and 'allows more flexible kernel choices' is too broad. Conditioning on H nearest neighbours is a conditional-independence assumption. It is faithful only when the kernel's dependence is dominantly local. For a periodic kernel with period longer than the neighbour window, or any global/long-range coupling, the nearest neighbours are nearly independent while far-away points carry the strong correlation. All the benchmarks use RBF, Matern-3/2, or Cauchy kernels on temporal or spatial grids, so the failure regime is never exercised. The claim should be scoped to locally-correlated kernels.\n\nSecond, the HPA objective in Eq. (10) is presented as an ELBO for the hierarchical prior in Eqs. (6)-(9), but Appendix B.1 replaces the expectation over p(w) with a specific neighbour-selection rule without showing the resulting expression is still a lower bound. The last step of the derivation looks like an approximation, not an equality. That's a genuine gap. SPA's derivation is fine.\n\nThird, the abstract's 'outperforms other GPVAE variants in predictive performance' is contradicted by Table 1: on corrupted-frame MNIST, both HPA and SPA have worse NLL and RMSE than MGPV AE, and GPV AE-Band has lower RMSE. The 'computational efficiency' part of the claim holds, but the universal predictive outperformance does not. Minor point: the SPA ordering for spatial data is unspecified; for 2D/3D coordinates the chain ordering can matter, and it's worth stating.\n\nWho this is for: anyone working on scalable GPVAEs or structured latent models with local dependence. It deserves a serious referee; the empirical core is reproducible and the method is useful in its proper scope, but the flexible-kernel claim and the HPA bound need fixing before I'd trust it as stated.","headline":"A solid, reproducible application of Vecchia-style nearest-neighbour conditioning to GPVAE latent spaces, with an overstated flexible-kernel claim and a hand-wavy HPA derivation.","tokens_in":27722,"tokens_out":2802,"would_cite":false,"duration_ms":23173,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["68T07","62M30"],"pacs":[],"model":"deepseek-v4-flash","headline":"This paper claims that replacing the dense Gaussian-process prior over latent variables with a nearest-neighbour-truncated prior makes Gaussian-process variational autoencoders scalable and kernel-flexible, competitive with or better than…","keywords":["Gaussian process variational autoencoder","nearest neighbour Gaussian process","scalable variational inference","structured latent models","sparse precision approximation","hierarchical prior approximation","mini-batch training","data imputation"],"falsifier":"Train GPVAE-SPA and GPVAE-HPA on a synthetic temporal dataset whose latent trajectories are drawn from a GP with an RBF kernel whose lengthscale is several times the observation window, so every point correlates strongly with every other point. If, with $H$ fixed at a small value such as 10, the latent-trajectory reconstruction RMSE is far worse than a full-batch GPVAE, while the same models match the full-batch on short-lengthscale trajectories, the local-adjacency premise is false for long-range dependence. A quantitative precursor would be the KL divergence between the true $N$-dimensional GP prior and the $H$-neighbour approximation, which should stay large for small $H$ under a long lengthscale.","tokens_in":26751,"feed_emoji":"🧠","tokens_out":12602,"duration_ms":87866,"temperature":0.7,"pith_summary":"The paper sets out to remove the main scalability bottleneck of Gaussian Process Variational Autoencoders (GPVAEs): the dense $N\\times N$ covariance matrix of the GP prior, which makes exact inference $O(N^3)$ and blocks mini-batching. Its central claim is that a local-neighbourhood truncation of that prior — keeping only the $H$ nearest neighbours of each auxiliary location — preserves the latent correlations that matter for structured data while cutting per-step cost to $O(N_b H^3)$. Two concrete approximations are proposed and tested: Hierarchical Prior Approximation (HPA), which switches off interactions between non-neighbours through a latent selection variable, and Sparse Precision Approximation (SPA), which factorises the prior into conditional distributions on the $H$ nearest neighbours. On representation learning, imputation, and conditional generation tasks (including a spatial dataset with about 140,000 points), the paper reports that both variants match or improve on full-batch and inducing-point GPVAE baselines in predictive error and training speed.","feed_headline":"Local neighbours scale GP variational autoencoders to 100k points","feed_subtitle":"Cost scales with the neighbour count, not the dataset size, and accuracy matches or beats sparse-GP baselines.","key_machinery":"The load-bearing object is the $H$-neighbour truncation of the latent GP prior, in two complementary forms. HPA (Hierarchical Prior Approximation) imposes sparsity on the covariance: a binary selection vector $w$ activates only neighbourhood-sized blocks $D_w K_{XX} D_w$, so each mini-batch's KL uses an $H\\times H$ covariance. SPA (Sparse Precision Approximation) imposes sparsity on the precision: the joint GP prior is factorised via the probability chain rule and each factor $p(z_j \\mid z_{n(j)})$ conditions only on the $H$ nearest neighbours of $x_j$ among previously ordered points, giving a sparse precision/Cholesky structure. Both replace the dense $N\\times N$ covariance or precision with $H\\times H$ pieces, so the per-iteration cost is $O(L N_b H^3)$ once nearest neighbours are precomputed, and both reduce to the full GPVAE ELBO when $H=N$ and to a plain VAE when $H=0$. The same local conditioning is used at prediction time, so unseen locations are handled by their $H$ neighbours only.","core_discovery":"The paper's contention is that the GP prior over latent variables in a GPVAE can be approximated by a neighbour-driven sparsity pattern without losing the structured dependence the prior is there to encode. HPA multiplies the dense GP covariance by a binary diagonal selection matrix $D_w$, so that only latent variables belonging to the same selected neighbourhood interact, yielding an ELBO whose KL term splits into $H\\times H$ sub-blocks. SPA instead writes the prior as a product of conditional densities $p(z_j \\mid z_{n(j)})$, where $n(j)$ collects the $H$ nearest neighbours of $x_j$ among earlier points, which is a sparse Cholesky-type approximation of the precision matrix; its ELBO again separates into mini-batches. In both cases the training objective is amortised over mini-batches and accepts arbitrary kernels, and prediction at an unseen location $x_*$ conditions only on the $H$ nearest neighbours of $x_*$. The paper presents experiments on moving-ball trajectories, rotated handwritten-digit sequences, physics-simulation action series, and two geostatistical datasets, reporting that the two approximations approach or exceed the accuracy of full-batch GPVAE baselines and of inducing-point models that use comparable or larger budgets, while training faster.","pith_inferences":["The choice of distance metric for 'nearest' is a hidden design parameter; the paper uses Euclidean distance throughout. Swapping in a kernel-induced or manifold-aware distance could change both HPA and SPA behaviour on data whose local structure is curved, such as rotation manifolds, and would be a direct test of how much the neighbourhood definition matters.","Because both approximations interpolate between a factorised VAE at $H=0$ and the exact GP prior at $H=N$, the gap between their curves as $H$ grows is a measurable diagnostic of how locally structured a dataset actually is; datasets where small-$H$ performance is already near the full-batch result are precisely those where the method's premise holds.","The same trick should transfer to other structured priors: any latent model whose prior is a Gaussian Markov random field or a graph GP could use graph-adjacency neighbourhoods instead of coordinate distance, extending the scalability claim beyond spatial and temporal auxiliaries.","The consistent small edge of SPA over HPA in the paper's tables suggests that sparsifying the precision (conditional independence) may be a more faithful local approximation than sparsifying the covariance (zeroing non-neighbour covariances); comparing the two on a long-range-kernel dataset would show whether that ordering generalises."],"forward_implications":["Training scales to large structured datasets: with $H$ fixed at, e.g., 20, the per-mini-batch cost $O(N_b H^3)$ is independent of total $N$, and the paper demonstrates a 140,000-point spatial imputation running in seconds per epoch.","Arbitrary kernels become usable in the latent prior, including RBF, Matérn, and Cauchy kernels, because the approximation needs only pairwise kernel evaluations on neighbour sets and avoids the low-rank or state-space kernel restrictions of earlier GPVAE variants; the paper's experiments lean on this flexibility.","Inducing points are no longer the resolution mechanism: on the moving-ball data, the paper reports that the sparse-precision variant reaches good reconstruction with about one sixth of the trajectory points as neighbours, where an inducing-point baseline needed roughly half the trajectory as inducing inputs to approach the full-batch result.","With $H=N$ the objectives recover the full-batch GPVAE ELBO, and with $H=0$ the SPA objective degenerates to a standard VAE, so the method forms a continuous family from factorised to fully structured latent priors.","Conditional generation at an unseen timestamp or location is cheap: prediction at $x_*$ uses only the $H$ nearest neighbours of $x_*$, a small encoder pass, and one $H\\times H$ linear solve."],"supporting_citations":[{"why":"Introduces the GPVAE and its ELBO with a non-factorised KL over the full GP prior — the computational target this paper replaces.","marker":"(Casale et al., 2018)"},{"why":"Origin of the conditional-factorisation approximation of spatial likelihoods that SPA adapts to the latent GP prior.","marker":"(Vecchia, 1988)"},{"why":"Neighbour-based likelihood approximation for large spatial datasets; cited as the basis of the local-adjacency intuition.","marker":"(Stein et al., 2004)"},{"why":"Hierarchical nearest-neighbour GP construction whose sparse covariance structure inspires HPA's selection-variable mechanism.","marker":"(Datta et al., 2016a)"},{"why":"Review of nearest-neighbour GP models for massive spatial data; establishes the precedents for neighbour-driven GP scalability.","marker":"(Datta et al., 2016b)"},{"why":"Sparse-within-sparse variational GP with neighbour information and a binary indicator prior; the closest analogue to HPA.","marker":"(Tran et al., 2021b)"},{"why":"Variational nearest-neighbour GP (VNNGP) with sparse precision over inducing variables; supplies the SPA mechanism and is a direct experimental baseline.","marker":"(Wu et al., 2022)"},{"why":"SVGPAE, the principal inducing-point baseline; the paper's running comparison for how many inducing points are needed versus H neighbours.","marker":"(Jazbec et al., 2021)"},{"why":"Full-batch GPVAE-Diag and GPVAE-Band baselines for missing-pixel imputation; Diag is the gold-standard full-batch analogue of the proposed models.","marker":"(Fortuin et al., 2020)"},{"why":"Markovian GPVAE (MGPVAE) baseline using state-space Matérn kernels; compared on rotated handwritten-digit and physics-simulation conditional generation.","marker":"(Zhu et al., 2023)"}],"fun_headline_variants":["Nearest neighbours scale GP variational autoencoders","Neighbour-driven sparsity scales GPVAEs without inducing points","Local neighbours cut GPVAE cost to fit 100k points","GPVAEs go large with neighbour-based approximations"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise is that the $H$ nearest neighbours of each point in the auxiliary space carry essentially all of the GP-prior dependence that matters, so that correlations to everything further away can be zeroed out without distorting the latent structure.","fun_headline_variants_meta":{"raw":{"variants":["Nearest neighbours scale GP variational autoencoders","Neighbour-driven sparsity scales GPVAEs without inducing points","Local neighbours cut GPVAE cost to fit 100k points","GPVAEs go large with neighbour-based approximations"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000614,"raw_usage":{"total_tokens":2859,"prompt_tokens":955,"completion_tokens":1904,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":571,"completion_tokens_details":{"reasoning_tokens":1838}},"tokens_in":571,"tokens_out":1904,"duration_ms":13084,"temperature":1.0,"reasoning_tokens":1838,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-07T14:59:57.571508+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Train GPVAE-SPA and GPVAE-HPA on a synthetic temporal dataset whose latent trajectories are drawn from a GP with an RBF kernel whose lengthscale is several times the observation window, so every point correlates strongly with every other point. If, with $H$ fixed at a small value such as 10, the latent-trajectory reconstruction RMSE is far worse than a full-batch GPVAE, while the same models match the full-batch on short-lengthscale trajectories, the local-adjacency premise is false for long-range dependence. A quantitative precursor would be the KL divergence between the true $N$-dimensional GP prior and the $H$-neighbour approximation, which should stay large for small $H$ under a long lengthscale.","supporting_citations":[{"cited_title":"P., Dalca, A., Saglietti, L., Listgarten, J., and Fusi, N","cited_arxiv_id":null,"evidence_quote":"Introduces the GPVAE and its ELBO with a non-factorised KL over the full GP prior — the computational target this paper replaces."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Origin of the conditional-factorisation approximation of spatial likelihoods that SPA adapts to the latent GP prior."},{"cited_title":"L., Chi, Z., and Welty, L","cited_arxiv_id":null,"evidence_quote":"Neighbour-based likelihood approximation for large spatial datasets; cited as the basis of the local-adjacency intuition."},{"cited_title":"Scalable Gaussian process variational autoencoders","cited_arxiv_id":null,"evidence_quote":"SVGPAE, the principal inducing-point baseline; the paper's running comparison for how many inducing points are needed versus H neighbours."},{"cited_title":"Gp-vae: Deep probabilistic time series imputation","cited_arxiv_id":null,"evidence_quote":"Full-batch GPVAE-Diag and GPVAE-Band baselines for missing-pixel imputation; Diag is the gold-standard full-batch analogue of the proposed models."},{"cited_title":"Markovian Gaussian process variational autoencoders","cited_arxiv_id":null,"evidence_quote":"Markovian GPVAE (MGPVAE) baseline using state-space Matérn kernels; compared on rotated handwritten-digit and physics-simulation conditional generation."}],"review_version":1}