{"id":"d65f41c8-55c6-484e-ae0d-0e39082eb17a","arxiv_id":"2512.18066","paper_version":2,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":4,"one_line_summary":"Gradient-enhanced deep Gaussian processes, built by sampling latent warpings and their derivatives via MCMC with a chain rule, outperform gradient-enhanced GPs and standard DGPs on nonstationary test functions.","lead":"This paper upgrades deep Gaussian process emulators so they can use derivative information from expensive simulations and also output derivative predictions. That makes surrogate modeling of nonstationary simulators more accurate and supports optimization and active learning.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Eq. (16) omits the change-of-variables Jacobian for the transformation ∇_x y → ∇_w y, so the ESS sampler does not target the posterior of Eq. (11); this threatens the Bayesian validity of the central claim.","rationale":"The reader's weakest_assumption identifies exactly the missing Jacobian in Eq. (16) and the invertibility of ∇_x W. My analysis confirms that the change of variables from ∇_w y to ∇_x y requires a determinant factor, and because this determinant varies with the sampled W_all, omitting it means the chain does not target the posterior of Eq. (11). This is a genuine correctness defect in the core Bayesian inference, so the paper cannot be accepted as a fully valid Bayesian framework without addressing it. However, the issue is local and fixable: adding the Jacobian term and rerunning the benchmarks would determine whether the empirical conclusions survive. The reader's CONDITIONAL verdict is therefore appropriate; I do not see a reason to escalate to REJECT, since the method may still be practically effective despite the statistical misspecification. I also note the abstract promises real quantum mechanics experiments that do not appear in the body, and Yang et al. (2025) is not benchmarked, but these are secondary to the Jacobian problem. My verdict is UNCHANGED, meaning it stays CONDITIONAL pending the Jacobian fix and empirical re-check.","tokens_in":18682,"tokens_out":7585,"duration_ms":68651,"concrete_test":"Modify the ESS acceptance probability in Section 3.3 to include the additive term -Σ_i log|det(∇_x W_i)| in the log-likelihood of Eq. (16). Rerun the squiggle experiment (D=2, n=25, 30 Monte Carlo repeats) with both the original and corrected samplers, comparing predictive RMSE/CRPS for y and ∇_x y, and also compare the posterior distribution of sampled ∇_x W (e.g., its determinant) across the two versions. If the predictive metrics shift by more than a small tolerance (e.g., 10%) or the posterior of ∇_x W changes materially, the omitted Jacobian is consequential.","verdict_should_be":"UNCHANGED","load_bearing_attack":"In Section 3.3, the observed input-space gradients ∇_x y are converted to latent-space gradients ∇_w y by solving Eq. (13), i.e., ∇_x y = ∇_x W · ∇_w y. For each training location i, this is a linear bijection with Jacobian det(∇_x W_i). The correct likelihood for the observations (y, ∇_x y) given the latent warping W_all is therefore\n\np(y, ∇_x y | W_all) = p(y, ∇_w y | W_all) · ∏_i |det(∇_x W_i)|^{-1}.\n\nThe paper's Eq. (16) evaluates the outer-layer likelihood at the solved ∇_w y but omits the log |det(∇_x W_i)| term. Since ∇_x W is part of the sampled W_all, this determinant changes with every proposal (17). Consequently, the stationary distribution of the gradient-enhanced ESS chain is not the posterior from the model in Eq. (11); it is weighted by an extra ∏_i |det(∇_x W_i)| factor, biasing the posterior over warping functions toward regions of high |det(∇_x W)|. This directly undermines the claim that the method provides a comprehensive Bayesian treatment with valid uncertainty quantification. The paper also does not discuss the required invertibility of ∇_x W; while singularities have zero probability under a continuous prior, near-singularity could cause numerical instability and extreme solutions when solving Eq. (13).","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper develops a Bayesian deep Gaussian process (DGP) framework that incorporates gradient information. The model places gradient-enhanced Gaussian priors on both the outer layer (response and response gradients with respect to latent warped inputs) and the inner layer (latent warping nodes and their gradients with respect to original inputs). The multivariate chain rule is used to convert between input-space gradients and warped-space gradients. Gradient-enhanced DGP training is performed by augmenting elliptical slice sampling to sample the latent warping and its gradient; DGP gradient predictions are obtained by applying GP gradient prediction at both layers and multiplying by the inferred warping Jacobian. The methodology is implemented in the CRAN package deepgp and benchmarked on three nonstationary functions and a quantum mechanics-inspired example, comparing GP, DGP, gradient-enhanced GP, and gradient-enhanced DGP in terms of RMSE and CRPS for both responses and gradients. The central empirical claim is that on nonstationary surfaces the gradient-enhanced DGP outperforms gradient-enhanced GPs and non-enhanced DGPs.","tokens_in":19107,"tokens_out":8092,"duration_ms":85167,"significance":"If the proposed framework is statistically valid, it fills a real gap: gradient-enhanced DGPs and DGP gradient predictions are useful for expensive nonstationary simulators, especially for downstream optimization and active learning. The paper is also accompanied by open-source software and reproducible benchmark code, which is a strength. However, the validity of the Bayesian inference and the uncertainty quantification claims rests on the correctness of the likelihood used in the MCMC scheme and on the gradient predictive moment formulas. As detailed below, both contain technical gaps that need to be addressed before the results can be taken as Bayesian posterior inference.","major_comments":[{"comment":"The ESS acceptance probability uses the wrong observed-data likelihood. Observed input-space gradients are related to the warped-space gradients by A_i ∇_w y_i = ∇_x y_i, where A_i = ∇_x W_i is D×D. The density of the observations given W_all is therefore p(y, ∇_x y | W_all) = p(y, ∇_w y | W_all) ∏_{i=1}^n |det A_i|^{-1}. Eq. (16) evaluates the Gaussian density at the solved ∇_w y but omits the ∏ |det A_i|^{-1} factor. Because A_i changes at every ESS proposal (17), the stationary distribution is ∝ prior(W_all) p(y, ∇_w y | W_all) ∏ |det A_i|, not the model posterior. This biases the posterior over warpings toward regions with large |det ∇_x W| and can seriously distort uncertainty quantification. The required invertibility of A_i at every training location is also not discussed; near-singular A_i will make the solve in Eq. (13) unstable.","section":"Section 3.3, Eqs. (13) and (16)"},{"comment":"The stated variance of the DGP gradient prediction is not the variance of the linear combination. For a given posterior draw, ∂y/∂x_d = Σ_i (∂w_i/∂x_d) ∂y/∂w_i, so the conditional variance is (∇_x w_d)^T Cov(∂y/∂w_1,...,∂y/∂w_D) (∇_x w_d). The displayed formula V[∂y/∂x_d] = Σ_i (∂w_i/∂x_d)^2 V[∂y/∂w_i] drops all off-diagonal terms of the predictive covariance from Eq. (5). These cross-covariances are generally nonzero for the squared-exponential kernel, so the reported gradient CRPS may not reflect the actual posterior predictive distribution. If the software uses this formula for UQ, the gradient uncertainty statements are not reliable.","section":"Section 3.2, moment formulas"},{"comment":"The method requires the D×D matrix ∇_x W_i to be invertible at every training location i. Under a continuous GP prior, exact singularity has probability zero, but near-singularity is possible and can produce extreme solved values of ∇_w y, especially because the Gaussian kernel gradients can be near zero or near-duplicated. The manuscript does not acknowledge this assumption, discuss numerical safeguards, or report any diagnostics (e.g., condition numbers) in the benchmark experiments.","section":"Section 3.3, Eq. (13)"}],"minor_comments":[{"comment":"Eq. (16) is described as 'the gradient-enhanced version of the likelihood of the outer layer.' It is actually the density of the transformed pseudo-observations (y, ∇_w y), not the likelihood of the observed data (y, ∇_x y). The wording should be revised together with the Jacobian correction.","section":"Section 3.3"},{"comment":"In the sentence introducing Wall, 'column-binded' should be 'column-bound' or 'column-bind.' Similarly, 'intractible' in the introduction and 'surrgotes' in Section 6 are typos.","section":"Section 3.1"},{"comment":"The caption says 'For observed training locations, we sample these (Section 3.3)' referring to dw/dx. Section 3.3 actually solves for ∇_w y using Eq. (13); what is sampled is W_all. Please clarify the wording.","section":"Figure 2 caption"},{"comment":"The Vecchia ordering description 'then append each ∂y_i/∂x_d using the same ordering, in turn' is ambiguous. Does this mean all derivatives with respect to x_1 for all i, then all with respect to x_2, or derivative-by-derivative within each observation? A short pseudocode block or explicit index sequence would remove ambiguity.","section":"Section 4.1"}],"recommendation":"major_revision","confidential_remarks":"The paper fits the journal's scope and the empirical benchmarks are extensive, with open-source code. The missing Jacobian in the ESS likelihood is a substantive statistical error, but it is local and correctable; the gradient prediction variance issue is also correctable. Given both affect the main Bayesian/UQ claims, I recommend major revision rather than rejection. I did not see circularity in the benchmark design."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Bottom line: the idea is genuinely useful and the engineering is careful, but there is a load-bearing Bayesian gap. The ESS acceptance in Section 3.3 plugs the solved \\nabla_w y into the outer-layer Gaussian likelihood as though it were observed, omitting the |det(\\nabla_x W)|^{-1} Jacobian from the change of variables in Eq. (13). Since \\nabla_x W changes with every proposal, the chain does not target the posterior of Eq. (11); it targets a distribution weighted by the missing determinant, biasing warping toward large |det(\\nabla_x W)|. The paper never discusses invertibility of \\nabla_x W either, and near-zero local derivatives would blow up the solved \\nabla_w y. That is a real correctness problem for the \"comprehensive Bayesian framework\" claim, though the fix is conceptually simple: add the log-determinant term, or treat \\nabla_w y as latent and integrate it out.\n\nWhat it does well: the chain-rule decomposition of DGP gradients is clean, and separating gradient-enhancement from gradient prediction is useful. The benchmarks are solid—three nonstationary surfaces, 30 reps, proper LHS testing, RMSE and CRPS—and geDGP consistently beats geGP and DGP. Extending Vecchia to gradient kernels is a practical contribution. Code is shipped in deepgp with a public repo; that's real, reproducible work.\n\nSmaller issues: the abstract promises real quantum mechanics experiments that never appear in the body—only mock ignition and two synthetic functions. That needs to be reconciled or removed. The only DGP-gradient baseline, Yang et al. (2025), is cited but not benchmarked; given the direct comparison, it should be. Invertibility and conditioning of \\nabla_x W also deserve a paragraph.\n\nThis is for people doing nonstationary emulation with gradient sources—adjoint solvers, CFD, molecular energy/force calculations. They'll want to try the package. As a referee, I'd send it out but expect major revision: fix the likelihood, fix the abstract mismatch, ideally add the Yang comparison. The empirical core is probably salvageable; the Bayesian framing is not yet.\n\nRecommendation: serious referee, conditional.","headline":"Useful, well-engineered gradient-enhancement for deep GPs with a real Jacobian omission in the MCMC sampler; conditional accept after fixing the likelihood and the abstract mismatch.","tokens_in":19507,"tokens_out":4321,"would_cite":false,"duration_ms":39635,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["62F15","62M30","60G15"],"pacs":[],"model":"deepseek-v4-flash","headline":"This paper shows a two-layer deep Gaussian process can be upgraded to use and predict gradients via the chain rule, beating gradient-enhanced GPs and plain DGPs on nonstationary surfaces.","keywords":["deep Gaussian processes","gradient enhancement","surrogate modeling","computer experiments","nonstationary emulation","gradient predictions","Vecchia approximation","Bayesian MCMC"],"falsifier":"On a one-dimensional nonstationary example such as the step-like function in the paper, record det(∇_x W) at each training location across the MCMC iterations: if any determinant is zero or the linear solve in Eq. (13) is numerically unstable, gradient-enhancement is not computable there. Then run the sampler twice — once with the Eq. (16) acceptance ratio and once with a |det(∇_x W)| factor accounting for the change of variables — and compare the posterior predictive distributions of the response and its gradients; any visible divergence shows the omitted Jacobian term is not negligible.","tokens_in":18556,"feed_emoji":"📈","tokens_out":13740,"duration_ms":123968,"temperature":0.7,"pith_summary":"Deep Gaussian processes (DGPs) — two-layer compositions of Gaussian processes that warp a nonstationary input space into a smoother one — are popular surrogates for expensive computer experiments, but they had not been integrated with derivative information. This paper tries to establish that both gradient-related capabilities are achievable inside a fully Bayesian DGP: absorbing observed gradient data into training ('gradient enhancement') and returning posterior predictive distributions of gradients at unobserved inputs. The bridge is the multivariate chain rule, which factorizes the response's input-gradient into the latent warp's derivative matrix times the outer layer's gradient in the warped space, so the model can convert between observed input-gradients and latent-space gradients. The paper argues, with benchmarks on nonstationary test functions and, per the abstract, a real quantum-mechanics computer experiment simulating molecular energy and forces as a function of atomic position, that gradient-enhanced DGPs outperform gradient-enhanced GPs and non-enhanced DGPs, and that DGP gradient predictions are more effective than GP gradient predictions. A practitioner cares because many physics and engineering simulators already return gradients at little extra cost, so a surrogate that spends that information wisely delivers more accuracy per expensive simulation run.","feed_headline":"Gradient-aware deep GPs win on nonstationary surfaces","feed_subtitle":"A chain-rule upgrade lets deep GPs use derivative data and predict gradients for costly simulations.","key_machinery":"The load-bearing object is the multivariate chain rule written as the linear system ∇_x y = (∇_x W)(∇_w y) (Eq. 13), in which ∇_x W is the D×D matrix of partial derivatives of the latent warp W with respect to the original inputs and ∇_w y collects the outer layer's derivatives with respect to the warped coordinates. The identity does double duty: during training, observed input-gradients are pushed backwards through the warp to form pseudo-observations of the latent-space gradients that enter the outer-layer likelihood; at prediction, inferred latent-space gradients are pushed forwards to produce posterior predictive distributions of the input-gradients. Around this identity the paper build","core_discovery":"The central claim is that a deep Gaussian process can be made fully gradient-aware while staying within exact MCMC-based Bayesian inference, and that doing so pays off precisely on the nonstationary surfaces where DGPs are chosen in the first place. Because the derivatives of a Gaussian process are themselves GPs, the paper upgrades both layers of the two-layer DGP to carry derivative information: the inner layer samples each latent node jointly with its input-gradients, and the outer layer carries the response together with its gradients with respect to the warped inputs. The two are linked by the multivariate chain rule ∇_x y = (∇_x W)(∇_w y), which converts observed input-gradients into l","pith_inferences":["The chain-rule bridge is compositional, so the same construction should extend to DGPs with more than one hidden layer by linking the gradient identities layer by layer; the paper restricts itself to a single latent layer, but nothing in the mechanism appears to require that.","Because converting observed input-gradients into latent-space gradients requires solving Eq. (13) at every training location, training designs that avoid near-singular latent warps — or a regularized fallback for the solve — could prove important in applications with flat or degenerate regions.","A direct testable extension: run the sampler with and without a |det(∇_x W)| Jacobian factor in the Eq. (16) acceptance ratio on a strongly warping function; any divergence in the resulting posterior predictive distributions would indicate the change of variables is not Jacobian-free.","The gradient predictions the paper develops could feed gradient-based acquisition functions in Bayesian optimization, where the DGP's learned warp might locate sharp features that stationary gradient surrogates smooth over — an application the paper points to but does not implement."],"forward_implications":["On nonstationary response surfaces, gradient-enhanced DGPs give lower RMSE and better-calibrated uncertainty than gradient-enhanced GPs and non-enhanced DGPs, with the largest gains when training data are scarce.","DGP posterior predictions of gradients are more accurate than GP gradient predictions on nonstationary functions, which should strengthen downstream tasks that rely on gradients, such as derivative-based optimization and active learning.","Gradient information can be folded into a DGP without abandoning full Bayesian inference: the latent warp and its gradients are sampled jointly by MCMC, avoiding variational or moment-matching shortcuts.","The Vecchia approximation built into the open-source software keeps gradient-enriched models computationally feasible, extending the method to larger training sets than the cubic-cost baseline would allow."],"fun_headline_variants":["Deep GP gets a derivative boost for nonstationary surrogates","Use gradients to make deep GPs smarter on tricky surfaces","Gradient-enhanced deep GPs improve nonstationary emulators","Chain-rule upgrade lets deep GPs use gradient data","Deep GP uses gradients to predict and enhance surrogate models"],"cache_read_input_tokens":2304,"weakest_assumption_plain":"Everything rests on the chain-rule conversion being sound: observed input-gradients are turned into latent-space gradients by solving ∇_x y = (∇_x W)(∇_w y), which requires the latent warp's derivative matrix ∇_x W to be invertible at every training location, while the MCMC acceptance probability in Eq. (16) treats the change of variables as carrying no Jacobian term.","fun_headline_variants_meta":{"raw":{"variants":["Deep GP gets a derivative boost for nonstationary surrogates","Use gradients to make deep GPs smarter on tricky surfaces","Gradient-enhanced deep GPs improve nonstationary emulators","Chain-rule upgrade lets deep GPs use gradient data","Deep GP uses gradients to predict and enhance surrogate models"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000983,"raw_usage":{"total_tokens":4040,"prompt_tokens":810,"completion_tokens":3230,"prompt_tokens_details":{"cached_tokens":256},"prompt_cache_hit_tokens":256,"prompt_cache_miss_tokens":554,"completion_tokens_details":{"reasoning_tokens":3147}},"tokens_in":554,"tokens_out":3230,"duration_ms":24997,"temperature":1.0,"reasoning_tokens":3147,"cache_read_input_tokens":256,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-03T15:05:15.721798+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"On a one-dimensional nonstationary example such as the step-like function in the paper, record det(∇_x W) at each training location across the MCMC iterations: if any determinant is zero or the linear solve in Eq. (13) is numerically unstable, gradient-enhancement is not computable there. Then run the sampler twice — once with the Eq. (16) acceptance ratio and once with a |det(∇_x W)| factor accounting for the change of variables — and compare the posterior predictive distributions of the response and its gradients; any visible divergence shows the omitted Jacobian term is not negligible.","supporting_citations":[],"review_version":1}