{"id":"dcb5a89a-948d-46d2-b9cf-e32e80f51803","arxiv_id":"2507.10952","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":3,"one_line_summary":"Heteroskedastic rational kriging is a fast, data-driven variance extension of rational kriging that improves active learning for computer experiments.","lead":"This paper introduces heteroskedastic rational kriging, a fast Gaussian process variant whose variance changes across the input space. Used for active learning, it concentrates new simulations in regions where the function changes rapidly, matching or beating slower deep Gaussian process methods on several test problems.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The n-parameter variance weight vector c is fit with no penalty or consistency analysis; if it overfits residuals, the acquisition function in Eq. 17 misallocates points and the central active-learning claim is unsupported.","rationale":"The reader's weakest assumption and my load-bearing concern are essentially the same: the high-dimensional, unpenalized fit of c is the point where the active-learning advantage could fail. I found no internal mathematical inconsistency in the derivation of the marginal likelihood, profile objective, or acquisition function; the weak spot is statistical robustness rather than algebra. Because this concern is real and unaddressed, a conditional verdict is appropriate. The paper does provide meaningful empirical support on several test functions and two real datasets, and the speed advantage over DeepGP is concrete, so I would not reject it. However, the lack of any overfitting or consistency analysis for c, combined with no released code, keeps the central claim from being fully established. Thus the reader's CONDITIONAL verdict should stand unchanged.","tokens_in":12864,"tokens_out":5417,"duration_ms":75010,"concrete_test":"Use a null simulator with no heteroskedasticity: draw 100 replicates of n=100 points from a stationary GP (e.g., Matern 5/2 with fixed lengthscales), fit HRK via Algorithm 1 on each, and compute the empirical standard deviation of estimated tau(x)/nu over a dense grid. Then run 30 replicates of active learning with HRK and OK on a stationary test function. If the median grid-CV of tau exceeds a pre-specified threshold (say 0.2) and HRK's RMSE/IS does not beat OK's, or if HRK places extra points in regions with no true variance increase, then the c fit is overfitting and the central claim is not established. If tau stays nearly flat and HRK matches OK, the concern is refuted.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central claim requires the learned heteroskedastic variance nu^2/(c0+r(x)'c)^2 to reflect genuine structure, not noise. In Algorithm 1 Step 2, after fixing theta at theta_RK and mu at mu_RK, the method minimizes g(c) over n nonnegative weights with only c'c <= 1 (Eq. 16). This is an unpenalized maximum-likelihood fit of a variance function with as many parameters as data points. The statement that initializing at c_RK 'enforces some sort of regularization' is not statistical regularization; the optimizer is free to move far from c_RK, and the term -2/n * sum(log d_i) rewards large d_i without any complexity penalty. There is no consistency, contraction, or overfitting analysis for c, nor for the fixed-theta approximation. If the fit exploits residuals, tau(x) will vary randomly, and the acquisition in Eq. 17 will dedicate design points to spurious high-variance regions. That would break the stated advantage over stationary OK, which at least remains space-filling. The empirical sections show HRK often helps, but they do not isolate c stability, and no code is provided, so this failure mode is not checked.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes Heteroskedastic Rational Kriging (HRK), a Gaussian-process model in which the response variance is modeled as tau(x) = nu / (c0 + r(x)'c), and uses the resulting posterior variance as an active-learning acquisition criterion. The model is estimated by profile maximum likelihood with theta fixed at the rational-kriging estimate, and the acquisition function in Eq. (17) seeks points with high ordinary-kriging variance relative to the estimated local variance scale. The authors compare HRK with ordinary kriging and deep Gaussian processes on six test functions and two real datasets, reporting lower RMSE and interval scores in most settings and computation times orders of magnitude faster than DeepGP.","tokens_in":13197,"tokens_out":5921,"duration_ms":72661,"significance":"If the variance-parameter estimation is stable, HRK is a valuable contribution: it gives a tractable likelihood and closed-form gradient for a nonstationary GP, an efficient sequential-design criterion, and empirical evidence of improved accuracy over OK with a large speed advantage over DeepGP. The paper is also commendable for including two real-data studies and for making the algorithmic steps (Algorithms 1 and 2) explicit. However, no code or data are provided, and the central variance weights c are fit without a penalty, so the paper's main claim rests on empirical demonstrations that do not directly address overfitting of the variance function. The proposed method is novel and the qualitative conclusions are plausible, but the missing stability analysis is load-bearing for the active-learning claim.","major_comments":[{"comment":"The unpenalized maximum-likelihood fit of c, an (n+1)-dimensional vector estimated from n observations with only the constraints c'c <= 1 and c >= 0, is the load-bearing step for the heteroskedastic variance. The statement in Section 3 that initializing at c_RK 'enforces some sort of regularization' is not a statistical penalty; the optimizer is free to move far from the initial value, and the objective in Eq. (16) contains no complexity term. If c overfits the residuals, tau(x) and the acquisition function in Eq. (17) will concentrate design points in spurious high-variance regions, which would undermine the central active-learning claim. The empirical sections do not isolate this failure mode: there are no plots or summaries of the estimated c, no re-randomization stability checks, and no experiments with a known tau(x). Please add either a penalized or otherwise regularized estimator with a consistency or shrinkage analysis, or a diagnostic study showing that c recovers a known variance function and is stable under re-randomization of the initial design.","section":"Section 3, Eq. (16) and Algorithm 1"},{"comment":"Theta is fixed at the rational-kriging estimate when optimizing c, and mu is fixed at mu_RK in Eq. (16). Since R and R^{-1} enter both the likelihood and the acquisition criterion, a poor theta estimate can distort the variance function even if c were otherwise well estimated. The paper does not report any sensitivity analysis for this profile approximation. Please justify fixing theta at theta_RK, for example by showing that the active-learning results are stable to a small grid of theta values or to a one-step joint update, or by stating precisely why theta_RK is sufficient for the heteroskedasticity to be captured.","section":"Algorithm 1, Step 2 and Section 3"},{"comment":"The abstract claims 'comparable or better performance relative to other non-stationary Gaussian process-based methods', but in Figure 7 DeepGP is run only for a few steps for most test functions, with the full trajectory shown only for the Gramacy-Lee function. The summary curves therefore compare HRK and OK over the full horizon but DeepGP over a shorter one, and the current wording is stronger than what the truncated comparison supports. Please report the exact number of active-learning steps completed for each function, state how the medians are computed when trajectories have different lengths, and qualify the comparison accordingly.","section":"Section 5, Figure 7"}],"minor_comments":[{"comment":"There is a typo in 'ALM startegy' that should read 'strategy'.","section":"Section 4, paragraph after Eq. (17)"},{"comment":"The notation theta_RK, mu_RK, and c_RK is used before being defined; please define these as the rational-kriging estimates from Joseph (2024).","section":"Algorithm 1, Step 1"},{"comment":"It would help to state explicitly that the constant factor nu^2 is omitted from the argmax in Eq. (17), since the full posterior variance is s^2(x) = nu^2 (1 - r'R^{-1}r) / (c0 + r'c)^2.","section":"Section 4, Eq. (17)"},{"comment":"The shaded bands are described as 5th and 95th quantiles; please clarify whether these are across the 10 repetitions or across a different number of replications, and state the number of repetitions in each panel.","section":"Section 5, Figure 7"},{"comment":"The description of the molecular descriptors is too brief to reproduce: please give details on how the RDKit descriptors were computed, any preprocessing or scaling, and how the five principal components were obtained.","section":"Section 6.1"},{"comment":"No code or data links are provided for the simulations or real-data analyses. Given the novelty of the MMA-based optimization and the central role of the estimated c vector, releasing code would substantially strengthen reproducibility.","section":"Global"}],"recommendation":"major_revision","confidential_remarks":"The paper makes a plausible and useful contribution, but I am not comfortable accepting it without a direct treatment of the high-dimensional variance-parameter estimation. The lack of regularization or consistency analysis for c, combined with the absence of code, leaves the central active-learning claim more fragile than the empirical sections suggest. The DeepGP comparison should also be reported with explicit step counts. I see no grounds for rejection: the method is novel, computationally attractive, and the reported improvements are consistent with the proposed mechanism."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Colleague,\n\nThe short version: this is a useful extension of rational kriging that actually works in practice. The authors treat the RK weights c as variance parameters, estimate them by empirical Bayes, and use the resulting variance in an ALM-style acquisition. The math in Section 3 checks out, the gradient is explicit, and the algorithm is genuinely fast because R^{-1} is computed once. The empirical section is the best part: six test functions plus two real data studies, with HRK beating OK consistently and often matching or beating DeepGP at a fraction of the cost.\n\nThe main soft spot is the estimation of c. It's an n+1 parameter vector fit by MLE with only the constraint c'c <= 1. That's a lot of freedom, and the paper doesn't provide any regularization or consistency analysis. The claim that initializing at c_RK 'enforces some sort of regularization' is thin. In principle, c could overfit residual noise and produce a spurious variance surface, which would then misdirect the acquisition function. The empirical results suggest this doesn't happen dramatically in the tested cases, but the paper doesn't isolate the stability of c, and no code is provided to check. A dedicated experiment with a known heteroskedastic function where you compare c recovery under increasing n would settle this.\n\nOther minor issues: theta is fixed at the RK estimate rather than jointly estimated, which is a pragmatic but asymmetric treatment; the baselines are OK and DeepGP only, leaving out treed GP or other nonstationary methods that might be more comparable. The authors note DeepGP is slow, which justifies the choice for timing, but a couple of extra points for accuracy comparison would help.\n\nIf you work on emulation or sequential design, this is worth a read. It's a solid contribution, not a breakthrough, but it deserves proper peer review. My recommendation: send it to a journal like Technometrics or JASA, and ask the authors to provide code and address the c overfitting question head-on.\n\nBest.","headline":"A solid, fast heteroskedastic emulator with convincing empirical results, but the high-dimensional variance weight estimation needs a closer look.","tokens_in":13630,"tokens_out":2029,"would_cite":true,"duration_ms":25100,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["62K05","62M30","62G08"],"pacs":[],"model":"deepseek-v4-flash","headline":"A heteroskedastic version of rational kriging lets active learning place new design points in regions where the response surface actually varies, yielding surrogates as accurate as deep Gaussian process methods at a fraction of the cost.","keywords":["active learning","heteroskedastic Gaussian process","rational kriging","sequential design","computer experiments","emulation","nonstationarity","surrogate model"],"falsifier":"Run HRK active learning on a synthetic function whose variance is known to be constant, and measure how far the estimated $\\tau(x)$ deviates from flat and whether ALM still redistributes points away from uniformity; if the variance function becomes strongly nonconstant on stationary data, the $c$-estimation is overfitting. On a function with a known narrow high-variability band, the method should place a large fraction of new points inside that band.","tokens_in":12675,"feed_emoji":"🎯","tokens_out":6773,"duration_ms":69730,"temperature":0.7,"pith_summary":"The paper proposes heteroskedastic rational kriging (HRK), a Gaussian process model whose variance is allowed to vary across the input space by making the rational kriging weight vector part of the likelihood. In active learning, the next design point is chosen where HRK's posterior variance is largest, and that criterion balances space filling against an explicit estimate of where the response is hard to predict. The paper reports that this concentrates runs in high-variability regions and gives predictions comparable to or better than deep Gaussian process-based active learning, while running orders of magnitude faster. If the claim holds, active learning for expensive computer models can focus simulations where they matter without paying the MCMC cost of other nonstationary methods.","feed_headline":"Heteroskedastic kriging targets regions that truly vary","feed_subtitle":"HRK beats stationary kriging and deep Gaussian process surrogates on accuracy, at orders of magnitude lower cost.","key_machinery":"The load-bearing object is the heteroskedastic rational kriging model $y(x) = \\mu + \\tau(x)Z(x)$ with $\\tau(x) = \\nu/(c_0 + r(x)'c)$, where $Z$ is a stationary Gaussian process and $c_0, r(x), c$ come from the rational kriging predictor. The denominator $c_0 + r(x)'c$ modulates the local variance; fitting $c$ by maximizing the marginal likelihood (with $c'c \\le 1$, $c \\ge 0$) lets the model learn where the response varies. The active learning criterion is the resulting posterior variance $s^2(x) \\propto (1 - r(x)'R^{-1}r(x))/(c_0 + r(x)'c)^2$, whose numerator preserves space-filling behavior and whose denominator directs new points into high-variance regions. The whole algorithm is fast because the correlation matrix and its inverse are computed once before the $c$-optimization.","core_discovery":"The central claim is that rational kriging's variance multiplier $\\tau(x) = \\nu/(c_0 + r(x)'c)$ already contains enough freedom to model heteroskedasticity, and that estimating the weight vector $c$ jointly with the other hyperparameters turns active learning from a space-filling exercise into a targeted search. The paper derives an empirical-Bayes objective, fixes the lengthscales from an initial rational kriging fit, and optimizes the resulting $(n+1)$-dimensional problem with a gradient method so that $R$ and $R^{-1}$ are computed only once. In the acquisition step, HRK's posterior variance is used as the ALM criterion; the numerator keeps the design spread out while the denominator concentrates points where the modeled variance is high. Experiments on six test functions and two real datasets show comparable or better accuracy than deep Gaussian process surrogates, with a speed advantage that is orders of magnitude.","pith_inferences":["The same acquisition function could be used inside Bayesian optimization objectives, such as expected improvement, to balance exploration and exploitation; the paper only mentions this direction in its conclusion.","The $(n+1)$-dimensional $c$ estimate is not regularized beyond the unit-ball and nonnegativity constraints, so the learned variance surface may overfit small initial designs; a simulation with known heteroskedasticity could quantify this.","The candidate-set approximation via down-sampled Latin hypercubes may miss narrow high-variance ridges in high dimensions; replacing it with adaptive Voronoi-style candidates is a natural stress test.","The speed comparison applies to point-estimate fitting; the method's uncertainty statements do not propagate uncertainty in $c$, so coverage intervals may be optimistic when the initial sample is small."],"forward_implications":["With HRK, active learning MacKay (ALM) no longer degenerates into a maximin space-filling design; the method places more runs near sharp transitions, as shown for the underdamped oscillator and Gramacy-Lee function.","HRK improves both root mean-squared error and interval score over ordinary kriging in nearly every simulation and both real-data studies.","The method matches or exceeds deep Gaussian process active learning on most test functions while being orders of magnitude faster, making it practical for sequential design.","Because only point estimates of hyperparameters are needed, the entire active learning loop can be run repeatedly at modest cost."],"supporting_citations":[{"why":"proposes rational kriging, supplying the predictor form and the variance multiplier that HRK generalizes","marker":"Joseph (2024)"},{"why":"gives the entropy-based active learning criterion ALM that HRK's acquisition adapts","marker":"MacKay (1992)"},{"why":"shows maximum-entropy designs approach maximin designs, motivating the need for heteroskedasticity","marker":"Johnson et al. (1990)"},{"why":"provides the related maximum-likelihood estimation for iterated kernel regression that HRK's c-estimation extends","marker":"Kang and Joseph (2016)"},{"why":"supplies the MMA optimizer used in Algorithm 1 to minimize the c-objective","marker":"Svanberg (2002)"},{"why":"defines the deep Gaussian process active-learning baseline and algorithmic settings HRK is compared against","marker":"Sauer et al. (2023b)"},{"why":"provides the Gramacy-Lee test function and the LGBB benchmark design used in the simulations and the real study","marker":"Gramacy and Lee (2009)"},{"why":"gives the twinning down-sampling used to build the candidate set for the acquisition search","marker":"Vakayil and Joseph (2022)"}],"fun_headline_variants":["Heteroskedastic kriging targets where functions vary","Kriging that adapts to variation for active learning","Fast kriging focuses samples on changing regions","Active learning via rational kriging that spots variation","Kriging that learns where to look, fast"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise is that the fitted weight vector $c$, constrained only by $c'c \\le 1$ and $c \\ge 0$, identifies the true local variance from the initial sample instead of overfitting it.","fun_headline_variants_meta":{"raw":{"variants":["Heteroskedastic kriging targets where functions vary","Kriging that adapts to variation for active learning","Fast kriging focuses samples on changing regions","Active learning via rational kriging that spots variation","Kriging that learns where to look, fast"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000307,"raw_usage":{"total_tokens":1708,"prompt_tokens":850,"completion_tokens":858,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":466,"completion_tokens_details":{"reasoning_tokens":780}},"tokens_in":466,"tokens_out":858,"duration_ms":8764,"temperature":1.0,"reasoning_tokens":780,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-06T17:20:18.147344+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run HRK active learning on a synthetic function whose variance is known to be constant, and measure how far the estimated $\\tau(x)$ deviates from flat and whether ALM still redistributes points away from uniformity; if the variance function becomes strongly nonconstant on stationary data, the $c$-estimation is overfitting. On a function with a known narrow high-variability band, the method should place a large fraction of new points inside that band.","supporting_citations":[{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"gives the entropy-based active learning criterion ALM that HRK's acquisition adapts"},{"cited_title":"E., Moore, L","cited_arxiv_id":null,"evidence_quote":"shows maximum-entropy designs approach maximin designs, motivating the need for heteroskedasticity"},{"cited_title":"and Joseph, V","cited_arxiv_id":null,"evidence_quote":"provides the related maximum-likelihood estimation for iterated kernel regression that HRK's c-estimation extends"},{"cited_title":"(2002), A class of globally convergent optimization methods based on conservative convex separable approximations, SIAM journal on optimization, 12, 555--573","cited_arxiv_id":null,"evidence_quote":"supplies the MMA optimizer used in Algorithm 1 to minimize the c-objective"},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"provides the Gramacy-Lee test function and the LGBB benchmark design used in the simulations and the real study"},{"cited_title":"and Joseph, V","cited_arxiv_id":null,"evidence_quote":"gives the twinning down-sampling used to build the candidate set for the acquisition search"}],"review_version":1}