Why and When Neural Networks Improve Local Approximation in Optimization

# Why and When Neural Networks Improve Local Approximation in Optimization
Source: [https://arxiv.org/html/2608.24963](https://arxiv.org/html/2608.24963)
Chengkuo Bian††thanks:University of California, Berkeley, CA 94720, USA\. Email:chengkuobian@gmail\.com\.Pengcheng Xie††thanks:Corresponding author\. Applied Mathematics and Computational Research Division, Lawrence Berkeley National Laboratory, University of California, 1 Cyclotron Road, Berkeley, CA 94720, USA\. Email:pxie@lbl\.gov,pxie98@gmail\.com\.

August 2026

###### Abstract

Published experience with neural surrogates in derivative\-free optimisation is contradictory: the same family of models that cuts the evaluation count of one solver leaves another unchanged, or makes it worse\. We show that the contradiction dissolves once three factors are stated, and that these, rather than the fit accuracy a training curve reports, are what delimit when a learned local model pays\.*Role*: a surrogate that proposes candidates the true objective must still approve helps, while one that replaces a gradient the solver depends on hurts\.*Radius*: a model fitted to an optimisation path is reliable only inside a bounded neighbourhood, and its error neither vanishes as that neighbourhood shrinks nor survives its growth\.*Room*: a surrogate can only accelerate progress the base method is still able to make\. We formalise radius\-aware local generalisation, relate it to the classical fully linear condition, and test each factor with the surrogate class, training pipeline and base method held fixed\. Over117117benchmark instances safeguarded assistance raises the instances solved to high accuracy from6767to8484while gradient replacement lowers them to6565; removing the gradient term from the training loss cuts surrogate acceptance from0\.7030\.703to0\.1480\.148; and10001000paired comparisons over ten noise levels show no noise threshold, only a base method that stops early\. The same factors bound the gain: a model\-based trust\-region solver, which leaves little room, drops from8888to8686when the identical surrogate is attached, and released interpolation software stays ahead at103103, and on a Monte\-Carlo inventory model repairing the acceptance interface is worth10\.4010\.40cost units against0\.000\.00for the surrogate\.

Keywords:simulation optimization; derivative\-free optimization; surrogate metamodels; neural networks; Sobolev training; stochastic oracles; local generalization; data profiles\.

## 1Introduction

We consider the unconstrained minimization problem

minx∈ℝn⁡f⁡\(x\),\\min\_\{x\\in\\mathbb\{R\}^\{n\}\}f\(x\),\(1\)where each value offfcomes from a simulation run rather than from a formula, and no derivative is available\. This is the standard situation in simulation\-based optimisation: a queueing, inventory, agent\-based or physics model is executed at a candidate designxx, and the resulting performance measure is all the optimiser gets to see\. Runs are slow, so the currency of the problem is the number of evaluations\. Derivative\-free optimisation \(DFO\) is the branch of the field built for this currency; the monograph of Conn et al\.\[[1](https://arxiv.org/html/2608.24963#bib.bib1)\]and the surveys\[[2](https://arxiv.org/html/2608.24963#bib.bib2),[3](https://arxiv.org/html/2608.24963#bib.bib3),[5](https://arxiv.org/html/2608.24963#bib.bib5)\]cover it\.

An optimisation run leaves behind a record: the points visited, their simulated values, the step sizes that worked, and any finite\-difference gradients that were paid for along the way\. Fitting a cheap metamodel to that record and using it to guide the search is an old idea in the simulation literature, and neural networks are a natural modern choice of metamodel\. We should say at the outset where our answer lands, because the title asks two questions and they have different answers\.*Why*a learned model improves local approximation is the easier one: fitting every evaluation the run has already paid for averages away noise that differencing amplifies, and gradient information in the training loss restrains the curvature a value\-only fit is free to invent\. Both effects are measurable: removing the gradient term from the loss cuts the rate at which surrogate steps survive validation by a factor of nearly five\.*When*that better approximation reaches the optimisation is the harder question\. Approximation quality is a precondition and not a sufficient one; what delimits whether it converts into fewer evaluations are three factors we identify empirically\. Two are properties of the surrounding algorithm: the surrogate’s role, and what the base method leaves, separated below into unused headroom and whether its acceptance test still functions\. The third, the radius at which the method works, is a property of the fit read at the scale the solver operates on, so it is not independent of approximation quality but is not what a training curve reports either\. Where any one fails, a better local model does not become a better optimiser, most sharply on the Monte\-Carlo model of Section[9\.6](https://arxiv.org/html/2608.24963#S9.SS6)\. Locating that boundary is the result of this paper rather than a caveat attached to it\.

Among DFO methods, model\-based trust\-region algorithms constitute one of the most successful frameworks\. These methods construct local surrogate models and compute trial steps by approximately minimizing the models within trust regions\. Classical polynomial interpolation models and their theoretical foundations are well established\[[4](https://arxiv.org/html/2608.24963#bib.bib4),[1](https://arxiv.org/html/2608.24963#bib.bib1),[5](https://arxiv.org/html/2608.24963#bib.bib5)\]\. Underdetermined quadratic interpolation, pioneered in Powell’s NEWUOA\[[6](https://arxiv.org/html/2608.24963#bib.bib6)\]and refined through later model\-updating strategies\[[7](https://arxiv.org/html/2608.24963#bib.bib7)\], forms the backbone of practical solvers\. Finite\-difference\-based methods offer a complementary line whose practical performance has been reassessed in recent studies\[[8](https://arxiv.org/html/2608.24963#bib.bib8),[9](https://arxiv.org/html/2608.24963#bib.bib9)\]\.

Despite these advances, classical model\-based approaches face fundamental limitations in high\-dimensional settings\. Constructing and maintaining polynomial models becomes increasingly expensive due to the growth in the number of parameters, and the quality of interpolation depends critically on the geometry of sample sets\. In practice, limited and potentially noisy data lead to fragile models and unreliable steps\. These challenges are closely related to the curse of dimensionality and motivate the exploration of richer approximation classes and learning\-based surrogates\.

Neural networks \(NNs\) provide a flexible alternative to classical surrogate models\. With smooth activation functions, NNs can approximate both function values and gradients while implicitly performing noise smoothing and feature extraction\. Recent studies, however, reveal a nuanced picture\. Giovannelli et al\. show that although NN surrogates can achieve strong approximation accuracy for function values and gradients, they may underperform in capturing second\-order information and do not necessarily improve optimization performance when directly replacing finite\-difference gradients in quasi\-Newton updates\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]\. In contrast, Taminiau et al\. demonstrate that NNs can significantly improve performance when used in a safeguarded manner, where surrogate\-based steps are accepted only if the true objective decreases sufficiently\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]; related learning\-augmented designs include the neural\-network\-accelerated implicit filtering of\[[13](https://arxiv.org/html/2608.24963#bib.bib13)\]\. These contrasting findings suggest that the effectiveness of NN surrogates depends on how the surrogate is embedded in the algorithm, not on approximation accuracy alone\.

From a theoretical perspective, the role of surrogate models in trust\-region methods is governed by uniform approximation properties, such as fully linear and fully quadratic conditions\[[1](https://arxiv.org/html/2608.24963#bib.bib1)\]\. While polynomial models can satisfy these conditions under suitable sampling geometry, it remains less understood when learned surrogates—particularly neural networks—can achieve comparable guarantees\. Recent advances in neural approximation theory provide partial insights\. For example, objective\-value\-change and shape\-based perspectives offer new ways to characterize local approximation quality beyond classical interpolation error\[[14](https://arxiv.org/html/2608.24963#bib.bib14)\]\. Additionally, the relationship between interpolation geometry and robustness, includingΛ\\Lambda\-poisedness and data irregularity, has been explored in\[[15](https://arxiv.org/html/2608.24963#bib.bib15)\]\. These developments indicate that learning\-based surrogates should be analyzed through a combination of approximation theory, data geometry, and algorithmic integration\.

Another important dimension is scalability and application context\. Modern optimization problems often involve large\-scale systems, distributed settings, and black\-box pipelines, and applications ranging from simulation\-based design to privacy\-constrained and uncertainty\-aware settings impose their own constraints on data access, noise, and computational cost\[[2](https://arxiv.org/html/2608.24963#bib.bib2)\]\. These trends further motivate the integration of learning\-based surrogates into DFO frameworks in a principled and scalable manner\.

##### Main message\.

This paper is organised around one claim:

> *Whether an NN surrogate improves a DFO method is determined primarily by \(i\) thealgorithmic rolein which the surrogate is embedded, replacement of a core quantity of the base method versus safeguarded assistance, and \(ii\) whether the method operates within the surrogate’sradius of reliable local generalization, with a third factor emerging from the experiments: what the base method leaves undone\. Pointwise approximation accuracy alone is a poor predictor of optimization benefit\.*

The claim explains the disagreement described above\. Under replacement, an error in the surrogate goes straight into a quantity the solver acts on, and from there into its search direction and its curvature estimates\. Under assistance, the same error is filtered by a test on the true objective, which the surrogate cannot fake\. The radius half of the claim is equally concrete: a model fitted to a handful of points collected along an optimisation path is accurate over a bounded region, and a solver working at a larger scale than that region will be misled no matter how well the model fits its training data\.

The contributions are these\.

1. 1\.Three conditions, each shown necessary\.A learned local model improves the optimisation only when its role is assistance rather than replacement, the working radius lies inside the region where it generalises, and the base method has budget left\. Two we fail on purpose and watch the benefit go: replacement drops below the base method it was meant to help,6565against6767; a base method with no room left gains nothing\. The third we establish at the level of approximation\.
2. 2\.A measured effect and a measured boundary\.Over117117instances safeguarded assistance lifts its base method from6767to8484instances solved to high accuracy, a sixth of the benchmark, and the margin grows as the tolerance tightens\. Released Py\-BOBYQA solves103103and ends two orders of magnitude lower at the median, so the effect is real and does not reach the state of the art\. Extended Wood is the exception, and we report it as an existence result rather than a rule\.
3. 3\.ARAS, with a convergence theorem tied to a computable gate\.The blended direction is provably descent and the method provably inherits the base method’s𝒪⁡\(n​ϵ−2\)\\mathcal\{O\}\(n\\epsilon^\{\-2\}\)complexity once the blending weight is bounded by\(1−θ\)/\(1\+κk\)\(1\-\\theta\)/\(1\+\\kappa\_\{k\}\), the gate statistic the framework already computes\. Running it shows the blended step costs more than it returns, which is what the calibration says from the other side, so the surviving instance of the framework is the one we benchmark\.
4. 4\.Radius\-aware local generalisation, with elementary uniform error bounds connecting sample coverage and training error to the fully linear condition\.
5. 5\.Three self\-corrections reported rather than buried:a crossover atn=64n=64that was a defect in our own interpolation code, two geometric diagnostics we proposed and refuted, and a network setting used throughout an earlier version that proves undertrained, worth nine instances of117117\.
6. 6\.A warning about reference sets\.Scoring the surrogate variants against each other puts assistance ahead of a model\-based method; adding a released solver reverses the ordering\. A Moré–Wild comparison that omits released software measures its own pool\.

One caveat belongs in front of all of this rather than buried in the experimental section\. This is a mechanistic study, not a solver paper\. Every variant we compare shares one network class, one training pipeline and one base method, precisely so that a difference between two of them can be attributed to the mechanism we varied and not to an implementation detail; the price of that design is that none of the variants is tuned to be competitive\. That design does not, however, excuse us from asking whether the surrogate’s benefit merely repairs a weak base method, so Section[9\.2](https://arxiv.org/html/2608.24963#S9.SS2)adds a model\-based trust\-region solver of thenewuoatype\[[6](https://arxiv.org/html/2608.24963#bib.bib6),[1](https://arxiv.org/html/2608.24963#bib.bib1)\]and answers the question directly\. It largely does: the strong solver alone outperforms assistance on the finite\-difference method, and attaching the same surrogate to the strong solver adds nothing\. We report this because it sharpens the claim rather than weakening it\. What a surrogate supplies is progress the base method cannot make on its own, so its value is measured by the gap it closes, not by the model’s accuracy\.

Two scope decisions follow\. Evaluations, not seconds, are the currency: one evaluation is a model run in the setting these experiments are about, so the wall\-clock figures we report come from an unoptimised pure\-NumPy implementation and locate a break\-even oracle cost rather than rank algorithms\. And the dimensions are small on purpose: isolating one mechanism at a time means running every variant on every instance, which atn≤16n\\leq 16we can do exhaustively and atn=128n=128once\.

Sections[2](https://arxiv.org/html/2608.24963#S2)–[6](https://arxiv.org/html/2608.24963#S6)reread the two source studies through the role distinction, Section[7](https://arxiv.org/html/2608.24963#S7)makes the radius precise, Section[8](https://arxiv.org/html/2608.24963#S8)gives the embeddings and ARAS, Section[9](https://arxiv.org/html/2608.24963#S9)the experiments, and Sections[10](https://arxiv.org/html/2608.24963#S10)–[11](https://arxiv.org/html/2608.24963#S11)what follows from them\.

## 2Background and Problem Setting

We consider the unconstrained black\-box optimization problem

minx∈ℝn⁡f⁡\(x\),\\min\_\{x\\in\\mathbb\{R\}^\{n\}\}f\(x\),\(2\)where the objective function is expensive to evaluate and derivative information is unavailable, unreliable, or too costly to compute\. In this regime, the central question is not only how to build a local approximation offf, but also how that approximation is embedded into the optimization loop\. This distinction is especially important for neural\-network surrogates: a model that fits local data well may still fail to improve the optimization method if it is used in a fragile way\.

### 2\.1Model\-Based Optimization

A standard approach in derivative\-free optimization is to build, at iterationkk, a local modelmkm\_\{k\}offfaround the current iteratexkx\_\{k\}and to use this model to compute a trial step\. In trust\-region methods, one typically solves a subproblem of the form

min‖s‖≤Δk⁡mk​\(xk\+s\),\\min\_\{\\\|s\\\|\\leq\\Delta\_\{k\}\}m\_\{k\}\(x\_\{k\}\+s\),\(3\)whereΔk\>0\\Delta\_\{k\}\>0is the trust\-region radius\. The trial point is then accepted or rejected by comparing the actual decrease inffwith the decrease predicted bymkm\_\{k\}\. This framework separates two issues:*model quality*and*acceptance control*\. The model should be accurate onℬ⁡\(xk,Δk\)\\mathcal\{B\}\\\!\\left\(x\_\{k\};\\Delta\_\{k\}\\right\), while the acceptance test prevents the algorithm from over\-trusting the model when the approximation is poor\.

In classical model\-based DFO,mkm\_\{k\}is often chosen from polynomial or radial\-basis\-function families\. These models are attractive because they admit explicit interpolation or regression constructions and can be analyzed through sampling geometry and trust\-region theory\. In finite\-difference\-based methods, the role of a model is different but closely related: the method still constructs local information from function evaluations, often through approximate gradients, and then applies a first\-order update or a quasi\-Newton step\. The standard forward\-difference estimator is

\[gh\(x\)\]i:=f⁡\(x\+h​ei\)−f⁡\(x\)h,i=1,…,n,\[g\_\{h\}\(x\)\]\_\{i\}:=\\frac\{f\(x\+he\_\{i\}\)\-f\(x\)\}\{h\},\\qquad i=1,\\dots,n,\(4\)which costsn\+1n\+1evaluations per gradient and whose stepsizehhmust balance truncation error against noise amplification\.

For the purposes of this paper, it is useful to distinguish two algorithmic roles for a learned surrogate\.

- •Replacement:the surrogate directly substitutes a core quantity used by the base method, such as a gradient, a search direction, or a local model minimized by the algorithm\.
- •Assistance:the surrogate is used only to propose additional candidate points or cheap exploratory steps, while the true objective remains responsible for final acceptance\.

This distinction clarifies the contrast between the two primary reference papers\. Giovannelli et al\. study a replacement mechanism, where the surrogate gradient is injected into a finite\-difference BFGS\-type step\. Taminiau et al\. study an assistance mechanism, where a surrogate trained from the accumulated dataset is optimized cheaply to generate extra candidates, but these candidates are retained only when the true objective decreases sufficiently\. The inferential status of this reading should be stated plainly, since the rest of the paper rests on it\. The two studies differ in their base method, their architecture, their training protocol, their data policy and their test set as well as in role, and the supplementary document tabulates those differences; a pair of observational studies differing in six ways cannot establish which one matters\. What the pair supplies is a hypothesis\. What Section[9](https://arxiv.org/html/2608.24963#S9)supplies is a test of it, holding the network class, the training pipeline and the base method fixed and varying role alone\.

The data available to a learned surrogate typically come from the optimization trajectory itself\. If

Fk=\{\(yi,f⁡\(yi\)\)\}i=1NkF\_\{k\}=\\\{\(y\_\{i\},f\(y\_\{i\}\)\)\\\}\_\{i=1\}^\{N\_\{k\}\}collects previously evaluated points and function values, then a value\-based surrogatemθm\_\{\\theta\}may be trained by minimizing an empirical loss overFkF\_\{k\}\. When approximate gradients are also available, one can augment the dataset with

Gk=\{\(zj,g⁡\(zj\)\)\}j=1Mk,g⁡\(zj\)≈∇f​\(zj\),G\_\{k\}=\\\{\(z\_\{j\},g\(z\_\{j\}\)\)\\\}\_\{j=1\}^\{M\_\{k\}\},\\qquad g\(z\_\{j\}\)\\approx\\nabla f\(z\_\{j\}\),and train the surrogate through a Sobolev\-type objective of the form

LFk,Gk​\(θ\):=1Nk​∑i=1Nk\(mθ​\(yi\)−f⁡\(yi\)\)2\+1Mk​∑j=1Mk‖∇mθ​\(zj\)−g⁡\(zj\)‖2\+λ​‖θ‖2\.L\_\{F\_\{k\},G\_\{k\}\}\(\\theta\):=\\frac\{1\}\{N\_\{k\}\}\\sum\_\{i=1\}^\{N\_\{k\}\}\\bigl\(m\_\{\\theta\}\(y\_\{i\}\)\-f\(y\_\{i\}\)\\bigr\)^\{2\}\+\\frac\{1\}\{M\_\{k\}\}\\sum\_\{j=1\}^\{M\_\{k\}\}\\\|\\nabla m\_\{\\theta\}\(z\_\{j\}\)\-g\(z\_\{j\}\)\\\|^\{2\}\+\\lambda\\\|\\theta\\\|^\{2\}\.\(5\)This formulation is important for our discussion because it makes the surrogate both a value approximator and a device that shapes gradients\. In particular, once a differentiable modelmθm\_\{\\theta\}has been trained, one can obtain cheap gradient\-based steps from∇mθ\\nabla m\_\{\\theta\}without incurring additional evaluations offf\.

### 2\.2Fully Linear and Fully Quadratic Models

The standard language for local model quality in trust\-region DFO is given by the notions of fully linear and fully quadratic models\. These are*uniform*approximation properties on a neighborhood, not merely pointwise fitting conditions\.

###### Definition 1\(Fully linear model\)\.

A modelmkm\_\{k\}is said to be*fully linear*onℬ⁡\(xk,Δk\)\\mathcal\{B\}\\\!\\left\(x\_\{k\};\\Delta\_\{k\}\\right\)if there exist constantsκf,κg\>0\\kappa\_\{f\},\\kappa\_\{g\}\>0, independent ofkk, such that for allx∈ℬ⁡\(xk,Δk\)x\\in\\mathcal\{B\}\\\!\\left\(x\_\{k\};\\Delta\_\{k\}\\right\),

\|f⁡\(x\)−mk​\(x\)\|≤κf​Δk2,\|f\(x\)\-m\_\{k\}\(x\)\|\\leq\\kappa\_\{f\}\\Delta\_\{k\}^\{2\},\(6\)and

‖∇f​\(x\)−∇mk​\(x\)‖≤κg​Δk\.\\\|\\nabla f\(x\)\-\\nabla m\_\{k\}\(x\)\\\|\\leq\\kappa\_\{g\}\\Delta\_\{k\}\.\(7\)

###### Definition 2\(Fully quadratic model\)\.

Assume thatffis twice continuously differentiable onℬ⁡\(xk,Δk\)\\mathcal\{B\}\\\!\\left\(x\_\{k\};\\Delta\_\{k\}\\right\)\. A modelmkm\_\{k\}is said to be*fully quadratic*onℬ⁡\(xk,Δk\)\\mathcal\{B\}\\\!\\left\(x\_\{k\};\\Delta\_\{k\}\\right\)if there exist constantsκf,κg,κH\>0\\kappa\_\{f\},\\kappa\_\{g\},\\kappa\_\{H\}\>0, independent ofkk, such that for allx∈ℬ⁡\(xk,Δk\)x\\in\\mathcal\{B\}\\\!\\left\(x\_\{k\};\\Delta\_\{k\}\\right\),

\|f⁡\(x\)−mk​\(x\)\|≤κf​Δk3,\|f\(x\)\-m\_\{k\}\(x\)\|\\leq\\kappa\_\{f\}\\Delta\_\{k\}^\{3\},\(8\)‖∇f​\(x\)−∇mk​\(x\)‖≤κg​Δk2,\\\|\\nabla f\(x\)\-\\nabla m\_\{k\}\(x\)\\\|\\leq\\kappa\_\{g\}\\Delta\_\{k\}^\{2\},\(9\)and

‖∇2f​\(x\)−∇2mk​\(x\)‖≤κH​Δk\.\\\|\\nabla^\{2\}f\(x\)\-\\nabla^\{2\}m\_\{k\}\(x\)\\\|\\leq\\kappa\_\{H\}\\Delta\_\{k\}\.\(10\)

These conditions explain why trust\-region methods can attach rigorous meaning to the phrase “good local model\.” They guarantee that minimizingmkm\_\{k\}inside the trust region gives useful information about the behavior offfin that same region\. They also show why local approximation in optimization is stronger than ordinary regression accuracy: what matters is not just that the surrogate fits sampled values, but that it does so uniformly over the region and with the right derivative accuracy\.

This viewpoint is particularly helpful for understanding learned surrogates\. In a replacement strategy, the method acts directly on surrogate derivatives or on steps computed from the surrogate, so errors in gradients or curvature information can immediately affect search directions, quasi\-Newton updates, and globalization logic\. In that case, fully linear or fully quadratic behavior is the natural benchmark\. In an assistance strategy, by contrast, the surrogate is filtered through an acceptance test based on the true objective, so the algorithm can benefit from a useful but imperfect model without requiring the same degree of uniform reliability at every iteration\. This is one of the main conceptual reasons why assistance can be substantially more robust than replacement\.

### 2\.3Limitations of Classical Polynomial Models

Classical polynomial models remain the theoretical backbone of model\-based DFO, but they face several structural limitations in the settings that motivate learned surrogates\.

First, their parameter count grows rapidly with dimension\. A full quadratic model innnvariables has

\(n\+1\)​\(n\+2\)2\\frac\{\(n\+1\)\(n\+2\)\}\{2\}coefficients, so building and updating such a model becomes increasingly expensive asnngrows\. This is not only a storage issue; it also affects the number of sample points required for stable interpolation or regression and the cost of the associated linear algebra\.

Second, polynomial models are highly sensitive to sampling geometry\. Their quality depends on the sample set being sufficiently well poised\. When the available points are clustered, nearly collinear, or concentrated on a low\-dimensional subset of the trust region, the resulting interpolation or regression system can become ill\-conditioned\. In practical algorithms, this situation occurs naturally because data are generated by the optimization path rather than by a carefully designed experimental design\. For example, points produced by backtracking line search often lie approximately on a single ray, which yields poor coverage of orthogonal directions\.

Third, the data available in modern black\-box optimization are often noisy, heterogeneous, and reused across iterations\. In such cases, a rigid polynomial basis may fail to exploit repeated structure effectively\. Neural networks, in contrast, can absorb oversampled or irregular datasets more flexibly, and once trained they provide smooth, cheap\-to\-evaluate gradients\. However, this flexibility should not be confused with automatic algorithmic improvement\. Giovannelli et al\. show that even when neural surrogates produce competitive value and gradient approximations, directly replacing the finite\-difference gradient in a sensitive quasi\-Newton framework does not reliably improve optimization performance\. This suggests that the main difficulty is not merely approximation in isolation, but the interaction between approximation error, data geometry, and the update mechanism of the base solver\.

These observations motivate the perspective adopted in the rest of the paper\. We do not ask only whether a neural network can approximate a local objective well\. Instead, we ask a more optimization\-relevant question:*when does a learned local model provide the type of information that the surrounding algorithm can safely exploit?*The answer depends jointly on local coverage, training objectives, smoothness of the surrogate, and crucially, whether the surrogate is used to replace a core step or merely to assist the search\.

## 3Neural Networks for Local Approximation

### 3\.1Motivation and Advantages

Two mechanisms explain why NN surrogates can improve local approximation in optimization\. In terms of the central claim of Section[1](https://arxiv.org/html/2608.24963#S1), these mechanisms explain why NNs*can*help; whether they*do*help is decided by the embedding role and the region radius studied in Sections 7–9\.

##### Smoothing and reduced oscillations\.

Finite\-difference gradients exhibit a fundamental bias–variance trade\-off\. Under smoothness, forward\-difference truncation error isO⁡\(h\)O\(h\); under noise, variance scales asO⁡\(h−2\)O\(h^\{\-2\}\)\. A basic model clarifies this\. Suppose we observef~​\(x\)=f​\(x\)\+ξ​\(x\)\\tilde\{f\}\(x\)=f\(x\)\+\\xi\(x\)where𝔼⁡\[ξ\]=0\\mathbb\{E\}\[\\xi\]=0andVar⁡\(ξ\)=σ2\\mathrm\{Var\}\(\\xi\)=\\sigma^\{2\}\. The forward\-difference estimator obeys \(approximately\)

Var⁡\(\[g~h​\(x\)\]i\)=Var⁡\(ξ⁡\(x\+h​ei\)−ξ⁡\(x\)h\)≈2​σ2h2\.\\mathrm\{Var\}\\big\(\[\\tilde\{g\}\_\{h\}\(x\)\]\_\{i\}\\big\)=\\mathrm\{Var\}\\left\(\\frac\{\\xi\(x\+he\_\{i\}\)\-\\xi\(x\)\}\{h\}\\right\)\\approx\\frac\{2\\sigma^\{2\}\}\{h^\{2\}\}\.\(11\)Thus reducing truncation error by using smallerhhincreases noise variance dramatically\. Learned surrogates can smooth because they fit a global \(or at least region\-wide\) function class to many data points, implicitly averaging noise and producing a more stable gradient field\.

##### Trend prediction and directional guidance\.

Beyond smoothing, the surrogate can infer a coherent descent direction \(trend\) from scattered samples\. This is especially plausible when the training objective includes derivative targets: Sobolev training matches both values and gradients, improving the surrogate’s local directional fidelity\[[16](https://arxiv.org/html/2608.24963#bib.bib16),[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\. In Taminiau et al\. the surrogate is designed to be continuously differentiable and is used via gradient descent with an Armijo\-like condition, directly exploiting this directional information\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\.

##### Why these mechanisms are not automatic\.

Giovannelli et al\. emphasize a key caveat: high\-quality approximation \(especially in function values or gradients\) might not improve optimization if the surrogate is inserted into a sensitive part of the algorithmic loop\. In particular, they report that using surrogate gradients \(including NNs\) inside a practical FD\-BFGS step within FLE \(the full\-low evaluation framework of\[[9](https://arxiv.org/html/2608.24963#bib.bib9)\]\) does not reliably improve performance across CUTEst problems\[[17](https://arxiv.org/html/2608.24963#bib.bib17),[10](https://arxiv.org/html/2608.24963#bib.bib10)\]\. This motivates careful embedding and safeguards\.

### 3\.2Approximation Properties and Smoothness Conditions

In optimization, we usually need surrogate derivatives\. This leads to explicit smoothness constraints\.

###### Assumption 1\(Differentiable surrogate family\)\.

The surrogatemθ:ℝn→ℝm\_\{\\theta\}:\\mathbb\{R\}^\{n\}\\to\\mathbb\{R\}is continuously differentiable for allθ\\thetaproduced by training\.

Taminiau et al\. exclude ReLU specifically because it does not yield a differentiable model in their setting and focus on smooth activations\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\. Giovannelli et al\. compare ReLU, ELU, SiLU, Sigmoid, and Tanh, and identify SiLU as the best\-performing activation in their approximation benchmarks\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]\.

##### First\-order vs second\-order approximation\.

Giovannelli et al\. show that interpolation/regression models often provide better accuracy for second\-order approximations, while NNs can be competitive for zero\- and first\-order approximations at high training cost\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]\. This is important because trust\-region convergence theory relies on fully linear/quadratic properties \(Section[2\.2](https://arxiv.org/html/2608.24963#S2.SS2)\)\. If an algorithm requires accurate Hessians \(fully quadratic\), NN surrogates may require additional structure \(e\.g\. special architectures, regularization targeting second derivatives, or explicit curvature information in training\); otherwise, they may be better suited as first\-order tools within safeguarded frameworks\.

### 3\.3Comparison with Polynomial Models

Giovannelli et al\. interpret activation functions as nonlinear feature maps and propose activation\-enriched polynomial bases\. Two representative constructions \(written here in a compact “feature list” form\) replace quadratic cross terms with activated nonlinear features\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]:

ϕ~​\(x\)\\displaystyle\\tilde\{\\phi\}\(x\)=\{1,x1,…,xn,x12/2,s⁡\(x1​x2\),…,s⁡\(xn−1​xn\),xn2/2\},\\displaystyle=\\Big\\\{1,\\ x\_\{1\},\\dots,x\_\{n\},\\ x\_\{1\}^\{2\}/2,\\ s\(x\_\{1\}x\_\{2\}\),\\dots,s\(x\_\{n\-1\}x\_\{n\}\),\\ x\_\{n\}^\{2\}/2\\Big\\\},\(12\)ϕ^​\(x\)\\displaystyle\\hat\{\\phi\}\(x\)=\{1,x1,…,xn,s⁡\(x1\),…,s⁡\(xn\),x12/2,…,xn2/2\},\\displaystyle=\\Big\\\{1,\\ x\_\{1\},\\dots,x\_\{n\},\\ s\(x\_\{1\}\),\\dots,s\(x\_\{n\}\),\\ x\_\{1\}^\{2\}/2,\\dots,x\_\{n\}^\{2\}/2\\Big\\\},\(13\)wheres⁡\(⋅\)s\(\\cdot\)is an activation function\. Their findings suggest that such enrichment can “waive the necessity” of including cross terms in some cases, yielding fewer parameters while maintaining approximation quality\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]\. From a DFO viewpoint, this is a valuable insight: ML\-inspired feature design can create models that sit between classical polynomials and full NNs, potentially balancing sample efficiency and stability\.

Taminiau et al\. also compare NNs with RBF surrogates trained by the same Sobolev objective \(gradient term for RBF reduces to least squares\); they setλ=0\\lambda=0for RBF and solve by a direct least\-squares solver, choosing the minimum\-norm solution when non\-unique\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\. This highlights that “ML helps” need not mean “NNs are best”: the key is the training objective and the embedding into a safeguarded loop\.

## 4Training Data Selection for Local Models

### 4\.1Sampling Strategies in the Trust Region

The surrogate model is only as good as its data\. In model\-based DFO, a common paradigm is to sample insideℬ⁡\(xk,Δk\)\\mathcal\{B\}\\\!\\left\(x\_\{k\};\\Delta\_\{k\}\\right\)in a way that ensures geometry and stability\[[1](https://arxiv.org/html/2608.24963#bib.bib1),[4](https://arxiv.org/html/2608.24963#bib.bib4)\]\. Giovannelli et al\. sample training and testing datasets uniformly in a ballℬ⁡\(x0,1\)\\mathcal\{B\}\\\!\\left\(x\_\{0\};1\\right\)around a CUTEst\-provided initial pointx0x\_\{0\}\(for approximation evaluation\)\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]\. They also normalize data by shifting by−x0\-x\_\{0\}and scaling byΔ=maxi⁡‖xi−x0‖\\Delta=\\max\_\{i\}\\left\\lVert x\_\{i\}\-x\_\{0\}\\right\\rVertso the data lie inℬ⁡\(0,1\)\\mathcal\{B\}\\\!\\left\(0;1\\right\)\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]\.

Taminiau et al\. derive data from the base method itself: each outer iteration adds a function\-value point toFFand adds a finite\-difference gradient target toGGat the current iterate\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\. This is a typical “DFO\-as\-data\-generator” loop: the algorithm’s exploration naturally populates the region of interest, and the learned surrogate is updated online\.

##### Practical sampling variants \(recommended\)\.

In addition to uniform sampling, we recommend considering: \(i\)*diversity\-aware sampling*\(maximize minimum distances\) to improve coverage, and \(ii\)*poisedness\-improving repairs*for polynomial/RBF surrogates\[[1](https://arxiv.org/html/2608.24963#bib.bib1)\]\. When using NNs, diversity\-aware sampling can reduce redundancy and improve generalization within a trust region\.

### 4\.2Geometry and Poisedness of Sample Sets

For polynomial interpolation, sample geometry is formalized by poisedness\[[1](https://arxiv.org/html/2608.24963#bib.bib1)\]\. We use a concise version\.

###### Definition 3\(Λ\\Lambda\-poisedness \(informal\)\)\.

Let𝒫\\mathcal\{P\}be a polynomial model space and letY=\{y0,…,yp\}⊂ℬ⁡\(x,Δ\)Y=\\\{y\_\{0\},\\dots,y\_\{p\}\\\}\\subset\\mathcal\{B\}\\\!\\left\(x;\\Delta\\right\)\. The setYYis*Λ\\Lambda\-poised*inℬ⁡\(x,Δ\)\\mathcal\{B\}\\\!\\left\(x;\\Delta\\right\)if the interpolation system is well\-conditioned and the corresponding Lagrange polynomials are uniformly bounded byΛ\\Lambdaon the ball\[[1](https://arxiv.org/html/2608.24963#bib.bib1)\]\.

Poisedness supports fully linear/quadratic guarantees for polynomial models\. NNs lack an analogous linear\-system criterion, but geometry is still decisive: if samples lie near a low\-dimensional set \(e\.g\. mostly along a line\-search ray\), then the learner’s coverage in orthogonal directions is poor, often limiting its local generalization radius\.

Giovannelli et al\. explicitly discuss this issue for their surrogate\-enhanced FD\-BFGS step: line\-search points in Armijo backtracking are aligned, which is problematic for interpolation/regression models\. They therefore add only the first line\-search point to the dataset for interpolation/regression surrogates, while for NN surrogates they add all such points\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]\. This motivates explicit redundancy diagnostics\.

### 4\.3Exploration–Accuracy Trade\-offs

Trust regions frame the key trade\-off: smallΔ\\Deltaimproves local approximation; largeΔ\\Deltasupports exploration but may exceed the surrogate’s reliable range\.

To operationalize this trade\-off, we recommend coupling the trust\-region radius \(or stepsize scale in FD methods\) with surrogate diagnostics\. In Section[7](https://arxiv.org/html/2608.24963#S7)we define an effective generalization radius; in Section[9](https://arxiv.org/html/2608.24963#S9)we estimate it empirically through region\-size sweeps\.

##### Redundancy diagnostics, and what became of them\.

An earlier version of this work proposed two cheap geometric diagnostics for redundancy in the training data, an overlap ratio between consecutive training sets and a step\-spacing statistic along the backtracking ray, and built the framework of Section[8\.4](https://arxiv.org/html/2608.24963#S8.SS4)around them\. Section[9\.7](https://arxiv.org/html/2608.24963#S9.SS7)tests both against32023202outer iterations and finds that neither predicts whether a surrogate proposal will be accepted, as does a singular\-value diagnostic for directional coverage\. We therefore gate on measured model quality instead, through the statistics of \([29](https://arxiv.org/html/2608.24963#S8.E29)\)\. The definitions and the full refutation are given in the supplementary document, since a reader building such a method would otherwise reach for the same quantities; what the main text needs from this subsection is only that the redundancy question was asked, answered negatively, and replaced\.

## 5Offline and Online Learning Strategies

### 5\.1Offline Pre\-trained Models

Offline pre\-training is attractive when the optimization problems arise from a task family with shared structure \(e\.g\. repeated design optimization, parametric PDEs\)\. In such settings, one can learn transferable representations and reuse them across runs\. Neither primary source assumes such a distributional setting: Giovannelli et al\. train surrogates per test problem on locally sampled datasets\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\], and Taminiau et al\. train and update surrogates online within a single run\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\.

Therefore, within the scope of our synthesis, offline models should be treated as: \(i\) initializations \(warm starts\), or \(ii\) priors/regularizers that bias the online fit toward smoothness\.

A caveat is that offline generalization across objectives is fundamentally different from local interpolation on a fixed objective, and it should be tested explicitly \(Section[7\.3](https://arxiv.org/html/2608.24963#S7.SS3)\); a systematic cross\-task study is beyond the scope of this paper\.

### 5\.2Online Adaptive Training

Online learning updates the surrogate as new points are evaluated\.

##### Giovannelli et al\. \(approximation pipeline\)\.

For approximation studies, Giovannelli et al\. train NN surrogates by minimizing a mean\-squared error over a training datasetDDsampled in a ball and evaluate on a separate testing dataset sampled in the same ball\. They train five times per problem to compare activations and report that SiLU achieves the best performance across their benchmark suite\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]\. They use a feedforward network with 2 hidden layers and4​n4nneurons per layer, train it with Adam for 300 epochs, and use ReduceLROnPlateau \(factor 0\.8, patience 15\) based on the*testing*empirical risk\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]\. These details matter because the surrogate’s approximation quality is tied to training protocol as much as to architecture\.

##### Giovannelli et al\. \(online within FLE\-S\)\.

In their surrogate\-enhanced FD\-BFGS step \(FLE\-S\), they train an NN surrogate on the accumulated datasetDkD\_\{k\}once\|Dk\|\\left\\lvert D\_\{k\}\\right\\rvertexceeds a threshold, and they use a very small number of epochs to limit overhead: 5 epochs at the first surrogate call, then 1 epoch in subsequent calls\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]\. They choose learning rates by grid search and setζ\\zeta\(the threshold multiplier\) differently for polynomial/RBF vs NN surrogates\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]\. These choices embody a core online\-learning principle: amortize training, limit per\-iteration cost, and accept that the surrogate may be imperfect\.

##### Taminiau et al\. \(online, warm\-start, capped\)\.

Taminiau et al\. train a shallow NN surrogate with one hidden layer of width5​n5nand smooth activations, solve the training problem by L\-BFGS, and warm\-start subsequent trainings from previous parameters\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\. To limit cost and prevent ill\-conditioning, they cap dataset sizes:N≤10​\(n\+1\)N\\leq 10\(n\+1\)andM≤10M\\leq 10, removing the oldest points first if thresholds are exceeded\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\. These design choices are tightly coupled to their safeguarded\-assistance embedding\.

### 5\.3Hybrid Learning Strategies

A hybrid strategy can combine offline priors with online guarantees: \(i\) pretrain a feature extractor on a task family \(if available\), \(ii\) normalize and fine\-tune online in local regions, \(iii\) accept surrogate steps only when the true objective decreases sufficiently \(assistance with safeguards\)\. This approach is conceptually aligned with the framework of\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\.

## 6Model Architecture and Complexity

### 6\.1Neural Network Structures for Local Modeling

We consider feedforward NNs\. A one\-hidden\-layer \(shallow\) model is

mθ​\(x\)=W2​ϕ​\(W1​x\+b1\)\+b2,m\_\{\\theta\}\(x\)=W\_\{2\}\\,\\phi\(W\_\{1\}x\+b\_\{1\}\)\+b\_\{2\},\(14\)withW1∈ℝq×nW\_\{1\}\\in\\mathbb\{R\}^\{q\\times n\},b1∈ℝqb\_\{1\}\\in\\mathbb\{R\}^\{q\},W2∈ℝ1×qW\_\{2\}\\in\\mathbb\{R\}^\{1\\times q\},b2∈ℝb\_\{2\}\\in\\mathbb\{R\}, and elementwise activationϕ\\phi\. Taminiau et al\. useq=5​nq=5n\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\. Giovannelli et al\. use a 2\-hidden\-layer network with width4​n4nper layer in their approximation experiments\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]\.

### 6\.2Model Size and Approximation Capacity

Model size controls expressivity and training cost\. The shallow model \([14](https://arxiv.org/html/2608.24963#S6.E14)\) has

nw=q​n\+q\+q\+1=q⁡\(n\+2\)\+1n\_\{w\}=qn\+q\+q\+1=q\(n\+2\)\+1\(15\)parameters\. Withq=5​nq=5n,nw≈5​n2\+10​n\+1n\_\{w\}\\approx 5n^\{2\}\+10n\+1\. A two\-hidden\-layer network with widthsq1=q2=4​nq\_\{1\}=q\_\{2\}=4nhas

nw=\(q1​n\+q1\)\+\(q2​q1\+q2\)\+\(q2\+1\)≈20​n2\+12​n\+1\.n\_\{w\}=\(q\_\{1\}n\+q\_\{1\}\)\+\(q\_\{2\}q\_\{1\}\+q\_\{2\}\)\+\(q\_\{2\}\+1\)\\approx 20n^\{2\}\+12n\+1\.\(16\)Thus a deeper/wider architecture has substantially higher parameter count, potentially improving approximation at the cost of larger datasets and higher training time\.

##### Complexity relative to oracle cost\.

In DFO, the relevant question is not training cost in isolation but training cost*relative to the cost of a function evaluation*\. If evaluations are extremely expensive \(e\.g\. PDE solves\), expensive training may still be worthwhile\. In many numerical benchmarks, however, evaluation costs are moderate and training overhead can dominate\. Giovannelli et al\. emphasize that NN surrogates may be competitive in function\-evaluation counts but at high training costs\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]; conversely, Taminiau et al\. explicitly restrict model size and dataset size to keep training overhead bounded\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\. This motivates reporting both evaluation counts and wall\-clock time \(Section[9](https://arxiv.org/html/2608.24963#S9)\)\.

### 6\.3Regularization and Stability

Regularization stabilizes training in low\-data regimes:

λ​‖θ‖22\.\\lambda\\left\\lVert\\theta\\right\\rVert\_\{2\}^\{2\}\.\(17\)Giovannelli et al\. also normalize data via shifting/scaling intoℬ⁡\(0,1\)\\mathcal\{B\}\\\!\\left\(0;1\\right\), which improves numerical conditioning and can mitigate training pathologies \(e\.g\. vanishing gradients\)\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]\. Taminiau et al\. use weight decayλ=10−4\\lambda=10^\{\-4\}and employ Sobolev training, which introduces a gradient\-mismatch penalty and can be interpreted as a curvature \(diagonal Hessian\) regularizer under finite differences \(Section[7\.3](https://arxiv.org/html/2608.24963#S7.SS3)\)\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\.

## 7Generalization and Error Analysis

### 7\.1Local Generalization in Trust Regions

This section develops the second half of our central claim: surrogate trustworthiness is a*radius*question\. We formalize a radius\-aware notion of generalization suited to local modeling\.

###### Definition 4\(Local uniform error profiles\)\.

Givenx∈ℝnx\\in\\mathbb\{R\}^\{n\}and radiusΔ\>0\\Delta\>0, define

ef​\(x,Δ\):=supu∈ℬ⁡\(x,Δ\)\|f⁡\(u\)−m⁡\(u\)\|,eg​\(x,Δ\):=supu∈ℬ⁡\(x,Δ\)‖∇f​\(u\)−∇m​\(u\)‖2\.e\_\{f\}\(x,\\Delta\):=\\sup\_\{u\\in\\mathcal\{B\}\\\!\\left\(x;\\Delta\\right\)\}\\left\\lvert f\(u\)\-m\(u\)\\right\\rvert,\\qquad e\_\{g\}\(x,\\Delta\):=\\sup\_\{u\\in\\mathcal\{B\}\\\!\\left\(x;\\Delta\\right\)\}\\left\\lVert\\nabla f\(u\)\-\\nabla m\(u\)\\right\\rVert\_\{2\}\.\(18\)

###### Definition 5\(Effective local generalization radius\)\.

Fix a floorϵg,0≥0\\epsilon\_\{g,0\}\\geq 0and a rateκg\>0\\kappa\_\{g\}\>0, and sete¯g​\(δ\):=ϵg,0\+κg​δ\\bar\{e\}\_\{g\}\(\\delta\):=\\epsilon\_\{g,0\}\+\\kappa\_\{g\}\\delta, withe¯f\\bar\{e\}\_\{f\}defined analogously fromϵf,0\\epsilon\_\{f,0\}andκf​δ2\\kappa\_\{f\}\\delta^\{2\}\. The*effective generalization radius*aroundxxis

rgen​\(x\):=sup\{Δ\>0:ef​\(x,δ\)≤e¯f​\(δ\)​and​eg​\(x,δ\)≤e¯g​\(δ\)​for every​0<δ≤Δ\}\.r\_\{\\mathrm\{gen\}\}\(x\):=\\sup\\bigl\\\{\\Delta\>0:\\ e\_\{f\}\(x,\\delta\)\\leq\\bar\{e\}\_\{f\}\(\\delta\)\\ \\text\{and\}\\ e\_\{g\}\(x,\\delta\)\\leq\\bar\{e\}\_\{g\}\(\\delta\)\\ \\text\{ for every \}0<\\delta\\leq\\Delta\\bigr\\\}\.\(19\)

The floor is not cosmetic\. Withϵg,0=0\\epsilon\_\{g,0\}=0the toleranceκg​δ\\kappa\_\{g\}\\deltavanishes asδ→0\\delta\\to 0whileeg​\(x,⋅\)e\_\{g\}\(x,\\cdot\)is nondecreasing and, for a trained network, bounded below by its own training error, so any surrogate with nonzero error atxxfails the test at every small enoughδ\\deltaand hasrgen​\(x\)=0r\_\{\\mathrm\{gen\}\}\(x\)=0exactly\. That is not an edge case here: Section[9\.3](https://arxiv.org/html/2608.24963#S9.SS3)measures the floor and finds it does not vanish\. Against the classical fully linear standard, which isϵg,0=0\\epsilon\_\{g,0\}=0, the radius of every network we fit is therefore zero, and the content of that section is how largeϵg,0\\epsilon\_\{g,0\}must be before the radius becomes positive and useful\. Quantifying over allδ≤Δ\\delta\\leq\\Deltais what makes the qualifying set an interval: both sides increase withδ\\delta, so testing atΔ\\Deltaalone would admit a surrogate that passes there and fails on a smaller ball inside\.

Definingrgenr\_\{\\mathrm\{gen\}\}makes explicit what is often implicit in surrogate\-based optimization: a surrogate may be “good enough” only within a limited region\. Giovannelli et al\. report that reducing ball radii can improve approximation accuracy\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\], consistent with finitergenr\_\{\\mathrm\{gen\}\}that depends on curvature and sample coverage\.

### 7\.2Error Bounds for Neural Approximation

The two propositions below connect training accuracy and sample coverage to uniform errors\. Neither uses anything specific to neural networks: the first is Lipschitz continuity and a fill distance combined by the triangle inequality, the second the standard inexact\-gradient descent condition\. We offer them as an elementary analytical interpretation of what the experiments measure, not as a theory of learned models, and we do not verify that a trained network satisfies their hypotheses\.

###### Assumption 2\(Smoothness\)\.

ffis continuously differentiable onℬ⁡\(x,Δ\)\\mathcal\{B\}\\\!\\left\(x;\\Delta\\right\)and∇f\\nabla fisLfL\_\{f\}\-Lipschitz onℬ⁡\(x,Δ\)\\mathcal\{B\}\\\!\\left\(x;\\Delta\\right\)\. The surrogatemmis continuously differentiable and∇m\\nabla misLmL\_\{m\}\-Lipschitz onℬ⁡\(x,Δ\)\\mathcal\{B\}\\\!\\left\(x;\\Delta\\right\)\.

###### Assumption 3\(Sample coverage via fill distance\)\.

A sample setY⊂ℬ⁡\(x,Δ\)Y\\subset\\mathcal\{B\}\\\!\\left\(x;\\Delta\\right\)satisfies

h⁡\(Y,ℬ⁡\(x,Δ\)\):=supu∈ℬ⁡\(x,Δ\)miny∈Y⁡‖u−y‖≤ε\.h\(Y;\\mathcal\{B\}\\\!\\left\(x;\\Delta\\right\)\):=\\sup\_\{u\\in\\mathcal\{B\}\\\!\\left\(x;\\Delta\\right\)\}\\min\_\{y\\in Y\}\\left\\lVert u\-y\\right\\rVert\\leq\\varepsilon\.\(20\)

###### Assumption 4\(Pointwise training accuracy\)\.

Training yields pointwise bounds onYY:

maxy∈Y⁡\|f⁡\(y\)−m⁡\(y\)\|≤ϵf,maxy∈Y⁡‖∇f​\(y\)−∇m​\(y\)‖≤ϵg\.\\max\_\{y\\in Y\}\\left\\lvert f\(y\)\-m\(y\)\\right\\rvert\\leq\\epsilon\_\{f\},\\qquad\\max\_\{y\\in Y\}\\left\\lVert\\nabla f\(y\)\-\\nabla m\(y\)\\right\\rVert\\leq\\epsilon\_\{g\}\.\(21\)

###### Proposition 1\(Uniform gradient error bound\)\.

Under Assumptions[2](https://arxiv.org/html/2608.24963#Thmassumption2)–[4](https://arxiv.org/html/2608.24963#Thmassumption4),

eg​\(x,Δ\)≤ϵg\+\(Lf\+Lm\)​ε\.e\_\{g\}\(x,\\Delta\)\\leq\\epsilon\_\{g\}\+\(L\_\{f\}\+L\_\{m\}\)\\varepsilon\.\(22\)

###### Proof\.

Fixu∈ℬ⁡\(x,Δ\)u\\in\\mathcal\{B\}\\\!\\left\(x;\\Delta\\right\)and choosey∈Yy\\in Ywith‖u−y‖≤ε\\left\\lVert u\-y\\right\\rVert\\leq\\varepsilon\. Then

‖∇f​\(u\)−∇m​\(u\)‖≤‖∇f​\(u\)−∇f​\(y\)‖\+‖∇f​\(y\)−∇m​\(y\)‖\+‖∇m​\(y\)−∇m​\(u\)‖\.\\left\\lVert\\nabla f\(u\)\-\\nabla m\(u\)\\right\\rVert\\leq\\left\\lVert\\nabla f\(u\)\-\\nabla f\(y\)\\right\\rVert\+\\left\\lVert\\nabla f\(y\)\-\\nabla m\(y\)\\right\\rVert\+\\left\\lVert\\nabla m\(y\)\-\\nabla m\(u\)\\right\\rVert\.Apply Lipschitz bounds and Assumption[4](https://arxiv.org/html/2608.24963#Thmassumption4)\. ∎

###### Proposition 2\(Descent direction under inexact gradients\)\.

Letgm​\(x\)=∇m​\(x\)g\_\{m\}\(x\)=\\nabla m\(x\)\. If‖gm​\(x\)−∇f​\(x\)‖≤ϵg\\left\\lVert g\_\{m\}\(x\)\-\\nabla f\(x\)\\right\\rVert\\leq\\epsilon\_\{g\}, then

⟨∇f​\(x\),−gm​\(x\)⟩≤−‖∇f​\(x\)‖2​\(‖∇f​\(x\)‖2−ϵg\)\.\\left\\langle\\nabla f\(x\),\-g\_\{m\}\(x\)\\right\\rangle\\leq\-\\left\\lVert\\nabla f\(x\)\\right\\rVert\_\{2\}\\big\(\\left\\lVert\\nabla f\(x\)\\right\\rVert\_\{2\}\-\\epsilon\_\{g\}\\big\)\.\(23\)In particular, if‖∇f​\(x\)‖2\>ϵg\\left\\lVert\\nabla f\(x\)\\right\\rVert\_\{2\}\>\\epsilon\_\{g\}, then−gm​\(x\)\-g\_\{m\}\(x\)is a strict descent direction forffatxx\.

###### Proof\.

Use⟨∇f,−gm⟩=−‖∇f‖2−⟨∇f,gm−∇f⟩\\left\\langle\\nabla f,\-g\_\{m\}\\right\\rangle=\-\\left\\lVert\\nabla f\\right\\rVert^\{2\}\-\\left\\langle\\nabla f,g\_\{m\}\-\\nabla f\\right\\rangleand Cauchy–Schwarz\. ∎

Theϵg\\epsilon\_\{g\}of Proposition[2](https://arxiv.org/html/2608.24963#Thmproposition2)is a bound on the gradient error at the single pointxx, whereas theϵg\\epsilon\_\{g\}of Assumption[4](https://arxiv.org/html/2608.24963#Thmassumption4)bounds it on the training set\. Any valid upper bound may be substituted, and the one we can measure is the uniform erroreg​\(x,Δ\)e\_\{g\}\(x,\\Delta\)of Definition[4](https://arxiv.org/html/2608.24963#Thmdefinition4), which dominates the pointwise error on the ball; the substitution is therefore conservative, and it is how Section[9\.3](https://arxiv.org/html/2608.24963#S9.SS3)reads the proposition against data\.

### 7\.3Robustness under Noise

Taminiau et al\. provide a key structural insight: Sobolev learning with finite\-difference gradients induces a curvature penalty\.

###### Proposition 3\(Sobolev learning as curvature regularization\)\.

\(Adapted from Proposition 3\.1 in\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\.\) Letmθ:ℝn→ℝm\_\{\\theta\}:\\mathbb\{R\}^\{n\}\\to\\mathbb\{R\}be twice continuously differentiable\. Fixz∈ℝnz\\in\\mathbb\{R\}^\{n\}andh\>0h\>0such that

mθ\(z\)=f\(z\),mθ\(z\+heℓ\)=f\(z\+heℓ\),ℓ=1,…,n\.m\_\{\\theta\}\(z\)=f\(z\),\\qquad m\_\{\\theta\}\(z\+he\_\{\\ell\}\)=f\(z\+he\_\{\\ell\}\),\\ \\ell=1,\\dots,n\.\(24\)Definegh​\(z\)g\_\{h\}\(z\)by forward differences \([4](https://arxiv.org/html/2608.24963#S2.E4)\)\. Then there existζ1,…,ζn∈\[0,1\]\\zeta\_\{1\},\\dots,\\zeta\_\{n\}\\in\[0,1\]such that

‖gh​\(z\)−∇mθ​\(z\)‖22=h24​∑ℓ=1n\(\[∇2mθ​\(z\+h​ζℓ​eℓ\)\]ℓ​ℓ\)2\.\\left\\lVert g\_\{h\}\(z\)\-\\nabla m\_\{\\theta\}\(z\)\\right\\rVert\_\{2\}^\{2\}=\\frac\{h^\{2\}\}\{4\}\\sum\_\{\\ell=1\}^\{n\}\\left\(\[\\nabla^\{2\}m\_\{\\theta\}\(z\+h\\zeta\_\{\\ell\}e\_\{\\ell\}\)\]\_\{\\ell\\ell\}\\right\)^\{2\}\.\(25\)

Proposition[3](https://arxiv.org/html/2608.24963#Thmproposition3)suggests why the gradient\-matching term helps: it penalises diagonal curvature along coordinate directions when the model interpolates the finite\-difference stencil\. The hypothesis is worth naming, because our training does not satisfy it\. Equation \([27](https://arxiv.org/html/2608.24963#S8.E27)\) is a regression with weight decay, not an interpolation, so \([24](https://arxiv.org/html/2608.24963#S7.E24)\) holds only approximately and the identity becomes an approximate one whose error we do not control\. The proposition is therefore an idealised mechanism consistent with the acceptance\-rate measurement of Section[9\.4](https://arxiv.org/html/2608.24963#S9.SS4), not a causal account of it\.

##### On noise in gradient targets\.

Gradient targetsg⁡\(z\)g\(z\)are computed from noisy function evaluations and inherit noise that is amplified byh−1h^\{\-1\}\. Weight decay, dataset capping, and warm\-start training \(all used in\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\) mitigate overfitting to noisy gradients, but the trade\-offs are problem dependent\. This motivates ablations that vary \(i\) the number of gradient pointsMM, \(ii\) the finite\-difference stepsize policy, and \(iii\) whether the value term is included in the loss\.

##### Cross\-task generalization\.

In some applications, one cares about generalization across objective functions \(tasks\) under similar sampling geometry\. This notion is distinct from local generalization on a single objective\. We recommend testing cross\-task generalization empirically by sampling multiple functions \(e\.g\. from CUTEst categories\) using the same sampling design and comparing achievableef,ege\_\{f\},e\_\{g\}profiles; such a study is left as future work\.

## 8Algorithm Design

### 8\.1Learning\-Augmented Trust\-Region Framework

This section develops the first half of our central claim, the role of the embedding\. We adopt a role\-based taxonomy that reconciles the two primary sources\. Let a base DFO method produce iterates and data; a surrogate model is trained on that data; the surrogate is used either to replace a core ingredient or to assist with safeguarded proposals\.

### 8\.2Neural Model Construction and Update

We recount the training objectives used in the primary sources\.

##### Value\-only training \(Giovannelli et al\.\)\.

Given a datasetD=\{\(xi,f⁡\(xi\)\)\}i=0ND=\\\{\(x\_\{i\},f\(x\_\{i\}\)\)\\\}\_\{i=0\}^\{N\}, Giovannelli et al\. trainfNN​\(⋅,w\)f\_\{\\mathrm\{NN\}\}\(\\cdot;w\)by minimizing the empirical risk

minw∈ℝnw⁡L⁡\(w,D\):=∑i=0N\(f⁡\(xi\)−fNN​\(xi,w\)\)2,\\min\_\{w\\in\\mathbb\{R\}^\{n\_\{w\}\}\}L\(w;D\):=\\sum\_\{i=0\}^\{N\}\\big\(f\(x\_\{i\}\)\-f\_\{\\mathrm\{NN\}\}\(x\_\{i\};w\)\\big\)^\{2\},\(26\)and assess generalization via a testing dataset sampled from the same ball\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]\. They normalize data by shifting/scaling:

xi↦xi−x0Δ,Δ=max1≤i≤N⁡‖xi−x0‖,x\_\{i\}\\mapsto\\frac\{x\_\{i\}\-x\_\{0\}\}\{\\Delta\},\\quad\\Delta=\\max\_\{1\\leq i\\leq N\}\\left\\lVert x\_\{i\}\-x\_\{0\}\\right\\rVert,so the transformed dataset lies inℬ⁡\(0,1\)\\mathcal\{B\}\\\!\\left\(0;1\\right\)\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]\.

##### Sobolev loss \(Taminiau et al\.\)\.

Taminiau et al\. defineF=\{\(yi,f⁡\(yi\)\)\}i=1NF=\\\{\(y\_\{i\},f\(y\_\{i\}\)\)\\\}\_\{i=1\}^\{N\}andG=\{\(zj,g⁡\(zj\)\)\}j=1MG=\\\{\(z\_\{j\},g\(z\_\{j\}\)\)\\\}\_\{j=1\}^\{M\}, withg⁡\(zj\)≈∇f​\(zj\)g\(z\_\{j\}\)\\approx\\nabla f\(z\_\{j\}\)from finite differences, and train a differentiable surrogatemθm\_\{\\theta\}by

minθ⁡LF,G​\(θ\):=1N​∑i=1N\(mθ​\(yi\)−f⁡\(yi\)\)2\+1M​∑j=1M‖∇mθ​\(zj\)−g⁡\(zj\)‖22\+λ​‖θ‖22\.\\min\_\{\\theta\}\\ L\_\{F,G\}\(\\theta\):=\\frac\{1\}\{N\}\\sum\_\{i=1\}^\{N\}\\big\(m\_\{\\theta\}\(y\_\{i\}\)\-f\(y\_\{i\}\)\\big\)^\{2\}\+\\frac\{1\}\{M\}\\sum\_\{j=1\}^\{M\}\\left\\lVert\\nabla m\_\{\\theta\}\(z\_\{j\}\)\-g\(z\_\{j\}\)\\right\\rVert\_\{2\}^\{2\}\+\\lambda\\left\\lVert\\theta\\right\\rVert\_\{2\}^\{2\}\.\(27\)They setλ=10−4\\lambda=10^\{\-4\}and solve \([27](https://arxiv.org/html/2608.24963#S8.E27)\) with low\-memory BFGS \(L\-BFGS\), stopping at iterationKKif

‖∇LF,G​\(θK\)‖2≤10−6​max⁡\{1,‖∇LF,G​\(θ0\)‖2\},\\left\\lVert\\nabla L\_\{F,G\}\(\\theta\_\{K\}\)\\right\\rVert\_\{2\}\\leq 10^\{\-6\}\\max\\\{1,\\left\\lVert\\nabla L\_\{F,G\}\(\\theta\_\{0\}\)\\right\\rVert\_\{2\}\\\},\(28\)or if a budget of 1000 iterations is reached\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\. They initialize the first surrogate with He initialization \(SoftPlus/SiLU\) or Glorot initialization \(sigmoid\) and warm\-start subsequent surrogates from the last trained model\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\. They cap dataset sizesN≤10​\(n\+1\)N\\leq 10\(n\+1\)andM≤10M\\leq 10, removing the oldest points first\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\.

##### Value\-removal variants and identifiability\.

If the value term is removed and one trains only on gradients, the learned model is identifiable only up to an additive constant\. This can be acceptable if only directions are used; however, any acceptance rule based on predicted decreases requires fixing the offset \(e\.g\. via anchoringmθ​\(xref\)=f⁡\(xref\)m\_\{\\theta\}\(x\_\{\\mathrm\{ref\}\}\)=f\(x\_\{\\mathrm\{ref\}\}\)\)\. This is relevant in workflows where one emphasizes gradient accuracy over value accuracy \(e\.g\. some shape\-optimization pipelines\)\.

### 8\.3Step Computation and Acceptance Criteria

To keep the main text light, we present pseudocode only for the safeguarded assistance procedure \(Algorithm[1](https://arxiv.org/html/2608.24963#alg1)\), which is the embedding at the heart of our central claim\. The base method, due to\[[11](https://arxiv.org/html/2608.24963#bib.bib11)\]and simplified by Taminiau et al\., its surrogate\-augmented driver, and the FLE\-S replacement step of Giovannelli et al\. are summarized in words below and reproduced in full in the supplementary documentsupplement\-pseudocode\.pdf, which accompanies thecode/directory\.

#### Assistance: safeguarded surrogate procedure

Algorithm 1Safeguarded surrogate steps \(assist\), adapted from Algorithm 1 in\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]1:

v∈ℝnv\\in\\mathbb\{R\}^\{n\}; oracle access to

ff; datasets

F=\{\(yi,f⁡\(yi\)\)\}i=1NF=\\\{\(y\_\{i\},f\(y\_\{i\}\)\)\\\}\_\{i=1\}^\{N\}and

G=\{\(zj,g⁡\(zj\)\)\}j=1MG=\\\{\(z\_\{j\},g\(z\_\{j\}\)\)\\\}\_\{j=1\}^\{M\}with

\(v,f⁡\(v\)\)∈F\(v,f\(v\)\)\\in F; constants

σ,ρ,λ,γ,ϵ\>0\\sigma,\\rho,\\lambda,\\gamma,\\epsilon\>0
2:

v\+v^\{\+\}\(best point found\),

t\+t^\{\+\}\(number of successful surrogate steps\), and

F\+F^\{\+\}\(newly evaluated points\)

3:Train a

C1C^\{1\}surrogate

mθm\_\{\\theta\}by approximately solving \([27](https://arxiv.org/html/2608.24963#S8.E27)\)

4:

v0←vv\_\{0\}\\leftarrow v,

L0←σL\_\{0\}\\leftarrow\\sigma,

F\+←∅F^\{\+\}\\leftarrow\\emptyset,

t←0t\\leftarrow 0
5:whiletruedo

6:Find smallest integer

ℓt≥0\\ell\_\{t\}\\geq 0such that with

v^t=vt−12ℓt​Lt∇mθ\(vt\),\\hat\{v\}\_\{t\}=v\_\{t\}\-\\frac\{1\}\{2^\{\\ell\_\{t\}\}L\_\{t\}\}\\,\\nabla m\_\{\\theta\}\(v\_\{t\}\),the surrogate decrease condition holds:

mθ​\(vt\)−mθ​\(v^t\)≥ρ2ℓt​Lt​‖∇mθ​\(vt\)‖22m\_\{\\theta\}\(v\_\{t\}\)\-m\_\{\\theta\}\(\\hat\{v\}\_\{t\}\)\\ \\geq\\ \\frac\{\\rho\}\{2^\{\\ell\_\{t\}\}L\_\{t\}\}\\,\\left\\lVert\\nabla m\_\{\\theta\}\(v\_\{t\}\)\\right\\rVert\_\{2\}^\{2\}
7:Evaluate

f⁡\(v^t\)f\(\\hat\{v\}\_\{t\}\)and set

F\+←F\+∪\{\(v^t,f⁡\(v^t\)\)\}F^\{\+\}\\leftarrow F^\{\+\}\\cup\\\{\(\\hat\{v\}\_\{t\},f\(\\hat\{v\}\_\{t\}\)\)\\\}
8:if

f⁡\(vt\)−f⁡\(v^t\)≥1γ​σ​ϵ2f\(v\_\{t\}\)\-f\(\\hat\{v\}\_\{t\}\)\\ \\geq\\ \\frac\{1\}\{\\gamma\\sigma\}\\epsilon^\{2\}then

9:Accept:

vt\+1←v^tv\_\{t\+1\}\\leftarrow\\hat\{v\}\_\{t\},

Lt\+1←2ℓt−1​LtL\_\{t\+1\}\\leftarrow 2^\{\\ell\_\{t\}\-1\}L\_\{t\},

t←t\+1t\\leftarrow t\+1
10:else

11:

t\+←tt^\{\+\}\\leftarrow t,

v\+←vtv^\{\+\}\\leftarrow v\_\{t\}andbreak

12:endif

13:endwhile

14:return

\(v\+,t\+,F\+\)\(v^\{\+\},t^\{\+\},F^\{\+\}\)

#### Taminiau et al\.: base method and surrogate\-augmented method

Taminiau et al\. integrate the surrogate procedure into a finite\-difference Armijo\-type base method \(a derivative\-free gradient method\); the pseudocode is reproduced from Algorithms 2 and 3 of\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\], with minimal notational changes, in the supplementary documentsupplement\-pseudocode\.pdf\. In brief, the base method computes a forward\-difference gradient whose stepsize is tied to the current scaleσk\\sigma\_\{k\}, accepts trial points by an Armijo\-type test, and doubles the scale upon rejection; the augmented variant calls the surrogate procedure \(Algorithm[1](https://arxiv.org/html/2608.24963#alg1)\) after every accepted base step and continues from the best point found\.

As reproduced there, these algorithms carry no termination test, and this is inherited from the source rather than an omission on our part: the analysis of\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]covers only the iterations before the hitting timeT⁡\(ϵ\)=inf\{k:‖∇f​\(xk\)‖2≤ϵ\}T\(\\epsilon\)=\\inf\\\{k:\\left\\lVert\\nabla f\(x\_\{k\}\)\\right\\rVert\_\{2\}\\leq\\epsilon\\\}, and Lemma 2\.2 there establishes that iterationkkis well defined precisely when‖∇f​\(xk\)‖2\>ϵ\\left\\lVert\\nabla f\(x\_\{k\}\)\\right\\rVert\_\{2\}\>\\epsilon\. Near a stationary point the inner loop can fail both tests indefinitely, shrinkinghih\_\{i\}without ever returning\. The pseudocode should therefore be read as describing the iterationsk<T⁡\(ϵ\)k<T\(\\epsilon\)only; any executable version needs an external stopping rule\. Our implementation \(Section[9](https://arxiv.org/html/2608.24963#S9)\) caps the inner index ati≤40i\\leq 40, treats exhaustion of that cap as an approximate\-stationarity or stall exit, and additionally bounds every run by the evaluation budget\.

#### Replacement: surrogate gradient within FD\-BFGS in FLE\-S

Giovannelli et al\. propose FLE\-S by modifying the FD\-BFGS \(Full\-Eval\) step: once the datasetDkD\_\{k\}exceeds a threshold, they build a surrogate and set the gradient approximation to the surrogate gradient\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]\. They also add a random point inℬ⁡\(xk,0\.1\)\\mathcal\{B\}\\\!\\left\(x\_\{k\};0\.1\\right\)before training/building the surrogate, and they updateDkD\_\{k\}differently depending on surrogate type \(add\-first for poly/RBF; add\-all for NN\)\. The corresponding pseudocode is reproduced in the supplementary documentsupplement\-pseudocode\.pdf\.

### 8\.4From Diagnostics to Design: the ARAS Framework

The analysis so far suggests a design recipe\. Monitor how far the current surrogate can be trusted, and let that decide how much of the step it is allowed to influence, rather than fixing the surrogate’s role in advance\. We call this the*Adaptive Radius\-Aware Surrogate*framework \(ARAS\)\.

The open question in such a design is which quantity to gate on\. Our first version used the geometric diagnostics of Section[4\.3](https://arxiv.org/html/2608.24963#S4.SS3), the overlap ratiorkr\_\{k\}and the step\-spacing statisticNspN\_\{\\mathrm\{sp\}\}\. Section[9\.7](https://arxiv.org/html/2608.24963#S9.SS7)tests those against data and finds their effect on acceptance an order of magnitude too small to gate on, so we do not gate on them\. The same section tests a second family of candidates, computed from quantities the method already holds, and identifies two that do predict whether the next surrogate proposal will be accepted: the ratio‖∇mk​\(xk\)‖/‖gkFD‖\\\|\\nabla m\_\{k\}\(x\_\{k\}\)\\\|/\\\|g\_\{k\}^\{\\mathrm\{FD\}\}\\\|between the surrogate and finite\-difference gradient norms, and the training residual of the surrogate on its own data\. Sorted into quartiles by the first of these, the empirical acceptance rate over the31963196outer iterations at which both statistics are defined falls monotonically through0\.4840\.484,0\.3430\.343,0\.1060\.106and0\.0450\.045, an elevenfold drop from the lowest quartile ofκk\\kappa\_\{k\}to the highest\. The framework therefore gates on

κk=‖∇mk​\(xk\)‖‖gkFD‖,ϱk=\(1\|Fk\|​∑\(y,f\)∈Fk\(mk​\(y\)−f\)2\)1/2,\\kappa\_\{k\}\\;=\\;\\frac\{\\\|\\nabla m\_\{k\}\(x\_\{k\}\)\\\|\}\{\\\|g\_\{k\}^\{\\mathrm\{FD\}\}\\\|\},\\qquad\\varrho\_\{k\}\\;=\\;\\Bigl\(\\tfrac\{1\}\{\|F\_\{k\}\|\}\\textstyle\\sum\_\{\(y,f\)\\in F\_\{k\}\}\(m\_\{k\}\(y\)\-f\)^\{2\}\\Bigr\)^\{1/2\},\(29\)a scale disagreement and a fit residual, both cheap and both validated\.

The blended direction admits a guarantee, and the quantity it needs is not the surrogate’s accuracy – which Section[9\.3](https://arxiv.org/html/2608.24963#S9.SS3)shows is not boundable below any useful floor – but its scale relative to the finite\-difference gradient, which is exactly theκk\\kappa\_\{k\}the gate already computes at no oracle cost\.

Bounding the blending weight by

0≤αk≤1−θ1\+κk,θ∈\(0,1\),0\\ \\leq\\ \\alpha\_\{k\}\\ \\leq\\ \\frac\{1\-\\theta\}\{1\+\\kappa\_\{k\}\},\\qquad\\theta\\in\(0,1\),\(30\)gives three things, proved in the supplementary document\. The blended direction satisfies⟨gk,dk⟩≤−θ​‖gk‖2\\left\\langle g\_\{k\},d\_\{k\}\\right\\rangle\\leq\-\\theta\\left\\lVert g\_\{k\}\\right\\rVert^\{2\}and‖dk‖≤2​‖gk‖\\left\\lVert d\_\{k\}\\right\\rVert\\leq 2\\left\\lVert g\_\{k\}\\right\\rVert\(Proposition\)\. Backtracking on the true objective terminates in⌈log2⁡\(4​L​β¯/\(\(1−c\)​θ\)\)⌉\+\\lceil\\log\_\{2\}\(4L\\bar\{\\beta\}/\(\(1\-c\)\\theta\)\)\\rceil\_\{\+\}halvings withβk≥min⁡\{β¯,\(1−c\)​θ/\(8​L\)\}\\beta\_\{k\}\\geq\\min\\\{\\bar\{\\beta\},\(1\-c\)\\theta/\(8L\)\\\}, provided the working scale satisfies

2ik​σk≥L\(1−c\)​θ,c∈\(0,1\),2^\{i\_\{k\}\}\\sigma\_\{k\}\\ \\geq\\ \\frac\{L\}\{\(1\-c\)\\theta\},\\qquad c\\in\(0,1\),\(31\)which is equivalent to the finite\-difference error obeyingek≤\(1−c\)​θ4​‖gk‖e\_\{k\}\\leq\\frac\{\(1\-c\)\\theta\}\{4\}\\left\\lVert g\_\{k\}\\right\\rVert\(Lemma\)\. And ARAS then requires at most𝒪⁡\(n​ϵ−2\)\\mathcal\{O\}\(n\\epsilon^\{\-2\}\)evaluations to reach‖∇f​\(xk\)‖2≤ϵ\\left\\lVert\\nabla f\(x\_\{k\}\)\\right\\rVert\_\{2\}\\leq\\epsilon, the base method’s order, with a constant no worse: atαk≡0\\alpha\_\{k\}\\equiv 0the iteration is the base method itself \(Theorem\)\.

Three points\. Condition \([30](https://arxiv.org/html/2608.24963#S8.E30)\) is computable:κk\\kappa\_\{k\}costs nothing beyond the surrogate gradient already formed\. It is not vacuous: over the32023202logged iterations the cap has median0\.2470\.247atθ=1/4\\theta=1/4and permitsαk≥0\.1\\alpha\_\{k\}\\geq 0\.1at63%63\\%of them\. And \([31](https://arxiv.org/html/2608.24963#S8.E31)\) is about the base method rather than the network: withhi=2​ϵ/\(5​n​2i​σk\)h\_\{i\}=2\\epsilon/\(5\\sqrt\{n\}\\,2^\{i\}\\sigma\_\{k\}\)the forward\-difference error obeysek≤L​ϵ/\(5⋅2ik​σk\)e\_\{k\}\\leq L\\epsilon/\(5\\cdot 2^\{i\_\{k\}\}\\sigma\_\{k\}\)while the inner loop exits only with‖gk‖≥4​ϵ/5\\left\\lVert g\_\{k\}\\right\\rVert\\geq 4\\epsilon/5, soek/‖gk‖≤L/\(4⋅2ik​σk\)e\_\{k\}/\\left\\lVert g\_\{k\}\\right\\rVert\\leq L/\(4\\cdot 2^\{i\_\{k\}\}\\sigma\_\{k\}\); the condition asks the working scale to have grown past a fixed multiple ofLL, which is what theσk\\sigma\_\{k\}update does\. What the theorem does not give is an*improvement*in rate\.

The gate is calibrated rather than guessed, and the answer is not the expected one: on a cheap deterministic oracle the calibrated gate never fires\. A rejected proposal costs one evaluation, an accepted one saves then\+1n\+1of a finite\-difference gradient, so gating pays only where acceptance falls below1/\(n\+1\)1/\(n\+1\), and on this benchmark it does not \(Section[9\.7](https://arxiv.org/html/2608.24963#S9.SS7)\)\.

##### Running the blended step\.

A theorem about an algorithm nobody has run is worth little, so we implemented the framework \(pseudocode in the supplementary document\) withθ=1/4\\theta=1/4,c=1/2c=1/2andαmax=1/2\\alpha\_\{\\max\}=1/2, taking the blended step at every outer iteration and following it with the same safeguarded loop as elsewhere\. Over the117117instances it solves7474atτ=10−5\\tau=10^\{\-5\}against9999for plain safeguarded assistance and6868for the base method, with464464of13391339blended steps accepted at a meanαk\\alpha\_\{k\}of0\.1990\.199\. The blended step is convergent, as the Theorem says, and it is not worth its price: each attempt spends backtracking evaluations on the true objective, and the two thirds that fail buy nothing the surrogate loop would not have found more cheaply\.

The gate calibration says the same from the other direction, and the two together settle what ARAS should be: the instance surviving both tests hasαk≡0\\alpha\_\{k\}\\equiv 0and is exactly the safeguarded assistance benchmarked throughout\. We keep the theorem because it certifies that opening the gate cannot destroy convergence, and report the experiment because it says when to leave it shut\.

### 8\.5Convergence Analysis

##### Assistance paradigm: finite termination and surrogate gain\.

Taminiau et al\. prove that the surrogate procedure terminates in finite time under a lower\-bounded objective assumption, and they derive evaluation\-complexity improvements quantified by a surrogate\-gain factor\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\. We restate the essential components\.

###### Assumption 5\(Lower bounded objective\)\.

ffis bounded below byflow∈ℝf\_\{\\mathrm\{low\}\}\\in\\mathbb\{R\}\.

Algorithm[1](https://arxiv.org/html/2608.24963#alg1)then terminates after finitely many successful surrogate steps: each accepted step decreasesffby at leastϵ2/\(γ​σ\)\\epsilon^\{2\}/\(\\gamma\\sigma\), so infinitely many acceptances would driveffbelowflowf\_\{\\mathrm\{low\}\}\(Lemma 2\.1 of\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\)\.

###### Theorem 1\(Surrogate gain factor; Theorem 2\.4 of\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\)\.

Letffbe bounded below \(Assumption[5](https://arxiv.org/html/2608.24963#Thmassumption5)\) withLL\-Lipschitz gradient, let

T⁡\(ϵ\):=inf\{k∈ℕ:‖∇f​\(xk\)‖2≤ϵ\},T\(\\epsilon\):=\\inf\\\{k\\in\\mathbb\{N\}:\\ \\left\\lVert\\nabla f\(x\_\{k\}\)\\right\\rVert\_\{2\}\\leq\\epsilon\\\},and lettkt\_\{k\}be the number of successful surrogate steps taken at outer iterationkk\. Define the*average*number of successful surrogate steps per outer iteration over the firstTTiterations,

S⁡\(T\):=1T​∑k=0T−1tk,S\(T\):=\\frac\{1\}\{T\}\\sum\_\{k=0\}^\{T\-1\}t\_\{k\},\(32\)and the surrogate gain

η⁡\(S\)=1\+S2​\(n\+1\)1\+S\.\\eta\(S\)=\\frac\{1\+\\frac\{S\}\{2\(n\+1\)\}\}\{1\+S\}\.\(33\)Then the numberF​E​\(ϵ\)FE\(\\epsilon\)of evaluations offfrequired by the surrogate\-augmented method satisfies

F​E​\(ϵ\)≤4​η​\(S⁡\(T⁡\(ϵ\)\)\)​\(n\+1\)​Cmax​\(f⁡\(x0\)−flow\)​ϵ−2\+log2⁡\(σmaxσ0\)​\(n\+1\)\+T⁡\(ϵ\),FE\(\\epsilon\)\\ \\leq\\ 4\\,\\eta\\\!\\left\(S\(T\(\\epsilon\)\)\\right\)\(n\+1\)\\,C\_\{\\max\}\\,\(f\(x\_\{0\}\)\-f\_\{\\mathrm\{low\}\}\)\\,\\epsilon^\{\-2\}\\;\+\\;\\log\_\{2\}\\\!\\left\(\\tfrac\{\\sigma\_\{\\max\}\}\{\\sigma\_\{0\}\}\\right\)\(n\+1\)\\;\+\\;T\(\\epsilon\),\(34\)withσmax\\sigma\_\{\\max\}andCmaxC\_\{\\max\}the problem constants defined in\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\.

The statement should be read with three qualifications\. First,S⁡\(T\)S\(T\)in \([32](https://arxiv.org/html/2608.24963#S8.E32)\) is an*average per outer iteration*, not a total count:η\\etais driven down only if the surrogate succeeds repeatedly at a typical iteration, andS⁡\(T\)≥nS\(T\)\\geq nis what yields the favourable regimeη≤3/\(2​\(n\+1\)\)\\eta\\leq 3/\(2\(n\+1\)\)\. Second, the gain multiplies only the leadingϵ−2\\epsilon^\{\-2\}term of \([34](https://arxiv.org/html/2608.24963#S8.E34)\); the remaining two terms are unaffected, so the improvement is in the constant of the dominant term rather than in the𝒪⁡\(n​ϵ−2\)\\mathcal\{O\}\(n\\epsilon^\{\-2\}\)rate itself\. Sinceη⁡\(0\)=1\\eta\(0\)=1, the bound never degrades relative to the base method\. Third, and this limits how the result may be used,S⁡\(T⁡\(ϵ\)\)S\(T\(\\epsilon\)\)is measured on the very run the bound describes\. The statement is therefore a posteriori: before running, the only value ofη\\etaone is entitled to assume isη⁡\(0\)=1\\eta\(0\)=1, at which the bound is the base method’s own\. It certifies that observed surrogate success translates into a smaller constant, not that success will occur\.

##### Replacement paradigm: why guarantees are harder\.

In replacement methods, surrogate gradient errors propagate through quasi\-Newton updates and line\-search decisions: an error ingkg\_\{k\}contaminates both the directionpk=−Hk​gkp\_\{k\}=\-H\_\{k\}g\_\{k\}and the curvature pair\(sk,yk\)\(s\_\{k\},y\_\{k\}\)that updatesHkH\_\{k\}, so the damage accumulates\. Giovannelli et al\. report that surrogate modeling “hardly” improves their state\-of\-the\-art FLE method when used to approximate gradients, even with NN surrogates\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]\. Without uniform region\-wise error control or acceptance safeguards tied to the true objective, replacement can underperform even where pointwise approximation looks good\.

## 9Numerical Experiments

This section reports the numerical evidence\. The experiments are small, and they hold the network class, the training pipeline and the base method fixed across variants, so a difference between two of them cannot be blamed on implementation\. That is weaker than varying one mechanism at a time, and we do not claim the stronger thing: two embeddings differ in more than their nominal role, since changing it also changes the evaluation cost per iteration, the number of proposals, the fallback behaviour and the trajectory along which the surrogate is trained, and those in turn change when a run effectively stops\. The deterministic runs of Sections[9\.2](https://arxiv.org/html/2608.24963#S9.SS2)to[9\.4](https://arxiv.org/html/2608.24963#S9.SS4)use a single run per instance\. What follows is therefore evidence consistent with the mechanisms we have described and a measurement of how large their effects are in a controlled setting, not a causal isolation of the role variable and not a benchmark against state\-of\-the\-art model\-based DFO software\. The complete implementation, the31463146raw result files behind every number reported here, and the scripts that turn one into the other are available at[https://github\.com/chengkuobian/neural\-surrogates\-dfo](https://github.com/chengkuobian/neural-surrogates-dfo): plain Python with NumPy and Matplotlib only, the surrogate, its analytic parameter gradients, the Adam trainer and all four solvers implemented from scratch, with aREADMEgiving the exact commands, parameter values, caps and fallback rules, and a table mapping each script to the table or figure it produces\.

### 9\.1Setup

##### Problems\.

We use thirteen problems from the Moré–Garbow–Hillstrom collection\[[19](https://arxiv.org/html/2608.24963#bib.bib19)\]: sphere, chained Rosenbrock, Powell singular, trigonometric, Broyden tridiagonal, Broyden banded, variably dimensioned, discrete boundary value, discrete integral equation, Brown almost\-linear, penalty I, extended Freudenstein–Roth and extended Wood\. Each is taken in dimensionsn∈\{4,8,16\}n\\in\\\{4,8,16\\\}from three starting points, the standard one, the standard one scaled by ten, and a random perturbation of it, giving117117instances\. These are117117instances of thirteen functions, not117117independent problems; counts across them are not independent samples, which is why the intervals reported below resample problems rather than instances\. The budget is100100simplex gradients, that is100​\(n\+1\)100\(n\+1\)evaluations, the convention of\[[20](https://arxiv.org/html/2608.24963#bib.bib20)\]\.

##### Methods\.

Six solvers are compared\.

1. 1\.FD\-base: the derivative\-free quadratic regularization method of\[[11](https://arxiv.org/html/2608.24963#bib.bib11)\], in the simplified finite\-difference Armijo\-type form\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]use as their base method, with the inner loop capped for practicality,ϵ=10−3\\epsilon=10^\{\-3\}andσ0=1\\sigma\_\{0\}=1\.
2. 2\.NN\-assist: the safeguarded assistance embedding \(Algorithm[1](https://arxiv.org/html/2608.24963#alg1), driven as in\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\) with a Sobolev\-trained surrogate\.
3. 3\.NN\-assist\-val: identical to NN\-assist except that the gradient term is removed from the loss \([27](https://arxiv.org/html/2608.24963#S8.E27)\) \(value\-only training\)\.
4. 4\.DFO\-TR: a model\-based trust\-region method with minimum\-Frobenius\-norm quadratic interpolation on2​n\+12n\+1points, exact trust\-region subproblem solves, and a geometry\-improving step when the interpolation set is poorly poised\. This is the standard scheme behindnewuoaandbobyqa\[[6](https://arxiv.org/html/2608.24963#bib.bib6)\]and the DFO\-TR family of\[[1](https://arxiv.org/html/2608.24963#bib.bib1)\]; implemented alongside the rest because our compute environment has no network access; its interpolation conditions and subproblem solutions are checked against direct computation, and Section[9\.10](https://arxiv.org/html/2608.24963#S9.SS10)reports a comparison against released Py\-BOBYQA run on the authors’ own machine\.
5. 5\.DFO\-TR\+NN: the same method with the assistance loop attached\. It refits the surrogate once everyn\+1n\+1accepted steps rather than at every step, because a trust\-region iteration costs one evaluation where a finite\-difference outer iteration costsn\+1n\+1; this equalises training cost per evaluation so the comparison is not confounded by refit frequency\.
6. 6\.NN\-replace: surrogate\-gradient replacement with line search and finite\-difference fallback\. Once\|F\|≥2​\(n\+1\)\|F\|\\geq 2\(n\+1\)the FD gradient is*replaced*by∇mθ​\(xk\)\\nabla m\_\{\\theta\}\(x\_\{k\}\), savingnnevaluations per iteration; backtracking still acts on the true objective and repeated failure triggers an FD fallback\. It is not an unsafeguarded method, so the contrast with assistance is not safeguarded versus unsafeguarded: the two also differ in evaluation cost per iteration and in how often the surrogate is consulted\. Keeping the skeleton fixed removes quasi\-Newton curvature updating as a confounder; how a surrogate gradient would behave inside an FD\-BFGS step such as FLE\-S is a separate question that this experiment does not address, and that\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]studies directly\.

The surrogate is a one\-hidden\-layer softplus network of width6​n6nwith inputs and targets normalized per training, weight decayλ=10−4\\lambda=10^\{\-4\}, trained from a fresh He initialization by full\-batch Adam \(900900iterations, the setting Section[9\.8](https://arxiv.org/html/2608.24963#S9.SS8)selects\); datasets are capped atN≤10​\(n\+1\)N\\leq 10\(n\+1\)andM≤10M\\leq 10with oldest\-first removal, following\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\], as areσ0=1\\sigma\_\{0\}=1andλ=10−4\\lambda=10^\{\-4\}\. We depart from it in five settings, all listed here because two of them change how readily a surrogate step is taken: Adam for900900iterations rather than L\-BFGS, width6​n6nrather than5​n5n,ϵ=10−3\\epsilon=10^\{\-3\}rather than10−510^\{\-5\},σmin=10−10\\sigma\_\{\\min\}=10^\{\-10\}rather than10−210^\{\-2\}, and, in the acceptance machinery itself,ρ=0\.25\\rho=0\.25rather than10−410^\{\-4\}andγ=100\\gamma=100rather than25/225/2\. The last two make a surrogate step harder to accept than in the original, so our acceptance rates are lower than that design would give and the advantage we measure for assistance is conservative\. All of them apply identically to every surrogate variant, so none affects the role comparison, which is a difference between variants rather than an absolute level\. All solver runs are deterministic; the region\-size study below uses three random seeds and reports medians\.

##### Metrics\.

The primary metric throughout is the number of true function evaluations, for the reason the problem class supplies: in the simulation\-based setting these experiments are about, one evaluation is a model run and everything else is free by comparison, so sample efficiency is what a method is judged on\. We report \(i\) evaluations to reach a target accuracy, \(ii\) the success rate across problems as data profiles, and \(iii\) the surrogate acceptance rate\.

Wall\-clock time appears in Sections[9\.9](https://arxiv.org/html/2608.24963#S9.SS9)and[10](https://arxiv.org/html/2608.24963#S10)and should be read for what it is\. Every solver here is unoptimised pure NumPy, so absolute times are properties of this implementation rather than of the algorithms and we draw no ranking from them\. What they are used for is a question implementation speed does not affect: at what oracle cost the evaluations a surrogate saves outweigh the time it takes to fit, and how that threshold moves withnnwhen the parameter count grows quadratically and the training set does not\. Probe\-set errors and redundancy measures are likewise reported only where they explain a primary metric\. Accuracy is measured in the standard way\[[20](https://arxiv.org/html/2608.24963#bib.bib20)\]: withxNx\_\{N\}the best point afterNNevaluations,x0x\_\{0\}the start andfLf\_\{L\}the smallest value any solver in the pool attains, a run reaches accuracyτ\\tauat the firstNNwithf⁡\(x0\)−f⁡\(xN\)≥\(1−τ\)​\(f⁡\(x0\)−fL\)f\(x\_\{0\}\)\-f\(x\_\{N\}\)\\geq\(1\-\\tau\)\\bigl\(f\(x\_\{0\}\)\-f\_\{L\}\\bigr\), andNa,p=\+∞N\_\{a,p\}=\+\\inftyif the budget runs out first\. The data profileδa​\(β\)\\delta\_\{a\}\(\\beta\)is the fraction of problems withNa,p≤β⁡\(np\+1\)N\_\{a,p\}\\leq\\beta\(n\_\{p\}\+1\), that is, solved withinβ\\betagroups ofnp\+1n\_\{p\}\+1evaluations\[[20](https://arxiv.org/html/2608.24963#bib.bib20)\]; performance profiles\[[21](https://arxiv.org/html/2608.24963#bib.bib21)\]give the same ranking on this test set and are omitted\.

### 9\.2Role of the Embedding: Replacement versus Assistance

Table[1](https://arxiv.org/html/2608.24963#S9.T1)reports how many of the117117instances each method solves at three accuracy levels, and Figure[1](https://arxiv.org/html/2608.24963#S9.F1)shows the corresponding data profiles\.

Table 1:Extended benchmark,117117instances solved to accuracyτ\\tauwithin100100simplex gradients, withfLf\_\{L\}taken over all eight variants below, released Py\-BOBYQA included\. The last column counts instances on which a method returns the lowest final value of the eight\. Surrogate variants use the width\-6​n6n,900900\-iteration network of Section[9\.8](https://arxiv.org/html/2608.24963#S9.SS8); the last row repeats safeguarded assistance at the width\-3​n3n,300300\-iteration setting to show what undertraining costs\.Methodτ=10−1\\tau=10^\{\-1\}10−310^\{\-3\}10−510^\{\-5\}lowestffFD base method10993678NN assist \(Sobolev\)113100847NN assist \(value\-only\)10994661NN replace110100650DFO\-TR \(ours\)1131028833DFO\-TR\+\+NN assist1131038610Py\-BOBYQA \(released\)11711610355NN assist, width3​n3n,300300it\.110100753Figure 1:Data profilesδa​\(β\)\\delta\_\{a\}\(\\beta\)over the117117instances, withfLf\_\{L\}taken over all six methods, all surrogate variants at the width\-6​n6n,900900\-iteration setting\. Safeguarded assistance separates from its own base method as the tolerance tightens and tracks the model\-based trust\-region solver, gradient replacement is the weakest embedding, and attaching the surrogate to the trust\-region solver does not help\.Table[1](https://arxiv.org/html/2608.24963#S9.T1)scores every variant against released Py\-BOBYQA\[[6](https://arxiv.org/html/2608.24963#bib.bib6)\]rather than our own interpolation code, and the choice matters\. The Moré–Wild criterion asks a solver to reachf0−τ⁡\(f0−fL\)f\_\{0\}\-\\tau\(f\_\{0\}\-f\_\{L\}\), so a released solver ending two orders of magnitude lower movesfLf\_\{L\}enough to reorder the methods above it: scored among the surrogate variants alone, assistance leads our trust\-region method9999to9494; adding the released solver reverses this to8484against8888\. We report the second, becausefLf\_\{L\}should be the best value anyone attains, and record the first because it is the comparison a paper that never installs the released code would make\.

What the surrogate does for the method it is attached to is not in doubt at either reference\. Safeguarded assistance lifts its own base method from6767instances to8484atτ=10−5\\tau=10^\{\-5\}, from9393to100100at10−310^\{\-3\}and from109109to113113at10−110^\{\-1\}: the margin is largest exactly where the base method has stalled and least where it succeeds anyway, which is the smoothing argument of Section[3\.1](https://arxiv.org/html/2608.24963#S3.SS1)and the reason a stalled finite\-difference method still has a dataset worth fitting\. Seventeen instances of117117is the effect the rest of this paper is about\.

The last row of the table is the same embedding at the width\-3​n3n,300300\-iteration network an earlier version of this paper used, and it matters for how the rest should be read: that setting reaches7575where the trained one reaches8484\. A third of what looked like a limit of the mechanism was a limit of the fit\. Every surrogate number reported from here on uses the trained setting\.

##### How much of that is the base method’s weakness?

Almost all of it\. Both interpolation methods are ahead of assistance at every tolerance, our own by four instances atτ=10−5\\tau=10^\{\-5\}and the released one by nineteen\. The last column sharpens it: Py\-BOBYQA returns the lowest of the eight final values on5555instances and our trust\-region method on3333, against77for the surrogate\. The gap in final values is larger than the gap in counts: the pairedlog10\\log\_\{10\}ratio between Py\-BOBYQA and the surrogate has median−2\.321\-2\.321with interval\[−5\.189,−0\.168\]\[\-5\.189,\-0\.168\], so the released solver typically ends more than two orders of magnitude lower, and it returns the lowest value of the four on5757of the117117instances against the surrogate’s88\. A trained surrogate on a weak base method is a large improvement on that base method and is not a substitute for a well\-engineered interpolation model\.

Nor does it help one\. Attaching the same trained network to our trust\-region solver gives8686against the8888it reaches alone, which is the room condition of Section[9\.5](https://arxiv.org/html/2608.24963#S9.SS5)seen from the other side: a method that already builds a quadratic from the same evaluations has nothing left for a surrogate to claim\.

##### One problem where it does win\.

The aggregate hides a case worth isolating, because it is the only one in this set where a learned local model beats a well\-engineered interpolation model outright\. On extended Wood the surrogate solves77of its99instances atτ=10−5\\tau=10^\{\-5\}against44for Py\-BOBYQA,44for our trust\-region method and44for its own base: it is the only one of the eight variants that gets past four, and it returns the lowest final value of all eight on33of the nine, once by2\.32\.3orders of magnitude\. Extended Wood is a coupled quartic whose minimum sits at the end of a narrow curved valley, the regime in which a minimum\-Frobenius\-norm quadratic is a poor model at any radius the method works at, and it is exactly the regime the smoothing argument of Section[3\.1](https://arxiv.org/html/2608.24963#S3.SS1)points to\.

We resist generalising, for two reasons the same data supply\. Chained Rosenbrock and extended Freudenstein–Roth have curved valleys of the same kind and the surrogate is behind on both\. And the approximation\-level statistic that ought to predict the win does not: the ratio of the quadratic’s uniform gradient error to the network’s is0\.890\.89atn=16n=16, the instance won by the largest margin, so by that measure the quadratic is the better local model exactly where the surrogate wins\. One problem is an existence result, not a rule\. The set of objectives on which a safeguarded neural surrogate beats released interpolation software is not empty, and finding a diagnostic that identifies it in advance is the question this paper leaves in the best shape for someone to answer\.

##### Does the picture change with dimension?

The obvious reply is thatn≤16n\\leq 16is where interpolation is at its strongest, since a minimum\-Frobenius\-norm quadratic needs only2​n\+1≤332n\+1\\leq 33points, and that a learned model should come into its own where that number grows\. This deserves a direct test rather than a limitation note, so we repeated the comparison atn=32n=32,n=64n=64andn=128n=128, where a quadratic needs6565,129129and257257points, on all thirteen problems under a budget of50​\(n\+1\)50\(n\+1\), and added released Py\-BOBYQA\[[6](https://arxiv.org/html/2608.24963#bib.bib6)\]at the two largest so that the comparison does not rest on our own interpolation code\.

Adding it changed the answer, and the episode bears on how such comparisons should be read\. Our own trust\-region implementation degrades with dimension: atτ=10−5\\tau=10^\{\-5\}it leads safeguarded assistance8888to8484atn≤16n\\leq 16, leads1111to88atn=32n=32, and atn=64n=64solves55of1313against77for assistance\. Each of these pairs is scored against the solvers available at that dimension, eight atn≤16n\\leq 16, four atn=32n=32andn=64n=64and four atn=128n=128, so a pair may be compared internally but the counts should not be read across dimensions as a single series\. Its last improvement on the median instance atn=64n=64comes at evaluation20242024of32503250against29772977\. Taken alone this looks like the crossover the motivation for learned surrogates predicts, arriving where the bookkeeping of a129129\-point interpolation set begins to dominate the budget\.

Released Py\-BOBYQA shows the reading is wrong: on the same thirteen problems atn=64n=64it solves1212of1313, far ahead of assistance at77and of our implementation at55, and1212against99and88at the looserτ=10−3\\tau=10^\{\-3\}, so its lead is not an artefact of demanding high accuracy\. The degradation is a property of our code, not of interpolation, since a solver that maintains its model by rank\-one updates and manages geometry carefully loses nothing at this dimension\. Our implementation is an adequate stand\-in atn≤16n\\leq 16, where its medianlog10\\log\_\{10\}difference from the released solver is0\.0000\.000\(Section[9\.10](https://arxiv.org/html/2608.24963#S9.SS10)\), and is not one atn=64n=64\. We report this rather than quietly dropping the earlier claim, because the same trap awaits anyone who benchmarks against a reimplementation without checking it at the dimension of interest\.

Atn=128n=128, where a quadratic needs257257points and the budget allows64506450evaluations, the released solver solves1010of1313atτ=10−5\\tau=10^\{\-5\}and returns the lowest value on88; assistance solves88and returns the lowest on33; the finite\-difference base method solves66\. Training the network properly is worth two instances here, the same embedding at the width\-3​n3nsetting solving66and never returning the best value, but it does not close the gap\. The pattern across the four dimensions is therefore not that the surrogate’s prospects improve as2​n\+12n\+1grows\. It is that a trained surrogate falls behind a well\-engineered interpolation method at every dimension, by nineteen instances of117117atn≤16n\\leq 16and by two of thirteen atn=128n=128\.

The conclusion is uniform across every dimension we can test\. A trained surrogate on a weak base method is not a weak method: it lifts that base by seventeen instances of117117, and on extended Wood it beats released interpolation software outright\. What it does not do is match that software in aggregate, descend as far on most instances, or add anything to a method already building a quadratic from the same evaluations\.

Second, replacement is behind the base method where the tolerance separates them,6565against6767atτ=10−5\\tau=10^\{\-5\}, and it is behind at the trained setting, so the loss is not a training artefact\. It also returns the lowest of the eight final values on no instance at all\. Across the benchmark it accepted no surrogate\-gradient step that survived the true\-objective test and triggered678678finite\-difference fallbacks\. The fallback count is a direct measurement and the natural explanation for the loss, though it is an observed correlate rather than an independently manipulated variable\.

Third, at the loose tolerance the methods are within five instances of each other\. Safeguarding makes assistance close to free: a rejected proposal costs one evaluation, so on instances where the base method already succeeds quickly there is nothing to gain and little to lose\.

### 9\.3Region Size and the Generalization Radius

To probe the radius half of the claim outside the optimization loop, we sampleN=10​\(n\+1\)N=10\(n\+1\)training points uniformly inℬ⁡\(x0,Δ\)\\mathcal\{B\}\\\!\\left\(x\_\{0\};\\Delta\\right\), add forward\-difference gradient targets atM=10M=10of them, train surrogates with the Sobolev and the value\-only loss, and measure the uniform errorsef,ege\_\{f\},e\_\{g\}of Definition[4](https://arxiv.org/html/2608.24963#Thmdefinition4)on300300independent probe points in the same ball\. The sweep covers all thirteen problems atn∈\{4,8,16\}n\\in\\\{4,8,16\\\}and seven radii from10−210^\{\-2\}to3\.23\.2with five seeds each,273273cells in all; reference gradients are central differences ath=10−6h=10^\{\-6\}, accurate to about10−910^\{\-9\}here\. Since the problems differ by orders of magnitude in scale we reportege\_\{g\}relative to the median‖∇f‖\\left\\lVert\\nabla f\\right\\rVertover the probe points, which makes the cells comparable, and intervals resample problems\.

Figure 2:Region\-size sweep,273273cells\. \(a\) Uniform gradient error relative to the median gradient norm in the same ball, median over cells with a bootstrap band over problems; the dotted line is the slope a fully linear model would have\. \(b\) Distribution over the3939problem–dimension pairs of the fitted exponentβ\\betaineg∼Δβe\_\{g\}\\sim\\Delta^\{\\beta\};69%69\\%fall below the fully linear rateβ=1\\beta=1, and those above it still carry a nonzero floor at smallΔ\\Delta\.Figure[2](https://arxiv.org/html/2608.24963#S9.F2)shows the two effects predicted by Section[7](https://arxiv.org/html/2608.24963#S7)\. The relative gradient error rises from0\.0550\.055, with interval\[0\.050,0\.063\]\[0\.050,0\.063\], atΔ=10−2\\Delta=10^\{\-2\}to1\.571\.57, with interval\[0\.89,2\.60\]\[0\.89,2\.60\], atΔ=3\.2\\Delta=3\.2: at the largest radius the surrogate’s gradient error exceeds the gradient it is estimating, which by Proposition[2](https://arxiv.org/html/2608.24963#Thmproposition2)is exactly where descent can no longer be guaranteed\. With a fixed sample budget the fill distanceε\\varepsilonscales withΔ\\Deltaand the coverage term of Proposition[1](https://arxiv.org/html/2608.24963#Thmproposition1)takes over, so the surrogate has a finite effective generalization radius \(Definition[5](https://arxiv.org/html/2608.24963#Thmdefinition5)\)\.

The second effect is the more consequential, and the larger sweep sharpens it\. Shrinking the radius does not drive the error to zero\. Over the smallest three radii, a factor of seven inΔ\\Delta, the relative error moves only from0\.0550\.055to0\.0750\.075; fittingeg∼Δβe\_\{g\}\\sim\\Delta^\{\\beta\}on each problem–dimension pair gives a median exponent ofβ=0\.87\\beta=0\.87with interquartile range\[0\.59,1\.01\]\[0\.59,1\.01\], and69%69\\%of the3939pairs grow more slowly than linearly\. A fully linear model requireseg=O⁡\(Δ\)e\_\{g\}=O\(\\Delta\), that isβ≥1\\beta\\geq 1with the constant carrying no floor; what we observe instead is the training\-error floorϵg\\epsilon\_\{g\}of Proposition[1](https://arxiv.org/html/2608.24963#Thmproposition1)dominating at small radii\. An NN surrogate trained on an optimization\-path\-sized dataset is therefore*not*automatically fully linear at small radii, precisely the regime a replacement embedding needs it to be, consistent with the caution of\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]\. Sobolev training reduces the relative error by a factor of1\.161\.16to1\.661\.66across radii, most at moderateΔ\\Delta, but does not move either the small\-radius floor or the large\-Δ\\Deltawall\.

##### Does depth move the floor?

The architecture is fixed at one hidden layer throughout, so the floor could in principle be a property of that choice\. Repeating the sweep with two hidden layers, the depth\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]use, at the same width as our single\-layer network, over the same273273cells, moves it the wrong way: the ratio of relative gradient errors, two layers over one, is1\.0691\.069with interval\[1\.04,1\.15\]\[1\.04,1\.15\]atΔ=10−2\\Delta=10^\{\-2\}and1\.0421\.042with\[1\.01,1\.12\]\[1\.01,1\.12\]atΔ=2\.6×10−2\\Delta=2\.6\\times 10^\{\-2\}, both excluding one, and it falls to0\.9630\.963with\[0\.93,0\.99\]\[0\.93,0\.99\]only atΔ≈1\.2\\Delta\\approx 1\.2\. Depth buys a few percent at radii where the error already exceeds the gradient it is estimating and costs a few percent exactly where full linearity would be needed\. Both networks are trained value\-only here, since the Sobolev loss’s analytic parameter gradients are derived for one hidden layer and the second\-order terms a second layer introduces are a source of silent error we chose not to risk; Section[9\.4](https://arxiv.org/html/2608.24963#S9.SS4)measures what dropping the gradient term costs, so the two effects can be read together\. Within that qualification, the small\-radius floor is not an artefact of shallowness\.

That comparison holds the sample set fixed, so a second layer could still pay once the surrogate steers the iteration\. Run inside the assisted method on all117117instances it does not: atτ=10−3\\tau=10^\{\-3\}the depths are indistinguishable,111111instances against112112, and atτ=10−5\\tau=10^\{\-5\}, where the floor binds, the deeper network solves9797against101101and ends lower on2424instances\. Acceptance barely moves,0\.1990\.199\(591/2967591/2967\) against0\.2040\.204\(722/3532722/3532\), so the extra layer proposes more steps at the same success rate rather than steps the test rejects less often\. Depth is not the missing ingredient at either level, the conclusion Section[9\.8](https://arxiv.org/html/2608.24963#S9.SS8)reaches for width and training length\.

### 9\.4Sobolev versus Value\-Only Training

The loss ablation connects the two previous experiments to the training signal\. Across the117117instances the surrogate acceptance rate falls from0\.7030\.703\(7488/106507488/10650\) with the Sobolev loss to0\.1480\.148\(496/3354496/3354\) with value\-only training, and the gap is not something more training closes: both arms use the width\-6​n6n,900900\-iteration network\. The effect on the optimisation is as large as the effect of the embedding role\. Value\-only assistance solves6666instances atτ=10−5\\tau=10^\{\-5\}against8484for Sobolev assistance, below the base method’s6767rather than anywhere between the two\. Because the two arms differ in one term of the loss and nothing else, the comparison can be made instance by instance rather than through a pooled reference set: the median oflog10\\log\_\{10\}of the ratio of final objective values, value\-only over Sobolev, is\+0\.099\+0\.099with bootstrap interval\[\+0\.019,\+0\.214\]\[\+0\.019,\+0\.214\], and Sobolev is lower on7777of the117117instances against2525the other way,1515tied \(p=2\.5×10−7p=2\.5\\times 10^\{\-7\}, sign test\)\. The typical instance gains about a fifth of a decade, which is modest; the count is what the claim rests on\. This matches Proposition[3](https://arxiv.org/html/2608.24963#Thmproposition3)and the trend\-prediction discussion of Section[3\.1](https://arxiv.org/html/2608.24963#S3.SS1): without gradient information the surrogate may fit the sampled values well and still fail to produce descent directions that survive the true\-objective test\. On the mildly nonlinear trigonometric and Broyden problems value\-only assistance does still help, the safeguard converting even a mediocre surrogate into a no\-worse\-than\-base option\.

### 9\.5Noise, and the Room the Base Method Leaves

Simulation output is random, so the experiments above miss the feature that matters most in practice\. We repeat the comparison with an oracle returningf​\(x\)​\(1\+σ​z\)f\(x\)\\,\(1\+\\sigma z\), wherez∼𝒩⁡\(0,1\)z\\sim\\mathcal\{N\}\(0,1\)is drawn afresh at every call: the relative\-noise model used in the noisy variants of the Moré–Wild collection\[[20](https://arxiv.org/html/2608.24963#bib.bib20)\]\. Progress is scored on the*true*ffat the best point returned, so a solver gains nothing from a lucky read\. The sweep covers ten problems atn∈\{4,8\}n\\in\\\{4,8\\\}and ten noise levels fromσ=0\\sigma=0toσ=10−2\\sigma=10^\{\-2\}, with five paired seeds per cell,200200cells,10001000paired comparisons and20002000solver runs in all; at a given seed both methods meet the same noise stream, so every comparison is paired\. We useϵ=10−1\\epsilon=10^\{\-1\}and a budget of200​\(n\+1\)200\(n\+1\)evaluations; at theϵ=10−3\\epsilon=10^\{\-3\}of the deterministic study every solver, the base method included, stalls immediately at any nonzero noise, which is itself part of what follows\.

Figure 3:Noise sweep,200200cells\. \(a\) The fraction of problem–dimension pairs on which safeguarded assistance still ends at least1\.26×1\.26\\timeslower, with a bootstrap band over problems; the decline is real but the band is wide throughout\. \(b\) On the same runs, the fraction of the evaluation budget consumed before the run stops, and the acceptance rate of surrogate proposals\. The runs get shorter; the proposals do not get worse\.##### There is no threshold in the noise\.

Figure[3](https://arxiv.org/html/2608.24963#S9.F3)\(a\) shows the fraction of the twenty problem–dimension pairs that are helped falling from0\.700\.70atσ=0\\sigma=0and0\.750\.75atσ=10−6\\sigma=10^\{\-6\}to0\.350\.35atσ=10−4\\sigma=10^\{\-4\}and0\.150\.15atσ=10−2\\sigma=10^\{\-2\}\. The bootstrap intervals overlap throughout, and the per\-problem picture explains why: the largest noise level at which a pair is still helped ranges over four decades, from10−610^\{\-6\}on extended Freudenstein–Roth atn=4n=4to10−210^\{\-2\}on the trigonometric function and the discrete integral equation, with a median of3×10−53\\times 10^\{\-5\}; three pairs are helped only in the noiseless case and one, variably dimensioned atn=8n=8, at no level at all\. Where a pair is helped the effect is substantial, the median gain among helped cells being0\.210\.21inlog10\\log\_\{10\}with an interquartile range of\[0\.08,0\.65\]\[0\.08,0\.65\], but which pairs those are depends on the problem far more than on the noise\. A clean cliff atσ≈3×10−5\\sigma\\approx 3\\times 10^\{\-5\}appears if one looks at two problems only; it is the threshold of the Powell singular function, which the wider sweep reproduces exactly, and not a property of the method\.

##### What does change is the base method\.

Figure[3](https://arxiv.org/html/2608.24963#S9.F3)\(b\) locates the effect elsewhere\. The fraction of the evaluation budget consumed before a run stops falls monotonically from62\.3%62\.3\\%atσ=0\\sigma=0to21\.7%21\.7\\%atσ=10−2\\sigma=10^\{\-2\}, and the base method and the assistance method track each other to within three percentage points at every level: they stop at the same place\. Surrogate proposals across the sweep fall from1868218682to152152, a factor of123123\. Their acceptance rate falls too, from0\.8430\.843atσ=0\\sigma=0to0\.5200\.520atσ=10−2\\sigma=10^\{\-2\}, but by a factor of1\.61\.6against the proposals’123123\. Two orders of magnitude separate the two effects, and only the smaller of them is about the surrogate: proposals disappear because the run has already stopped, not because they are being rejected\.

The cause is the base method’s Armijo test, which compares two noisy evaluations and, once the noise exceeds the decrease being tested for, stops accepting anything; the inner loop then shrinks the finite\-difference stepsize indefinitely and the run stalls\. Every surrogate proposal is validated by a test of exactly this kind, so a safeguard inherits the noise tolerance of the test it relies on and cannot exceed it\. Unsafeguarded replacement is actively harmful: atσ=10−6\\sigma=10^\{\-6\}on Powell it returns a median of5\.45\.4against0\.0660\.066for the base method\.

##### What would make a surrogate pay at this noise level?

None of it is a matter of training the network better\. The binding quantity is the variance entering the acceptance test, which the surrogate inherits twice, through gradient targets that are finite differences of noisy runs and through the test that validates its proposals\. Gradient targets would have to be estimated rather than differenced, by common random numbers, batch regression, or an unbiased estimator where the model admits one; the acceptance test would have to become a statistical comparison at a sample size the method controls; and the surrogate would have to be trained on replication\-averaged targets of known noise level\. Each is a change to the interface between model and oracle rather than to the model, which is what the ablation below finds from the other side\.

##### The third condition\.

This qualifies the role and radius conditions rather than replacing them\. A surrogate accelerates a base method that is still making progress and cannot accelerate one that has stopped, so the budget the base method leaves unspent is the room available to any acceleration; under noise that room closes well before the surrogate becomes the limiting factor\. The consequence is an ordering: make the acceptance test survive the noise first, through replication, common random numbers or a statistical test in place of a deterministic sufficient\-decrease inequality, and only then ask what a surrogate adds\. Section[9\.6](https://arxiv.org/html/2608.24963#S9.SS6)puts that ordering to the test on a real simulation model\.

### 9\.6A Simulation Testbed, and What Breaks on It

Section[9\.5](https://arxiv.org/html/2608.24963#S9.SS5)adds noise to analytic functions, which isolates a mechanism but is not a simulation\. The ordering it suggests — repair the acceptance test first, then ask about the surrogate — is a claim about practice, and we test it here on a genuine Monte\-Carlo model\.

##### The model\.

Four products share a replenishment channel under a periodic\-review\(s,S\)\(s,S\)policy, with Poisson demand of meanλ=\(5,10,15,20\)\\lambda=\(5,10,15,20\): when the inventory position of productjjfalls tosjs\_\{j\}or below it is raised toSjS\_\{j\}\. A period in which at least one product orders incurs a joint setup costK=60K=60, and each ordering product incurs a minor setup cost of88, a unit purchase costcjc\_\{j\}, a holding costhjh\_\{j\}and a backorder penaltypjp\_\{j\}\. The joint setup couples the eight decision variablesx=\(s1,d1,…,s4,d4\)x=\(s\_\{1\},d\_\{1\},\\dots,s\_\{4\},d\_\{4\}\), whereSj=sj\+max⁡\(dj,0\)S\_\{j\}=s\_\{j\}\+\\max\(d\_\{j\},0\), so the problem is not separable\. One oracle call runs one replication over6060periods and returns the realised average cost; the single\-replication standard deviation at the starting design is3\.53\.5against a mean of222\.3222\.3, a relative noise of1\.6%1\.6\\%\. Reported costs are expected costs of the returned design, estimated separately under common random numbers, so no solver gains from a lucky read\. A random search under the same20002000\-replication budget reaches206\.8206\.8; we use it as a rough indication of what is reachable, not as an optimum\.

##### The methods of Sections[9\.2](https://arxiv.org/html/2608.24963#S9.SS2)–[9\.5](https://arxiv.org/html/2608.24963#S9.SS5)do not work here\.

The comparison uses3030seeds at each of three starting designs, costing222\.3222\.3,239\.3239\.3and226\.1226\.1, for9090runs per method with a budget of20002000replications each\. A seed fixes the demand stream, so every method meets the same randomness at the same start; comparisons are paired by construction and we report medians of paired differences with bootstrap95%95\\%confidence intervals\. With a deterministic oracle interface the finite\-difference base method returns its starting design on the median run, and so do both surrogate variants: the median paired difference against the base method is0\.000\.00for safeguarded assistance and\+2\.32\+2\.32for replacement, neither significant\. The failure is complete rather than marginal, and it does not originate in the surrogate\.

The diagnosis follows from the stepsize rulehi=2​ϵ/\(5​n​2i​σk\)h\_\{i\}=2\\epsilon/\(5\\sqrt\{n\}\\,2^\{i\}\\sigma\_\{k\}\), which ties the difference stepsize to the backtracking index\. That coupling is right in a deterministic setting, where shrinkinghhwith the step reduces truncation error, but under a stochastic oracle the noise in a difference quotient scales ash−1h^\{\-1\}, so each backtrack*doubles*the gradient noise\. At the starting designh0≈0\.14h\_\{0\}\\approx 0\.14puts the noise on each gradient component near3535while the true partial derivatives are of order one, and the computed “gradient” is a random vector of norm about100100\. The Armijo test then asks for a decrease of18​α​‖g‖2≈1250\\tfrac\{1\}\{8\}\\alpha\\\|g\\\|^\{2\}\\approx 1250from a cost of222222, which is unattainable, so the method backtracks, worsens its own gradient estimate, raises its own acceptance threshold and repeats until the inner loop is exhausted\. Raising the stepsize multiplier does not repair this, since the required decrease grows with‖g‖2\\\|g\\\|^\{2\}whileϵ\\epsiloncontrols the stepsize, the stationarity test and the acceptance threshold at once\. This is the mechanism of Section[9\.5](https://arxiv.org/html/2608.24963#S9.SS5)at a noise level,1\.6×10−21\.6\\times 10^\{\-2\}, far past where the synthetic sweep shows runs stalling\.

##### A noise\-aware interface\.

We therefore modify the oracle interface, leaving the surrogate machinery untouched: the difference stepsize is decoupled from the backtracking index and set from an estimateσ^ε\\hat\{\\sigma\}\_\{\\varepsilon\}of the replication standard deviation, balancing truncation against a noise term of orderσ^ε/h\\hat\{\\sigma\}\_\{\\varepsilon\}/h; every evaluation entering a gradient averagesrrreplications, so gradient noise falls asr−1/2r^\{\-1/2\}; and acceptance requires the observed decrease to exceed a standard Armijo term and, optionally,zzstandard errors of it,

f^​\(x\)−f^​\(x\+\)≥c​α​‖g‖2\+z​SE^,SE^=σ^ε​2/r,\\hat\{f\}\(x\)\-\\hat\{f\}\(x^\{\+\}\)\\ \\geq\\ c\\,\\alpha\\\|g\\\|^\{2\}\\;\+\\;z\\,\\widehat\{\\mathrm\{SE\}\},\\qquad\\widehat\{\\mathrm\{SE\}\}=\\hat\{\\sigma\}\_\{\\varepsilon\}\\sqrt\{2/r\},\(35\)withc=10−4c=10^\{\-4\},z=1z=1andr=10r=10\. Nothing here is new; the same ingredients underpin adaptive\-sampling trust\-region methods for simulation optimisation, of which ASTRO\-DF\[[18](https://arxiv.org/html/2608.24963#bib.bib18)\]is the reference example, and our point is precisely that they are*necessary*before a question about surrogates can even be posed\.

##### Result\.

Table[2](https://arxiv.org/html/2608.24963#S9.T2)reports the comparison\. Repairing the interface is worth a median paired improvement of10\.4010\.40cost units, with bootstrap interval\[8\.84,12\.74\]\[8\.84,12\.74\], the repaired method being ahead on7373of9090runs\. Adding safeguarded surrogate assistance on top of the repair is worth nothing measurable: the median paired difference is0\.000\.00with interval\[0\.00,0\.18\]\[0\.00,0\.18\], and the surrogate variant leads on only3333of9090runs\.

A five\-seed pilot had suggested a gain of about0\.70\.7cost units; it does not survive atN=90N=90, and we report the correction\. The repaired base method reaches211\.6211\.6against the206\.8206\.8a random search of the same budget reaches, capturing roughly nine tenths of the achievable reduction on its own, so little is left for any acceleration to claim\. This is not a budget artefact: raising the budget from20002000to80008000replications changes nothing, because both methods stop when backtracking fails rather than when the budget runs out, and at budgets from200200to12001200the paired difference never favours assistance either\. The ordering the rest of the paper is about survives in its weaker form: assistance does its base method no systematic damage at either interface, its median paired difference being zero on the deterministic interface and statistically indistinguishable from zero on the noise\-aware one, though individual runs go both ways, whereas replacement trails the base method on8585of9090runs\.

##### Why the surrogate contributes nothing: it is almost never invoked\.

A null result of this kind is uninformative unless one knows what produced it, so we instrumented the surrogate loop\. Across thirty runs on the four\-product model it issued3131proposals in all, against1065010650over the117117deterministic instances, and on thirteen of the thirty runs it issued none\. Five were accepted, a rate of0\.1610\.161against0\.7030\.703on the deterministic side\. The surrogate is not adding a small benefit our sample size cannot resolve; it is barely being consulted, because the base method stops before it accumulates the data a surrogate needs\.

This is a stronger statement than a null effect, and it has two consequences\. The first concerns the mechanism\. The surrogate’s value targets are clean, being averages of ten replications, but its gradient targets are finite\-difference gradients formed from those same noisy runs, and the diagnosis above shows how badly differencing amplifies replication noise\. Section[9\.4](https://arxiv.org/html/2608.24963#S9.SS4)established on the deterministic side that acceptance depends on gradient information in the loss, the rate falling from0\.7030\.703to0\.1480\.148when it is removed; under simulation noise the gradient targets degrade further and the rate falls to0\.1610\.161\. The three measurements line up: the surrogate loop needs trustworthy gradient information, and simulation\-grade noise degrades it\.

The second concerns the safeguard\. A surrogate that contributes nothing measurable cost the method3131evaluations out of60,00060\{,\}000replications, one part in two thousand, which is what the role distinction of Section[2](https://arxiv.org/html/2608.24963#S2)predicts: under assistance a useless model costs one evaluation per rejected proposal, whereas under replacement a bad model cost entire sequences of failed backtracks \(735735fallbacks across the benchmark\)\. The inventory model is a limiting case of the thesis rather than evidence against it\.

##### Is the null result only a ceiling?

One explanation for the null result is that the repaired base method already reaches211\.6211\.6against the206\.8206\.8of a random search, some2\.3%2\.3\\%away, so a surrogate may simply have had nothing left to win\. We tested that explanation rather than assuming it\. Enlarging the model to eight products,n=16n=16, leaves8\.3%8\.3\\%of the starting cost available, four times the headroom, and the repaired base method there returns a median cost of421\.4421\.4against a random\-search best of416\.9416\.9\. Over9090paired runs the surrogate again adds nothing: the median difference is\+0\.00\+0\.00, with bootstrap interval\[−0\.25,\+0\.91\]\[\-0\.25,\+0\.91\], and assistance leads on3737of9090runs\. The null result is therefore not an artefact of a saturated small model\. Whatever governs it is the noise, not the amount of room\.

Taken with Sections[9\.2](https://arxiv.org/html/2608.24963#S9.SS2)and[9\.5](https://arxiv.org/html/2608.24963#S9.SS5), a sequence emerges\. On the deterministic benchmark, where the base method spends its whole budget, the surrogate lifts it from6767instances to8484while a released interpolation solver reaches103103without one\. Across the noise sweep the base method’s runs shorten to a fifth of their budget and the surrogate’s contribution shortens with them\. On two Monte\-Carlo models the base method does not move at all until its oracle interface is repaired, that repair is worth10\.4010\.40cost units, and the surrogate on top is worth0\.000\.00on both\. What moves along the sequence is the per\-step decrease relative to the noise in the test that accepts it, which suggests an ordering of effort: fix the sampling and acceptance machinery first, and only then ask whether a surrogate is worth its cost\.

Table 2:Joint\-replenishment\(s,S\)\(s,S\)system: expected cost of the returned design over9090runs \(3030seeds×\\times33starting designs\)\. Paired differences are taken against the FD base method at the same seed and start, with bootstrap95%95\\%intervals for the median; “wins” counts runs in which the method beats that base\.MethodmedianIQRpaired diff\. vs FD basewinsFD base method222\.3\[219\.3,226\.1\]\[219\.3,226\.1\]——NN assist \(Sobolev\)222\.3\[218\.4,226\.1\]\[218\.4,226\.1\]\+0\.00\+0\.00\[\+0\.00,\+0\.00\]\[\+0\.00,\+0\.00\]21/90NN replace226\.1\[222\.3,239\.3\]\[222\.3,239\.3\]\+2\.32\+2\.32\[\+0\.00,\+5\.69\]\[\+0\.00,\+5\.69\]5/90Noise\-aware base211\.6\[209\.2,215\.0\]\[209\.2,215\.0\]−10\.40\\mathbf\{\-10\.40\}\[−12\.74,−8\.84\]\\mathbf\{\[\-12\.74,\-8\.84\]\}73/90Noise\-aware\+\+NN assist212\.2\[208\.6,217\.1\]\[208\.6,217\.1\]−10\.47\-10\.47\[−12\.24,−8\.28\]\[\-12\.24,\-8\.28\]71/90Adaptive\-sampling TR211\.4\[210\.6,212\.4\]\[210\.6,212\.4\]——
##### How far is this from the state of the art?

We implemented the mechanism that defines the reference class, adaptive sampling: replications are increased until the standard error falls belowκ​Δ2\\kappa\\Delta^\{2\}, with model, step and ratio test from the trust\-region machinery of Section[9\.2](https://arxiv.org/html/2608.24963#S9.SS2)\. This is the idea behind ASTRO\-DF\[[18](https://arxiv.org/html/2608.24963#bib.bib18)\]; that solver is not installable here, so we run a reproduction of its sampling rule rather than the published code\. Over the same9090runs it returns a median of211\.4211\.4, interquartile range\[210\.6,212\.4\]\[210\.6,212\.4\], against211\.6211\.6for our noise\-aware base method: indistinguishable at this sample size\. The method the surrogate fails to improve is therefore level with an adaptive\-sampling reference rather than a strawman, and all three sit within five cost units of the206\.8206\.8a random search attains\.

##### Which part of the repair does the work?

The interface changes three things at once, so we ablate them\. The result is not the one we expected, and only one of the three components survives the larger sample\. Replication is decisive: weakening it to a single replication per evaluation costs12\.7212\.72cost units, with interval\[9\.49,15\.07\]\[9\.49,15\.07\], and even three replications cost4\.764\.76, with interval\[1\.84,6\.81\]\[1\.84,6\.81\]\. Neither of the other two changes has a detectable effect\. Recoupling the difference stepsize to the backtracking index costs0\.290\.29with interval\[0\.00,1\.28\]\[0\.00,1\.28\], and removing the standard\-error guard of \([35](https://arxiv.org/html/2608.24963#S9.E35)\) is likewise indistinguishable from keeping it\. We had expected that guard to be the essential ingredient\. It is not\.

The surviving reading is therefore sharper and narrower than the one we started with\. What binds is the*variance of the estimates entering the test*, not the form of the test and not the stepsize schedule: ten replications shrink that variance enough for an ordinary Armijo comparison to become reliable, while a more careful rule layered on noisy estimates buys nothing\. This matches the design of adaptive\-sampling trust\-region methods\[[18](https://arxiv.org/html/2608.24963#bib.bib18)\], where the algorithm controls the sample size rather than the acceptance inequality, and it revises Section[9\.5](https://arxiv.org/html/2608.24963#S9.SS5): a safeguard is limited by its acceptance test, but the way to lift that limit is to reduce the noise in what the test compares, not to make the test cleverer\.

##### Caveats\.

We use the simplified base method of\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]withBk=0B\_\{k\}=0, not the quasi\-Newton version\. Practitioners of simulation optimisation would not apply single\-replication differencing in the first place, so the negative result should be read as a statement about transplanting a published deterministic method onto a stochastic oracle unchanged, not as a claim that the method class is unusable\. The noise\-aware variant is a minimal repair rather than a competitive solver; a proper comparison against adaptive\-sampling methods such as ASTRO\-DF is future work\.

### 9\.7Do the Redundancy Diagnostics Predict Anything?

Section[4\.3](https://arxiv.org/html/2608.24963#S4.SS3)defines two diagnostics, the overlap ratiorkr\_\{k\}and the step\-spacing statisticNspN\_\{\\mathrm\{sp\}\}\. Both are natural candidates for deciding when a surrogate should be trusted, and neither had been validated against data\. We do that here, on the same117117instances as Section[9\.2](https://arxiv.org/html/2608.24963#S9.SS2)rather than on a subset of them: at every outer iteration of the assistance solver we log the two diagnostics, the numerical rank of the normalized displacements of the training inputs, and four further candidates computed from quantities the method already holds, and correlate each with the fraction of surrogate proposals accepted at that iteration\. This gives32023202outer iterations across the116116instances on which the solver completes at least one; the remaining instance exhausts its budget inside the first finite\-difference gradient\.

Two methodological choices affect the result\. Outer iterations are clustered within instances, so intervals computed as though the32023202rows were independent are far too narrow; every interval in Table[3](https://arxiv.org/html/2608.24963#S9.T3)comes from a bootstrap that resamples*instances*, which is the level at which the design varies\. And several candidates are heavily tied:NspN\_\{\\mathrm\{sp\}\}equals11at77%77\\%of all iterations\. Rank correlations must therefore use mid\-ranks, since ordinal ranking breaks ties by position in the log and so manufactures correlation out of the order in which iterations were written\. Ordinal ranking makesNspN\_\{\\mathrm\{sp\}\}look predictive for precisely that reason\.

Table 3:Candidate gate statistics against the surrogate acceptance rate:32023202outer iterations on116116instances at the width\-6​n6n,900900\-iteration setting\. Spearmanρ\\rhowith mid\-ranks; intervals from800800bootstrap replicates resampling instances rather than iterations\.The results separate the coverage diagnostics from the model\-quality ones, and the separation is one of size rather than of existence\. Two of the three coverage diagnostics are statistically distinguishable from zero at this sample size,−0\.097\-0\.097forrkr\_\{k\}and−0\.063\-0\.063forNspN\_\{\\mathrm\{sp\}\}, and the sign is the intuitive one: more redundancy, slightly worse odds\. But the effects are five to eight times smaller than those of the model\-quality statistics,−0\.492\-0\.492for the gradient\-scale ratio and\+0\.465\+0\.465for the previous acceptance rate, and a diagnostic that moves the acceptance odds by a hundredth is not a gate\. The directional rank is not distinguishable from zero at all, and for it there is a plain explanation: it is nearly always full, averaging3\.833\.83of44,7\.547\.54of88and14\.4214\.42of1616and attaining the maximum at92%92\\%,90%90\\%and81%81\\%of iterations respectively, so there is little variation for it to explain\. Whatever limits the surrogate here, it is not a training set collapsing onto a subspace\.

The model\-quality diagnostics do work\. The strongest is the ratioκk=‖∇mk‖/‖gkFD‖\\kappa\_\{k\}=\\\|\\nabla m\_\{k\}\\\|/\\\|g\_\{k\}^\{\\mathrm\{FD\}\}\\\|between the surrogate and finite\-difference gradient norms atρ=−0\.492\\rho=\-0\.492, then the acceptance rate at the previous outer iteration at\+0\.465\+0\.465and the training residualϱk\\varrho\_\{k\}at−0\.394\-0\.394\. Each carries the sign a mechanism would predict, an over\-confident or poorly fitted model being the one whose proposals fail, and the effect sizes are usable rather than marginal\. Sorting outer iterations byκk\\kappa\_\{k\}, the acceptance rate falls monotonically across quartiles from0\.4840\.484to0\.3430\.343,0\.1060\.106and0\.0450\.045; split atκk=2\\kappa\_\{k\}=2it is0\.4150\.415below against0\.0760\.076above, a gap of0\.3390\.339with interval\[0\.292,0\.384\]\[0\.292,0\.384\]\. The previous\-iteration acceptance rate separates comparably,0\.4750\.475against0\.1330\.133\. One entry in Table[3](https://arxiv.org/html/2608.24963#S9.T3)looked anomalous and turned out to be mislabelled, so we record what it actually measures\. The decrease ratio atk−1k\-1correlates−0\.155\-0\.155with interval\[−0\.234,−0\.065\]\[\-0\.234,\-0\.065\], the opposite of the sign a trust\-region argument predicts\. The explanation is that the quantity logged is not the trust\-region ratio of the previous*accepted*step but that of the previous*proposal*, and since the surrogate loop exits on its first rejection, that proposal is the rejected one at99\.2%99\.2\\%of iterations; its median value is−0\.202\-0\.202and it exceeds11at eleven of32023202iterations\. Read correctly it is a measure of how badly the last failed proposal overshot, and a worse overshoot atk−1k\-1predicting a worse model atkkis the expected sign, not a reversed one\. Recomputing the statistic over accepted steps only gives\+0\.309\+0\.309, the sign theory predicts, on the10531053of32023202iterations where the previous loop accepted anything at all\. That last figure is the reason we do not gate on it: a trust\-region ratio is undefined at two thirds of the iterations where a gate would have to fire, whereasκk\\kappa\_\{k\}andϱk\\varrho\_\{k\}are defined at all of them\.

Two consequences follow\. The concern raised aboutNspN\_\{\\mathrm\{sp\}\}in Section[4\.3](https://arxiv.org/html/2608.24963#S4.SS3)was understated: it measures no directional coverage and what it does predict is an order of magnitude too weak to act on, and the same holds forrkr\_\{k\}\. And since the framework of Section[8\.4](https://arxiv.org/html/2608.24963#S8.SS4)needs*some*gate,κk\\kappa\_\{k\}andϱk\\varrho\_\{k\}supply one on the strength of an effect five to eight times larger\.

##### Calibrating the gate\.

A predictive signal is not yet a decision rule, so we calibrate one\. Gating the loop off at an outer iteration saves theppevaluations its proposals would have cost and gives up theaadescent steps it would have supplied, each of which the base method must otherwise buy with a finite\-difference gradient costingn\+1n\+1evaluations; the net saving isp−a⁡\(n\+1\)p\-a\(n\+1\)\. Maximising that quantity overκk\\kappa\_\{k\}thresholds on half the instances, chosen at random, gives the same answer at every candidate value: it is negative throughout, and the optimum is not to gate at all\. Gating pays only where the acceptance rate falls below1/\(n\+1\)1/\(n\+1\), that is0\.2000\.200,0\.1110\.111and0\.0590\.059atn=4,8,16n=4,8,16, and even in the least trustworthyκk\\kappa\_\{k\}quartile the observed rates are0\.3550\.355,0\.3820\.382and0\.1140\.114\. Two conventions differ from those of Section[8\.4](https://arxiv.org/html/2608.24963#S8.SS4)and both matter\. These are pooled ratios, total accepted over total proposed, which is what the accounting needs, whereas the quartile series quoted there is a mean of per\-iteration ratios and is lower because90%90\\%of high\-κk\\kappa\_\{k\}iterations make a single proposal; and the quartiles here are taken within each dimension, since the break\-even1/\(n\+1\)1/\(n\+1\)is a per\-dimension quantity\. Every one clears its break\-even, though atn=16n=16by a factor of two rather than by the wide margins of the two smaller dimensions\.

##### Validating it\.

Running the gated solver at four thresholds against ungated assistance on all117117instances confirms the prediction\. The validation set contains the calibration half, which we report rather than conceal: the calibrated answer is a corner solution, not a fitted threshold, and what is being validated is a fresh solver run rather than a re\-scoring of the logs the threshold was chosen on\. Ungated assistance solves109109instances atτ=10−5\\tau=10^\{\-5\}; gating atκk\>1\\kappa\_\{k\}\>1,22,55and3030, which fires at71\.5%71\.5\\%,52\.5%52\.5\\%,41\.6%41\.6\\%and28\.9%28\.9\\%of iterations, solves9797,103103,108108and109109\. The aggressive settings lose outright, the permissive ones converge back to the ungated method as they stop firing, and the paired evaluation difference on jointly solved instances is zero throughout\. Gating atκk\>2\\kappa\_\{k\}\>2avoids11071107proposals but forfeits925925accepted steps, a net loss of about84008400evaluations at the meann\+1≈10\.3n\+1\\approx 10\.3here\. Training the network properly widens the loss rather than narrowing it: the same comparison at the undertrained setting cost six instances atκk\>1\\kappa\_\{k\}\>1where it now costs twelve, because a gate that fires on a model worth trusting forfeits more\.

This is a completed calibration with a negative outcome rather than an open question, and the outcome is informative: safeguarded assistance is already cheap enough that no gate on top of it pays for itself while evaluations are cheap and the model is merely imperfect\. The criterion also says where a gate would pay, namely once acceptance drops below1/\(n\+1\)1/\(n\+1\)\. The inventory model of Section[9\.6](https://arxiv.org/html/2608.24963#S9.SS6)does not supply the counter\-example we expected: withn=8n=8its break\-even is0\.1110\.111and the measured rate is5/31=0\.1615/31=0\.161, above it, so the same rule leaves the gate off there too\. What that model illustrates is the third condition rather than a use for the gate – a surrogate consulted3131times in60,00060\{,\}000replications is not something a gate can improve\. The correlations remain moderate, soκk\\kappa\_\{k\}ranks iterations usefully rather than separating them cleanly; what we can now state is the gate’s operating regime, and that this benchmark lies outside it\.

### 9\.8Is the Network the Problem?

The architecture used throughout is not a free choice, and this section is where it was made\. We varied the two settings controlling capacity and convergence on all thirteen problems atn=8n=8: width in\{n,3​n,6​n,12​n\}\\\{n,3n,6n,12n\\\}and training length in\{300,900\}\\\{300,900\\\}\. An earlier version of this paper reported every surrogate result at width3​n3nand300300iterations; the ablation showed that setting to be undertrained, and the whole benchmark, the loss ablation, the noise sweep and the high\-dimensional runs were repeated at the setting selected here\. Depth is varied separately, in Section[9\.3](https://arxiv.org/html/2608.24963#S9.SS3), and value\-only in both arms, since the analytic parameter gradients of the Sobolev loss in Section[9\.1](https://arxiv.org/html/2608.24963#S9.SS1)are derived for a single hidden layer\.

Training length matters and width barely does\. At300300iterations every width solves66of1313atτ=10−5\\tau=10^\{\-5\}, exactly matching the base method; at900900the same widths solve88,99and99for3​n3n,6​n6nand12​n12n, and the two best configurations, widths6​n6nand12​n12n, are indistinguishable\. Capacity is not the binding constraint – widthnnto12​n12nchanges nothing at300300iterations and little at900900– and convergence is\. We take6​n6nand900900forward\.

That correction is the largest in this paper: the same embedding solves7575instances atτ=10−5\\tau=10^\{\-5\}at width3​n3nand300300iterations against8484at the setting above, so a third of the measured gap between assistance and a model\-based solver was a property of the fit\. What survives is every factor above, at the trained setting: replacement still loses to its base,6565against6767; value\-only training leaves the surrogate below it at6666; the trust\-region method still gains nothing,8686against8888; released interpolation software is ahead at every dimension\. Better training moved the level and left the ordering\. We varied three settings one at a time and did not search the architecture jointly\.

### 9\.9What the Surrogate Costs in Time

Everything above counts objective evaluations, the usual currency in derivative\-free optimisation because in the applications that motivate it one evaluation dominates everything else\. That accounting hides the surrogate’s own cost, which is not negligible for cheap oracles\. We measure it directly\.

One training call, at the dataset caps of Section[9\.1](https://arxiv.org/html/2608.24963#S9.SS1)and with the300300Adam iterations used throughout, takes a median of4343ms atn=4n=4,5757ms atn=8n=8and111111ms atn=16n=16in our pure\-NumPy implementation\. Evaluating the trained surrogate’s gradient costs about18​μ18\\,\\mus and is negligible by comparison; training dominates\. The analytic test objectives, at5\.5​μ5\.5\\,\\mus per call, are some four orders of magnitude cheaper than a training call, so on the benchmark of Section[9\.2](https://arxiv.org/html/2608.24963#S9.SS2)the surrogate variants are far slower in wall\-clock time than the base method even where they are far cheaper in evaluations\. This is expected rather than a defect of the method: it is precisely why the evaluation count is the appropriate currency only for expensive oracles\.

It remains to quantify how expensive an oracle must be\. Over the117117instances the assistance solver performs32113211trainings, issues1065010650proposals and has74887488of them accepted, that is3\.323\.32proposals and2\.332\.33accepted steps per training\. An accepted surrogate step is a descent step obtained without a finite\-difference gradient, so it displaces then\+1n\+1evaluations that gradient would have cost, while every proposal costs one true evaluation whether or not it is accepted\. The net saving per training is therefore about2\.33​\(n\+1\)−3\.322\.33\(n\+1\)\-3\.32evaluations, namely8\.38\.3atn=4n=4,17\.717\.7atn=8n=8and36\.336\.3atn=16n=16, and the surrogate pays for its own time once a single evaluation costs more than16\.716\.7ms atn=4n=4,12\.412\.4ms atn=8n=8and14\.714\.7ms atn=16n=16\. Crediting an accepted step with a whole finite\-difference gradient is generous, so these are lower bounds on the oracle cost at which assistance becomes worthwhile\. Training the network properly raises both sides of that ledger: the yield per training rises by a factor of two and a half and the cost of one training by a factor of three to five, so the threshold barely moves atn=4n=4and rises atn=8n=8andn=16n=16; Section[10](https://arxiv.org/html/2608.24963#S10)shows it continuing to rise beyondn=16n=16\.

Ten milliseconds is a low threshold for the simulation models that motivate this work, since a single replication of the inventory model of Section[9\.6](https://arxiv.org/html/2608.24963#S9.SS6)already exceeds it at the replication counts used there, but a high one for an analytic benchmark function\. The wall\-clock case for a surrogate is therefore governed by the cost of the oracle rather than by the surrogate itself\. An optimised training implementation would lower the threshold roughly in proportion to the speed\-up it achieves\.

### 9\.10Scope and Limitations

One limitation deserves to be stated before the others, because it bounds every claim here\. Our main studies usen≤16n\\leq 16and the probe of Section[9\.2](https://arxiv.org/html/2608.24963#S9.SS2)reachesn=128n=128on thirteen problems with one start each\. Neural surrogates interest the field largely because they are expected to scale where interpolation models do not, and at the dimensions we reach that expectation is not met: a released interpolation solver is ahead at every one, and the margin atn=128n=128is the margin atn=4n=4\. What we cannot claim is that this continues\. The thirteen problems are smooth and unconstrained, one start each at the two largest dimensions, andnnin the thousands, where interpolation’sO⁡\(n2\)O\(n^\{2\}\)model bookkeeping is genuinely prohibitive and subspace or sketching methods take over, is out of reach here: a singlen=128n=128run of the released solver took us up to3333hours\. Whether the ordering we observe persists there is the most important thing this paper leaves open, and our own experience atn=64n=64is a caution against settling it with a reimplementation\. The experiments are limited in the usual ways as well\. The deterministic study uses one run per instance, which is exact because the solvers are deterministic there; the region\-size sweep covers273273cells and the noise sweep200200, each at five seeds; the inventory studies use3030seeds at three starting designs on two model sizes\. The diagnostic study of Section[9\.7](https://arxiv.org/html/2608.24963#S9.SS7)pools outer iterations across problems, so its correlations still mix within\-run and between\-problem variation even though its intervals are computed at the instance level\. A word on the reference methods\. Both are our implementations of the published schemes, so we checked the one that can be checked\. Running released Py\-BOBYQA\[[6](https://arxiv.org/html/2608.24963#bib.bib6)\]on all117117benchmark instances under the same budget and the same counting oracle, the medianlog10\\log\_\{10\}ratio of final objective values between our trust\-region method and the released solver is0\.0000\.000, with ours ahead on3838of117117; neither dominates and the typical difference is nil\. That agreement holds atn≤16n\\leq 16only\. Atn=64n=64it fails badly, our implementation solving55of1313against1212for the released solver, and Section[9\.2](https://arxiv.org/html/2608.24963#S9.SS2)reports what that cost us: a crossover we had measured and believed turned out to be an artefact of our own code\. A reimplementation validated at one problem size is not validated at another, and the dimensions where an interpolation method is under most pressure are exactly where careful engineering separates it from a straightforward version\. The adaptive\-sampling method of Section[9\.6](https://arxiv.org/html/2608.24963#S9.SS6)has no such check, no released ASTRO\-DF being available to us, and its numbers are a reference point rather than a benchmark; released implementations do exist in the simulation\-optimisation testbed literature and running against one is the obvious next step there\. Three things would settle what this paper leaves open\. Ourn=128n=128runs are incomplete, and since the argument for learned surrogates concerns the regime where a full quadratic is impractical,nnin the hundreds is where the question is decided\. A formal benchmark on the CUTEst and Moré–Wild suites, noisy variants included, should use released code throughout; our experience atn=64n=64is why we call that necessary rather than desirable\. And the assistance mechanism itself should be mounted on released software rather than on a base method we wrote: the role distinction predicts that a safeguarded loop attached to Py\-BOBYQA, or to a released adaptive\-sampling solver in the stochastic case, would cost little and gain little, and that prediction is worth testing where a failure would be visible\. All of this quantifies mechanisms rather than state\-of\-the\-art performance, and the wall\-clock figures of Section[9\.9](https://arxiv.org/html/2608.24963#S9.SS9)come from an unoptimised pure\-NumPy trainer and so bound the surrogate’s overhead from above\.

## 10Discussion

Our synthesis yields practical and theoretical implications\.

##### What the experiments support\.

The measurements in Section[9](https://arxiv.org/html/2608.24963#S9)line up with the mechanisms we have argued for\. Holding the surrogate class and the training pipeline fixed, safeguarded assistance improved evaluation efficiency and gradient replacement degraded it\. Uniform error grew quickly with the sampling radius and failed to shrink in proportion at small radii\. Removing the gradient term from the loss cut the acceptance rate from0\.7030\.703to0\.1480\.148\. The stochastic studies add a fourth observation, and it is the one a simulation reader should weigh most heavily\. What noise removes is not the surrogate’s accuracy but the base method’s ability to keep going: across the sweep its runs shorten until they use a fifth of their budget, and the surrogate’s contribution shortens with them: proposals fall by a factor of123123while the rate at which they are accepted falls by a factor of1\.61\.6\. On an actual Monte\-Carlo model every method with a deterministic oracle interface stalls at its starting design; replication repairs the interface and recovers most of the achievable improvement, after which the surrogate adds nothing\. Over9090paired runs the median difference is zero, it stays zero on a larger instance built to leave four times as much room, and the loop is invoked3131times in60,00060\{,\}000replications against1065010650over the benchmark\. The safeguard holds the cost of that idle model to one part in two thousand of the budget\. Assistance does no systematic damage, which is weaker than saying it does none: its median paired difference against the base method is zero, but run by run it is worse on2020of the9090and better on2121, the rest tied\. Replacement, by contrast, trails on8585of9090\. Read with the benchmark, where the same surrogate lifts high\-accuracy solutions from6767to8484of117117over its own base, stays behind the103103a released interpolation solver reaches alone, and lowers our own such solver from8888to8686when attached to it, and where that ordering holds atn=32n=32,n=64n=64andn=128n=128as well, the results bracket the regime in which learned surrogates earn their cost: they close gaps that the base method leaves open, and only those\. On smooth problems the gap they close is worth as much as a competent interpolation model atn≤16n\\leq 16and less than a well\-engineered one above it\. Our experiments therefore establish the mechanism by which a learned local model helps, and locate its usefulness where the base method is weak and the dimension is low enough that a derivative\-free budget can still buy the data the model needs\.

##### What the mechanism faces aboven=100n=100\.

Our experiments stop atn=128n=128, so the regime the motivation is really about is out of reach, but three of the surrogate’s own costs can be measured there without running a solver, and they point the wrong way\. Training one surrogate at the dataset caps used throughout takes140140ms atn=4n=4and3131s atn=128n=128, growing asn2\.0n^\{2\.0\}over the upper half of that range, because the parameter count at width6​n6ngrows as6​n26n^\{2\}\. Each finite\-difference gradient target costsn\+1n\+1evaluations, so the ten targets the training set holds consume a fifth of the evaluation budget at every dimension,12901290of64506450atn=128n=128\. Combining these with the yield measured in Section[9\.9](https://arxiv.org/html/2608.24963#S9.SS9), the oracle cost at which assistance pays for its own wall\-clock time is16\.716\.7,12\.412\.4and14\.714\.7ms atn=4n=4,88and1616, then25\.625\.6,56\.556\.5and104\.4104\.4ms atn=32n=32,6464and128128\. The surrogate’s economics improve with dimension only while the avoided finite\-difference gradient grows faster than the cost of fitting the model, and pastn≈8n\\approx 8they stop doing so\.

The arithmetic behind that turn is worth setting out, because it also explains why the training correction of Section[9\.8](https://arxiv.org/html/2608.24963#S9.SS8)does not carry to high dimension\. The number of training points is fixed by the evaluation budget at10​\(n\+1\)10\(n\+1\), since each gradient target costs a finite\-difference gradient; the parameter count at width6​n6nis6​n2\+12​n\+16n^\{2\}\+12n\+1\. The ratio therefore grows linearly, from55parameters per point atn=8n=8to3939atn=64n=64and7777atn=128n=128\. Widening the network and training it longer cures underfitting atn≤16n\\leq 16, where the ratio is single\-digit, and cannot help where the ratio is already large: the tuned setting is worth1818instances of117117atn≤16n\\leq 16and two of thirteen atn=128n=128\. What limits a learned local model in high dimension on this evidence is not its expressive power but the number of points a derivative\-free budget can buy\.

This is the opposite of the expectation that motivates learned surrogates\. A network’s advantage over interpolation is supposed to be that its parameter count need not track the number of interpolation conditions; ours grows quadratically because we scale width withnn, and a fixed\-width network would trade that against capacity\. The measurements indict our scaling choice, not the idea, which makes an architecture whose cost grows slowly without its local gradient accuracy collapsing the substantive open problem here, rather than the same experiment at largernn\.

The comparison method is not cheap atn=128n=128either: the thirteen released\-solver runs of Section[9\.2](https://arxiv.org/html/2608.24963#S9.SS2)took9090hours between them, from1616seconds on the sphere to3333hours on extended Wood, on an objective whose own evaluation costs microseconds\. Maintaining a257257\-point model for64506450evaluations is itself expensive\. Both figures are properties of Python implementations rather than of the two approaches, so we draw no ranking from them; the point is only that at this dimension the wall\-clock argument does not separate the methods the way the evaluation counts do\.

##### Smoothness vs trend prediction\.

Smoothing reduces oscillations induced by noisy or poorly tuned finite differences; trend prediction produces descent directions that correlate with the objective’s decrease over a neighbourhood\. Sobolev training encodes the latter and, under finite\-difference stencils, penalises surrogate curvature along coordinate directions \(Proposition[3](https://arxiv.org/html/2608.24963#Thmproposition3)\)\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\. With weight decay and dataset capping this stabilises local surrogate descent in a way raw finite differences do not\.

##### Why “replace” can fail even with good approximation\.

Giovannelli et al\. find that surrogate gradients do not reliably improve a strong FD\-BFGS\-based method\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\], which is consistent with the sensitivity of quasi\-Newton updates: inexact gradients affect both the descent direction and the curvature update, and the errors accumulate\. Without uniform error bounds or acceptance safeguards the algorithm over\-trusts the surrogate\.

##### Why “assist” can succeed\.

Assistance designs use the surrogate to propose steps but validate each step by sufficient decrease in the true objective \(Algorithm[1](https://arxiv.org/html/2608.24963#alg1)\)\. This decouples “surrogate imagination” from “true objective acceptance” and controls failure modes\. Taminiau et al\. formalize how repeated successful surrogate steps improve evaluation\-complexity constants via the gain factorη⁡\(S\)\\eta\(S\)\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\]\. Our propositions \(Section[7\.2](https://arxiv.org/html/2608.24963#S7.SS2)\) further clarify that, under bounded gradient error, surrogate steps are more likely to align with true descent, thereby passing the true\-decrease filter\.

##### Generalization as a radius question\.

Both papers implicitly highlight that surrogates are local: Giovannelli et al\. use ball sampling and observe radius effects; Taminiau et al\. cap datasets and operate within stepsize\-controlled neighborhoods\. Our effective generalization radius \(Definition[5](https://arxiv.org/html/2608.24963#Thmdefinition5)\) reframes generalization as “how far can we trust the surrogate?” This enables direct experimental estimation by radius sweeps and provides guidance on region\-size control in learning\-augmented loops\.

##### Speed and cost drawbacks\.

NN surrogates impose training overhead and hyperparameter sensitivity\. Giovannelli et al\. emphasize that NNs can be evaluation\-competitive at high training cost\[[10](https://arxiv.org/html/2608.24963#bib.bib10)\]\. Taminiau et al\. restrict model complexity \(shallow NN\) and use dataset caps and warm starts to limit overhead\[[12](https://arxiv.org/html/2608.24963#bib.bib12)\], but inference about the best surrogate \(NN vs RBF\) remains problem dependent\. Our proposed experiments emphasize measuring both evaluation efficiency and wall\-clock time\.

## 11Conclusion

The two questions of the title have answers, and they are not the same answer\. A neural surrogate improves local approximation because it fits at once every evaluation the run has already paid for, which averages the noise that differencing amplifies, and because gradient information in the training loss restrains the curvature a value\-only fit would invent: remove that term and the acceptance rate falls from0\.7030\.703to0\.1480\.148, taking the instances solved to high accuracy back to exactly the base method’s count\. It improves the*optimisation*only under conditions we delimit empirically rather than prove\. Three factors bound the gain, on evidence of three different kinds\.*Role*we fail on purpose and watch the benefit go: moving the surrogate from proposing candidates to supplying the gradient drops it from8484instances to6565, below the6767of the base method it was meant to help\.*Radius*we establish at the level of approximation rather than of the solver; that the same threshold binds inside a running solver follows from Proposition[2](https://arxiv.org/html/2608.24963#Thmproposition2), not from a separate experiment\. And what an earlier draft called*room*is two things we now keep apart:*headroom*, the performance the base method leaves on the table, which a model\-based solver has little of; and*viability*, whether its acceptance test still functions, which noise destroys\. The eight\-product inventory instance separates them, having four times the headroom and still no gain, so headroom is not sufficient once viability has gone\. Reading Giovannelli et al\. and Taminiau et al\. through the first of them, algorithmic role, accounts for their opposite conclusions without either being wrong: replacement is fragile because surrogate error enters the core updates, whereas safeguarded assistance preserves robustness and can improve the evaluation\-complexity constant through the surrogate\-gain mechanism\.

We proposed a radius\-aware notion of local generalisation with diagnostics for the quality of the local data, and tested it\. Assistance improved evaluation efficiency where replacement degraded it; surrogate reliability fell away quickly outside a bounded region and did not improve without limit inside it; Sobolev training was what made surrogate steps survive validation against the true objective\. Under a stochastic oracle no noise threshold separates the regimes: the largest level at which a problem is still helped ranges over four decades, and what shortens with noise is the base method’s run, not the surrogate’s success rate\.

That last finding is the one we would take forward\. Safeguarding keeps a learned model from damaging a solver, but it offers no protection that the underlying acceptance test does not already provide, and on a real simulation model that limit binds immediately: with a deterministic oracle interface, every variant we tried returned its starting design\. The minimal repair of Section[9\.6](https://arxiv.org/html/2608.24963#S9.SS6)removes the obstacle and restores the ordering observed in the deterministic experiments\. Its ablation corrects our initial guess twice over: what matters is replication, not the shape of the acceptance rule and not the stepsize schedule\. The practical rule we would draw from the three studies together is an ordering of effort\. Repair the sampling and acceptance machinery first, because nothing else matters until the acceptance test can tell a real decrease from noise\. Then ask whether a surrogate pays, and expect the answer to turn on how much headroom the base method leaves and on whether its acceptance test still functions: a great deal of room on the benchmark, essentially none on the inventory model\. Role, radius and room are cheap to check before any network is trained, and where they hold the improvement is real and reproducible; where they do not, no amount of approximation accuracy substitutes for them\. That is why a question about*when*admits a sharper answer than a question about*whether*\.

What we have established is a set of base mechanisms and the conditions under which each operates, measured where every variant can be run on every instance\. Three tasks follow\. The dimensions at which interpolation genuinely breaks down,nnin the thousands, are where Section[10](https://arxiv.org/html/2608.24963#S10)predicts the picture should change; our evidence says the obstacle there is the number of points a derivative\-free budget can buy rather than the network’s capacity, which points at architectures whose parameter count does not tracknn\. Extended Wood shows the set of objectives on which a learned local model beats released interpolation software is not empty, and a diagnostic identifying that set in advance would turn an existence result into a rule\. And the ordering of effort above should be tested against a released adaptive\-sampling solver\[[18](https://arxiv.org/html/2608.24963#bib.bib18)\]rather than the reproduction of its sampling rule we ran, on the CUTEst and Moré–Wild suites\.

## References

- \[1\]A\. R\. Conn, K\. Scheinberg, and L\. N\. Vicente\.*Introduction to Derivative\-Free Optimization*\.SIAM, 2009\.
- \[2\]J\. Larson, M\. Menickelly, and S\. M\. Wild\.Derivative\-free optimization methods\.*Acta Numerica*, 28:287–404, 2019\.
- \[3\]Z\. Zhang\.Derivative\-free optimization \(in Chinese\)\.In Y\.\-X\. Yuan, editor,*China Discipline Development Strategy: Mathematical Optimization*, chapter 7, pages 84–92\. Science Press, Beijing, 2021\. ISBN 978\-7\-03\-065866\-1\.
- \[4\]A\. R\. Conn, K\. Scheinberg, and L\. N\. Vicente\.Global convergence of general derivative\-free trust\-region algorithms to first\- and second\-order critical points\.*SIAM Journal on Optimization*, 20\(1\):387–415, 2009\.
- \[5\]C\. Audet and W\. Hare\.*Derivative\-Free and Blackbox Optimization*\.Springer Series in Operations Research and Financial Engineering\. Springer, Cham, 2017\.
- \[6\]M\. J\. D\. Powell\.The NEWUOA software for unconstrained optimization without derivatives\.In G\. Di Pillo and M\. Roma, editors,*Large\-Scale Nonlinear Optimization*, pages 255–297\. Springer, Boston, 2006\.
- \[7\]P\. Xie and Y\.\-x\. Yuan\.A derivative\-free method using a new underdetermined quadratic interpolation model\.*SIAM Journal on Optimization*, 35\(2\):1110–1133, 2025\. DOI: 10\.1137/23M1582023\.
- \[8\]H\.\-J\. M\. Shi, M\. Q\. Xuan, F\. Oztoprak, and J\. Nocedal\.On the numerical performance of finite\-difference\-based methods for derivative\-free optimization\.*Optimization Methods and Software*, 38\(2\):289–311, 2023\.
- \[9\]A\. S\. Berahas, O\. Sohab, and L\. N\. Vicente\.Full\-low evaluation methods for derivative\-free optimization\.*Optimization Methods and Software*, 38\(2\):386–411, 2023\.
- \[10\]T\. Giovannelli, O\. Sohab, and L\. N\. Vicente\.The limitation of neural nets for approximation and optimization\.*Journal of Global Optimization*, 94:347–376, 2026\. DOI: 10\.1007/s10898\-024\-01426\-9\. \(Also arXiv:2311\.12253\.\)
- \[11\]G\. N\. Grapiglia\.Worst\-case evaluation complexity of a derivative\-free quadratic regularization method\.*Optimization Letters*, 18:195–213, 2024\. DOI: 10\.1007/s11590\-023\-01984\-z\.
- \[12\]T\. Taminiau, G\. N\. Grapiglia, and E\. Massart\.Enhancing finite\-difference\-based derivative\-free optimization with machine learning\.*Optimization Letters*, 20:1159–1177, 2026\. DOI: 10\.1007/s11590\-026\-02281\-1\. \(Also arXiv:2502\.07435\.\)
- \[13\]B\. Irwin, E\. Haber, R\. Gal, and A\. Ziv\.Neural network accelerated implicit filtering: integrating neural network surrogates with provably convergent derivative\-free optimization methods\.In*Proceedings of the 40th International Conference on Machine Learning \(ICML\)*, pages 14376–14389\. PMLR, 2023\.
- \[14\]P\. Xie, Z\. Zhou, and Z\. Zhou\.Objective value change and shape\-based accelerated optimization for the neural network approximation\.arXiv:2508\.20290, 2025\.
- \[15\]Q\. Zhang and P\. Xie\.On the relationship betweenΛ\\Lambda\-poisedness in derivative\-free optimization and outliers in local outlier factor\.arXiv:2407\.17529, 2024\.
- \[16\]W\. M\. Czarnecki, S\. Osindero, M\. Jaderberg, G\. Swirszcz, and R\. Pascanu\.Sobolev training for neural networks\.In*Advances in Neural Information Processing Systems \(NeurIPS\)*, 2017\.\(Also available as arXiv:1706\.04859\.\)
- \[17\]N\. I\. M\. Gould, D\. Orban, and Ph\. L\. Toint\.CUTEst: a constrained and unconstrained testing environment with safe threads for mathematical optimization\.*Computational Optimization and Applications*, 60\(3\):545–557, 2015\.
- \[18\]S\. Shashaani, F\. S\. Hashemi, and R\. Pasupathy\.ASTRO\-DF: a class of adaptive sampling trust\-region algorithms for derivative\-free stochastic optimization\.*SIAM Journal on Optimization*, 28\(4\):3145–3176, 2018\.
- \[19\]J\. J\. Moré, B\. S\. Garbow, and K\. E\. Hillstrom\.Testing unconstrained optimization software\.*ACM Transactions on Mathematical Software*, 7\(1\):17–41, 1981\.
- \[20\]J\. J\. Moré and S\. M\. Wild\.Benchmarking derivative\-free optimization algorithms\.*SIAM Journal on Optimization*, 20\(1\):172–191, 2009\.
- \[21\]E\. D\. Dolan and J\. J\. Moré\.Benchmarking optimization software with performance profiles\.*Mathematical Programming*, 91\(2\):201–213, 2002\.