# Surrogate with derivatives **URL:** \ **Category:** Modelling & Simulations\ **Tags:** question, surrogates\ **Created:** [May 25, 2026, 1:59pm UTC](https://discourse.julialang.org/t/surrogate-with-derivatives/137277 "2026-05-25T13:59:17Z")\ **Posts on this page:** 7\ **Page:** 1

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\ **Post date:** [May 25, 2026, 1:59pm UTC](https://discourse.julialang.org/t/surrogate-with-derivatives/137277/1 "2026-05-25T13:59:17Z")

I am solving an implicit function f(x, y(x)) = 0 with x, y \in \mathbb{R}^n (typical n is 50–200), using iterative methods (trust region), at many x’s. This **question is about finding a good initial value** y\_0(x) for each x, good initial values speed up the solution. Suppose I have solved the problem for x\_i, i = 1, \dots, k, obtaining the corresponding y\_i, and I also have computed all the information I need for Jacobian-vector products vJ\_i and J\_i v with J\_i = \partial y/\partial x. Importantly, typically k \< n: I have fewer points than dimensions. I can explore a lot of nice models with the excellent Surrogates.jl, but AFAICT only [GEK(PLS)](https://docs.sciml.ai/Surrogates/stable/gekpls/) support derivatives, and my understanding is that I would have to go elementwise (recall: y(x): \mathbb{R}^n \to \mathbb{R}^n). Otherwise, all the models I am familiar with use just values, not derivatives. It feels like a waste not to use them since they are already computed. I am kind of fishing suggestions here that I could start exploring. Again, the surrogate does not have to be very good _per se_, it is just a starting point for an iterative method. A naive way I thought of is obtaining first-order approximations from _all_ points, y^{(i)}(x) = y\_i + \partial y/\partial x \cdot (x - x\_i) and averaging them out somehow. But this looks very _ad-hoc_, and putting discipline on this would be nice, surely someone has invented a good method for this already. ---

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\ **Post date:** [May 25, 2026, 2:55pm UTC](https://discourse.julialang.org/t/surrogate-with-derivatives/137277/2 "2026-05-25T14:55:00Z")

> [@Tamas\_Papp](#): > > Otherwise, all the models I am familiar with use just values, not derivatives. It feels like a waste not to use them since they are already computed. Another choice is the Gradient-Enhanced Neural Network (GENN), GEKPLS is generally the preferred method since pure NNs are rarely a good surrogate. > [@Tamas\_Papp](#): > > I can explore a lot of nice models with the excellent [Surrogates.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/Surrogates), but AFAICT only [GEK(PLS)](https://docs.sciml.ai/Surrogates/stable/gekpls/) support derivatives, and my understanding is that I would have to go elementwise (recall: y(x): \mathbb{R}^n \to \mathbb{R}^n). It should do vector output? If not, just open an issue. IIRC there isn’t anything in the math that doesn’t apply to a vector output. ---

> [@Tamas\_Papp](#): > > A naive way I thought of is obtaining first-order approximations from _all_ points, > > y^{(i)}(x) = y\_i + \partial y/\partial x \cdot (x - x\_i) > > and averaging them out somehow. But this looks very _ad-hoc_, and putting discipline on this would be nice, surely someone has invented a good method for this already. Or a Hermite interpolation in higher dimensions. ---

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\ **Post date:** [May 25, 2026, 7:04pm UTC](https://discourse.julialang.org/t/surrogate-with-derivatives/137277/4 "2026-05-25T19:04:13Z")

> [@Tamas\_Papp](#): > > with x, y \in \mathbb{R}^n (typical n is 50–200) […] Importantly, typically k \< n: I have fewer points than dimensions. Dimensions that high seem to rule out any kind of tensor-product method like Hermite polynomials or splines, and make interpolation difficult in general. For low-accuracy interpolation in high dimensions, neural networks such as the GENN methods Chris mentioned are one option. Another option is radial basis functions (RBFs), and there are “Hermite RBF” variants that exploit derivative information, e.g. [Vaillant (2013)](https://mobile.rodolphe-vaillant.fr/images/pdfs/hrbf_summary.pdf) or [Fashamiha and Salac (2025)](https://arxiv.org/abs/2503.05752v2). (Although RBFs are often written down for scalar-valued functions, the same formulations are trivially extensible to vector-valued functions simply by computing separate RBF coefficients for each component of the function, i.e. vector-valued coefficients. The cost scales linearly with the number of components, and it is [embarrassingly parallel](https://en.wikipedia.org/wiki/Embarrassingly_parallel).) There seems to be an ancient (Julia 0.5) implementation of the Vaillant HRBF algorithm [SpatialFields.jl](https://github.com/rdeits/SpatialFields.jl). However, RBFs are pretty straightforward to implement yourself, especially with an AD package like ForwardDiff.jl to compute the RBF derivatives for you in the HRBF case. You’ll get a 2(n+1)k \times 2(n+1)k system of equations, or rather n of them (one for each component of y), which in your case should be small enough to solve directly in a minute or less (each system being a few thousand by a few thousand, solved in \lesssim 1 s). Not sure if that is fast enough for you? ---

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\ **Post date:** [May 26, 2026, 8:41am UTC](https://discourse.julialang.org/t/surrogate-with-derivatives/137277/5 "2026-05-26T08:41:35Z")

> [@stevengj](#): > > Another option is radial basis functions (RBFs), and there are “Hermite RBF” variants that exploit derivative information, e.g. [Vaillant (2013)](https://mobile.rodolphe-vaillant.fr/images/pdfs/hrbf_summary.pdf) or [Fashamiha and Salac (2025)](https://arxiv.org/abs/2503.05752v2). GEKPLS is Gradient-Enhanced Kriging, i.e. a Hermite RBF. ---

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\ **Post date:** [May 26, 2026, 11:38am UTC](https://discourse.julialang.org/t/surrogate-with-derivatives/137277/6 "2026-05-26T11:38:47Z")

Thanks for all the answers! After experimenting a bit with my problem I think that my main problem will be with _extra_polation, not _inter_polation. That is, often the point at which I want a surrogate will be outside (the convex hull of) the existing points. What seems to work decently in practice (for the simplified problem I am using for testing) is \hat{y}(x) = (1-\alpha) \cdot \bar{y} + \alpha \cdot (y\_c + \frac{\partial}{\partial x}y(x\_c) \cdot (x - x\_c)) where \bar{x}, \bar{y} are the means of existing points, x\_c is the closest point to x, and \alpha \in [0,1] is a function of x that decreases in \| x - \bar{x} \| / \| x\_c - \bar{x} \|. That is, use a linear approximation from the closest point, but to avoid giving totally crazy estimates just because of a wacky local derivative, dampen that to the mean. ---

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\ **Post date:** [May 26, 2026, 11:42am UTC](https://discourse.julialang.org/t/surrogate-with-derivatives/137277/7 "2026-05-26T11:42:02Z")

> [@ChrisRackauckas](#): > > GEKPLS is Gradient-Enhanced Kriging, i.e. a Hermite RBF. Good point, although RBFs (unlike Kriging) are traditionally augmented with some low-order polynomials, and often use a wider variety of basis functions than Gaussians.