Non-Linear Surrogate Model Extension

• Non-Linear Surrogate Model Extension is a framework that augments conventional surrogate models with nonlinear representations to capture curved dynamics and regime changes.
• It leverages techniques such as quadratic Bézier trajectories, operator lifting, and local spectral expansions to improve surrogate fidelity in diverse applications.
• NL-SME is applied in federated learning gradient inversion, nonlinear model predictive control warm-starting, and system identification, offering enhanced accuracy and efficiency.

Non-Linear Surrogate Model Extension (NL-SME) denotes a class of surrogate-modeling strategies in which a baseline surrogate representation is extended so that nonlinear trajectories, nonlinear dynamics, local regimes, or nonlinear operator structure are represented more faithfully than with a single global linear or weakly nonlinear approximation. In the supplied literature, the term is explicit in the federated-learning attack model of (Xia et al., 26 Sep 2025), where a linear surrogate trajectory is replaced by a learnable quadratic Bézier curve. A broader reading of related work suggests the same extension principle in surrogate dynamics warm-starting for nonlinear model predictive control (Beja-Battais et al., 7 Apr 2026), quadratic-bilinear identification from input-output data (Karachalios et al., 2021), on-the-fly reconstruction of nonlinear dynamics from local linearizations (Olucha et al., 31 Mar 2025), local regime-aware aerodynamic surrogates (Dupuis et al., 2019), recursive spectral embeddings (Marelli et al., 2020), state-space Kriging for stochastic dynamics (Chenga et al., 2024), and conditional Gaussian surrogate or preconditioner constructions for complex nonlinear systems (Chen et al., 2021).

1. Defining pattern and scope

The common structural feature is that the surrogate is not merely fitted as a static approximation of outputs; instead, it is augmented so that nonlinear structure is introduced where a simpler surrogate is known to fail. In (Xia et al., 26 Sep 2025), this means replacing straight-line interpolation in parameter space with a quadratic Bézier path. In (Beja-Battais et al., 7 Apr 2026), it means learning surrogate dynamics only for an auxiliary optimization layer, while preserving the original full-scale NMPC solve. In (Karachalios et al., 2021), it means extending linear Loewner interpolation by inferring quadratic and bilinear operators after lifting. In (Olucha et al., 31 Mar 2025), it means reconstructing a nonlinear field from interpolated Jacobians rather than regressing the field directly. In (Dupuis et al., 2019) and (Marelli et al., 2020), it means replacing a single global surrogate by local or recursive pieces adapted to regime changes, discontinuities, or localized complexity. In (Chenga et al., 2024) and (Chen et al., 2021), it means moving from static emulation to dynamical or conditional-statistics formulations.

This suggests that NL-SME is best understood as an extension principle rather than a single algorithm. The extension may occur in trajectory space, operator space, state space, regime space, or conditional-statistics space. Across these settings, the extension is introduced because a baseline surrogate is too restrictive: the linear trajectory assumption in multi-step gradient inversion, global POD/GPR under transonic regime changes, global spectral expansions under discontinuities, or direct brute-force nonlinear simulation under partial observations and stochastic forcing.

2. Main formalizations of nonlinear extension

The supplied works instantiate the extension principle through several mathematically distinct constructions.

Formulation	Extension mechanism	Reported role
(Xia et al., 26 Sep 2025)	Quadratic Bézier surrogate trajectory with control point $P_1$	Multi-step gradient inversion in Federated Learning
(Beja-Battais et al., 7 Apr 2026)	Learned surrogate dynamics $m_\theta$ inside an auxiliary warm-start NLP	Acceleration of full-scale NMPC
(Karachalios et al., 2021)	Loewner fit of linear part plus least-squares recovery of $H$ and $N$	Quadratic-bilinear surrogate identification
(Olucha et al., 31 Mar 2025)	RBF interpolation of local linearizations plus line integral reconstruction	On-the-fly surrogation of nonlinear dynamics
(Dupuis et al., 2019)	Regime classifier plus local POD/GPR surrogates	Multi-regime aerodynamic prediction
(Marelli et al., 2020)	Recursive partitioning with local spectral residual expansions	Piecewise continuous surrogate for localized or non-stationary responses
(Chenga et al., 2024)	Componentwise Kriging of $\dot X=f(X,U)$ in state space	Emulation of stochastic nonlinear dynamics
(Chen et al., 2021)	Conditional Gaussian split $u=(X,Y)$ with closed conditional statistics	Cheap surrogate model and fast preconditioner

At the level of equations, the nonlinear extension takes markedly different forms. The explicit NL-SME attack paper uses

$$\hat{w}(t,P_1)=(1-t)^2w_0+2(1-t)tP_1+t^2w_T,\qquad t\in[0,1],$$

which strictly contains linear interpolation as the special case $P_1=(w_0+w_T)/2$ (Xia et al., 26 Sep 2025). The quadratic-bilinear identification framework uses

$$\dot{x}(t)=Ax(t)+H\bigl(x(t)\otimes x(t)\bigr)+Nx(t)u(t)+Bu(t),\qquad y(t)=Cx(t),$$

after lifting analytic nonlinearities into QB form (Karachalios et al., 2021). The on-the-fly surrogate of (Olucha et al., 31 Mar 2025) reconstructs nonlinear dynamics by

$\hat f(\hat x,u,\eta)=\hat F(\hat x,u,\eta)\begin{bmatrix}\hat x\u\end{bmatrix},\qquad \hat F(\hat x,u,\eta)=\int_0^1 \mathcal I^D(\lambda \hat x,\lambda u,\eta)\,d\lambda,$

so the surrogate is computed only at evaluation time. The S2K emulator writes the nonlinear stochastic system as $m_\theta$0 and learns each component of $m_\theta$1 by Kriging in state space (Chenga et al., 2024). The CGNS construction splits the state as $m_\theta$2 and imposes that, conditional on the history of $m_\theta$3, the hidden sector $m_\theta$4 is linear-Gaussian, which yields closed filter and smoother equations (Chen et al., 2021).

A plausible implication is that NL-SME is not tied to one surrogate family. Neural surrogates, RBF interpolants, Kriging, Loewner realizations, local POD/GPR blocks, spectral expansions, and conditional Gaussian closures all appear in the supplied literature.

3. Explicit NL-SME in multi-step gradient inversion

The most direct use of the term appears in (Xia et al., 26 Sep 2025). The paper addresses multi-step Gradient Inversion Attacks in FedAVG, where the attacker observes an initial model $m_\theta$5, a final model $m_\theta$6, and the local sample size $m_\theta$7, and seeks to reconstruct private client data from the aggregated update. The baseline SME models a surrogate parameter point by linear interpolation,

$m_\theta$8

and matches the gradient direction to $m_\theta$9. The paper argues that this is fundamentally limited because SGD trajectories on non-convex objectives are generally curved rather than linear, and it states a gradient-direction bias lower bound under Lipschitz assumptions.

NL-SME replaces the line segment by the quadratic Bézier path above and optimizes dummy data $H$0, the curve parameter $H$1, the control point $H$2, and a coordinate-wise scaling vector $H$3. Its objective is

$H$4

with total variation regularization, midpoint regularization

$H$5

and dvec scaling

$H$6

The initialization is conservative: $H$7 and the paper specifies hierarchical learning rates satisfying $H$8.

On FEMNIST with $H$9, the reported similarity loss is $N$0 for NL-SME versus $N$1 for SME, $N$2 for IG, and $N$3 for FEDLEAK, with PSNR/SSIM of $N$4. On CIFAR-100 with the same setting, NL-SME reaches $N$5, PSNR $N$6, and SSIM $N$7, compared with SME’s $N$8, $N$9, and $\dot X=f(X,U)$0. Runtime and memory remain close to SME: for CIFAR-100 at $\dot X=f(X,U)$1, NL-SME uses $\dot X=f(X,U)$2 GB and $\dot X=f(X,U)$3 min, while SME uses $\dot X=f(X,U)$4 GB and $\dot X=f(X,U)$5 min. The paper therefore treats nonlinear trajectory modeling as a privacy-relevant strengthening of surrogate attacks rather than as a mere architectural variation (Xia et al., 26 Sep 2025).

4. Surrogate dynamics extension in nonlinear control

In (Beja-Battais et al., 7 Apr 2026), the nonlinear extension is formulated as Surrogate Dynamics Optimization for full-scale NMPC. The stated motivation is that the practical bottleneck is often repeated evaluation of the embedded process model $\dot X=f(X,U)$6, especially for high-dimensional systems, multi-time-scale or stiff dynamics, DAE-based models, long-horizon control problems, and industrial settings where in-house solvers are used and advanced automatic differentiation or sparsity tools may not be available. Standard warm starts such as shift initialization are reported to be insufficient when the exogenous input trajectory $\dot X=f(X,U)$7, the cost function $\dot X=f(X,U)$8 or $\dot X=f(X,U)$9, or the constraints $u=(X,Y)$0 change abruptly.

The workflow has three stages. First, a surrogate dynamics model $u=(X,Y)$1 is trained offline on trajectories generated by the full-scale DAE solver $u=(X,Y)$2, using recursive multi-step loss over a horizon $u=(X,Y)$3: $u=(X,Y)$4 Second, at inference time the surrogate replaces $u=(X,Y)$5 in a lightweight single-shooting NLP, solved with a first-order method such as projected gradient descent. The surrogate stage is deliberately limited to a small part of the total computational budget, typically $u=(X,Y)$6, and the paper also reports significant acceleration even with $u=(X,Y)$7 budget. Third, the surrogate solution $u=(X,Y)$8 is used to warm-start the original full NMPC problem

$u=(X,Y)$9

so the final solution still comes from the trusted full-scale formulation.

The case study is 24-hour optimal load-following control of a pressurized water reactor. The paper evaluates MLP, RNN, and Transformer-based encoder-decoder surrogates and reports that the Transformer performs best because it better captures long-horizon temporal dependencies. The training setup uses $$\hat{w}(t,P_1)=(1-t)^2w_0+2(1-t)tP_1+t^2w_T,\qquad t\in[0,1],$$0 transients, each of 24 hours, with an $$\hat{w}(t,P_1)=(1-t)^2w_0+2(1-t)tP_1+t^2w_T,\qquad t\in[0,1],$$1 training and $$\hat{w}(t,P_1)=(1-t)^2w_0+2(1-t)tP_1+t^2w_T,\qquad t\in[0,1],$$2 validation split, hyperparameter optimization with Optuna over 100 trials, a maximum model size of $$\hat{w}(t,P_1)=(1-t)^2w_0+2(1-t)tP_1+t^2w_T,\qquad t\in[0,1],$$3 parameters, and training in PyTorch on an NVIDIA Tesla V100-SXM2-32GB. Relative to behavior cloning, surrogate training data generation took less than 1 hour, whereas generating BC training data required solving 300 full-scale OCPs, producing about 1 million trajectories and about 100 hours of cumulative compute. Under a main comparison budget of 20 full-scale iterations, about 20 minutes, SDO consistently outperforms shift initialization and BC in convergence speed and objective value after a fixed number of gradient steps, while preserving safety because the full solver remains the final decision stage (Beja-Battais et al., 7 Apr 2026).

5. Operator lifting and on-the-fly nonlinear reconstruction

The framework of (Karachalios et al., 2021) constructs quadratic-bilinear surrogates non-intrusively from input-output measurements. Its central claim is that many nonlinear systems with analytic nonlinearities can be rewritten exactly as QB systems by lifting, rather than approximately by Taylor truncation or Carleman linearization. The linear part is first identified from first-kernel data using the classical Loewner framework; then the supplementing quadratic and bilinear operators are recovered from second-kernel measurements by least squares, with

$$\hat{w}(t,P_1)=(1-t)^2w_0+2(1-t)tP_1+t^2w_T,\qquad t\in[0,1],$$4

and the requirement $$\hat{w}(t,P_1)=(1-t)^2w_0+2(1-t)tP_1+t^2w_T,\qquad t\in[0,1],$$5. Because the system may be ill-conditioned or rank-deficient, the paper recommends truncated SVD or a Moore–Penrose pseudoinverse in regularized form. In the nonlinear RC ladder benchmark, the reduced linear order is chosen as $$\hat{w}(t,P_1)=(1-t)^2w_0+2(1-t)tP_1+t^2w_T,\qquad t\in[0,1],$$6, the second-kernel approximation error is reported as around $$\hat{w}(t,P_1)=(1-t)^2w_0+2(1-t)tP_1+t^2w_T,\qquad t\in[0,1],$$7 in the infinity norm, and the time-domain response error is about $$\hat{w}(t,P_1)=(1-t)^2w_0+2(1-t)tP_1+t^2w_T,\qquad t\in[0,1],$$8.

The on-the-fly surrogation method of (Olucha et al., 31 Mar 2025) starts from local linearization snapshots rather than trajectories or output kernels. For the nonlinear system $$\hat{w}(t,P_1)=(1-t)^2w_0+2(1-t)tP_1+t^2w_T,\qquad t\in[0,1],$$9, it stores distinct operating points $P_1=(w_0+w_T)/2$0 and local Jacobian-based linearizations $P_1=(w_0+w_T)/2$1. An RBF interpolant $P_1=(w_0+w_T)/2$2 is built so that $P_1=(w_0+w_T)/2$3, optionally with a polynomial tail, and the nonlinear dynamics are reconstructed by integrating the interpolated Jacobian along the line segment from the origin to the query point. The model is “on-the-fly” because the global surrogate is never written out analytically; it is evaluated numerically only when needed. In the controlled Van der Pol example, exact reconstruction from perfect local linearizations reproduces the true nonlinear trajectories, and a surrogate built from only $P_1=(w_0+w_T)/2$4 snapshots still closely matches them. In the parameter-varying test with $P_1=(w_0+w_T)/2$5, $P_1=(w_0+w_T)/2$6 observations suffice for interpolation over parameter space. In the 10-state mass-spring-damper example, the method uses $P_1=(w_0+w_T)/2$7 random linearizations in an 11-dimensional space and reports mean simulation time about 4 seconds, an average of 390 online integral evaluations per simulation, and RMSE mostly below $P_1=(w_0+w_T)/2$8, with some outliers up to $P_1=(w_0+w_T)/2$9 (Olucha et al., 31 Mar 2025).

These two works illustrate two distinct NL-SME mechanisms: exact or lifted operator enrichment, and nonlinear field reconstruction from local differential information.

6. Local, hierarchical, and stochastic extensions

Several supplied papers extend surrogate modeling by introducing locality, recursive structure, or conditional-statistics structure. In aerodynamics, the Local Decomposition Method of (Dupuis et al., 2019) addresses discontinuous signals, moving shocks, bifurcations, and regime changes by combining a physical shock sensor, Gaussian Mixture Model clustering, Gaussian Process Classification of input space, entropy-based active resampling, and local POD/GPR surrogates. For the RAE2822 transonic airfoil, the reported global results improve from $$\dot{x}(t)=Ax(t)+H\bigl(x(t)\otimes x(t)\bigr)+Nx(t)u(t)+Bu(t),\qquad y(t)=Cx(t),$$0 and $$\dot{x}(t)=Ax(t)+H\bigl(x(t)\otimes x(t)\bigr)+Nx(t)u(t)+Bu(t),\qquad y(t)=Cx(t),$$1 to $$\dot{x}(t)=Ax(t)+H\bigl(x(t)\otimes x(t)\bigr)+Nx(t)u(t)+Bu(t),\qquad y(t)=Cx(t),$$2 and $$\dot{x}(t)=Ax(t)+H\bigl(x(t)\otimes x(t)\bigr)+Nx(t)u(t)+Bu(t),\qquad y(t)=Cx(t),$$3 for $$\dot{x}(t)=Ax(t)+H\bigl(x(t)\otimes x(t)\bigr)+Nx(t)u(t)+Bu(t),\qquad y(t)=Cx(t),$$4, and from $$\dot{x}(t)=Ax(t)+H\bigl(x(t)\otimes x(t)\bigr)+Nx(t)u(t)+Bu(t),\qquad y(t)=Cx(t),$$5 and $$\dot{x}(t)=Ax(t)+H\bigl(x(t)\otimes x(t)\bigr)+Nx(t)u(t)+Bu(t),\qquad y(t)=Cx(t),$$6 to $$\dot{x}(t)=Ax(t)+H\bigl(x(t)\otimes x(t)\bigr)+Nx(t)u(t)+Bu(t),\qquad y(t)=Cx(t),$$7 and $$\dot{x}(t)=Ax(t)+H\bigl(x(t)\otimes x(t)\bigr)+Nx(t)u(t)+Bu(t),\qquad y(t)=Cx(t),$$8 for $$\dot{x}(t)=Ax(t)+H\bigl(x(t)\otimes x(t)\bigr)+Nx(t)u(t)+Bu(t),\qquad y(t)=Cx(t),$$9. In the transonic regime, $\hat f(\hat x,u,\eta)=\hat F(\hat x,u,\eta)\begin{bmatrix}\hat x\u\end{bmatrix},\qquad \hat F(\hat x,u,\eta)=\int_0^1 \mathcal I^D(\lambda \hat x,\lambda u,\eta)\,d\lambda,$0 for $\hat f(\hat x,u,\eta)=\hat F(\hat x,u,\eta)\begin{bmatrix}\hat x\u\end{bmatrix},\qquad \hat F(\hat x,u,\eta)=\int_0^1 \mathcal I^D(\lambda \hat x,\lambda u,\eta)\,d\lambda,$1 improves from $\hat f(\hat x,u,\eta)=\hat F(\hat x,u,\eta)\begin{bmatrix}\hat x\u\end{bmatrix},\qquad \hat F(\hat x,u,\eta)=\int_0^1 \mathcal I^D(\lambda \hat x,\lambda u,\eta)\,d\lambda,$2 to $\hat f(\hat x,u,\eta)=\hat F(\hat x,u,\eta)\begin{bmatrix}\hat x\u\end{bmatrix},\qquad \hat F(\hat x,u,\eta)=\int_0^1 \mathcal I^D(\lambda \hat x,\lambda u,\eta)\,d\lambda,$3.

The Stochastic Spectral Embedding method of (Marelli et al., 2020) targets responses that are too nonlinear, localized, or non-stationary for a single global spectral surrogate. It recursively partitions the domain into subdomains $\hat f(\hat x,u,\eta)=\hat F(\hat x,u,\eta)\begin{bmatrix}\hat x\u\end{bmatrix},\qquad \hat F(\hat x,u,\eta)=\int_0^1 \mathcal I^D(\lambda \hat x,\lambda u,\eta)\,d\lambda,$4, fits local spectral residual expansions $\hat f(\hat x,u,\eta)=\hat F(\hat x,u,\eta)\begin{bmatrix}\hat x\u\end{bmatrix},\qquad \hat F(\hat x,u,\eta)=\int_0^1 \mathcal I^D(\lambda \hat x,\lambda u,\eta)\,d\lambda,$5, and assembles a piecewise continuous approximation

$\hat f(\hat x,u,\eta)=\hat F(\hat x,u,\eta)\begin{bmatrix}\hat x\u\end{bmatrix},\qquad \hat F(\hat x,u,\eta)=\int_0^1 \mathcal I^D(\lambda \hat x,\lambda u,\eta)\,d\lambda,$6

The paper reports that, at $\hat f(\hat x,u,\eta)=\hat F(\hat x,u,\eta)\begin{bmatrix}\hat x\u\end{bmatrix},\qquad \hat F(\hat x,u,\eta)=\int_0^1 \mathcal I^D(\lambda \hat x,\lambda u,\eta)\,d\lambda,$7, SSE achieves at least an order-of-magnitude smaller error than PCE on a 1D localized nonlinear function, and that it significantly outperforms PCE on a 100D analytical function and on a 3D snap-through truss with discontinuity.

The state space Kriging model of (Chenga et al., 2024) learns the instantaneous derivative map $\hat f(\hat x,u,\eta)=\hat F(\hat x,u,\eta)\begin{bmatrix}\hat x\u\end{bmatrix},\qquad \hat F(\hat x,u,\eta)=\int_0^1 \mathcal I^D(\lambda \hat x,\lambda u,\eta)\,d\lambda,$8 rather than a long lagged input-output mapping. It uses componentwise Kriging with a Matérn-5/2 kernel, sparse sample selection by the acquisition criterion

$\hat f(\hat x,u,\eta)=\hat F(\hat x,u,\eta)\begin{bmatrix}\hat x\u\end{bmatrix},\qquad \hat F(\hat x,u,\eta)=\int_0^1 \mathcal I^D(\lambda \hat x,\lambda u,\eta)\,d\lambda,$9

and tailored training trajectories generated by magnifying excitation variability. Reported errors include $m_\theta$00 to $m_\theta$01 for a nonlinear quarter-car benchmark, $m_\theta$02 and $m_\theta$03 for displacement and velocity in a Duffing oscillator, and $m_\theta$04, $m_\theta$05, and $m_\theta$06 for displacement, velocity, and hysteretic variable in a Bouc–Wen oscillator.

The Conditional Gaussian Nonlinear System of (Chen et al., 2021) provides a different nonlinear extension: it preserves nonlinear self-interactions among resolved variables $m_\theta$07, keeps $m_\theta$08 cross-couplings, and approximates only the self-interactions among unresolved variables so that, conditional on the history of $m_\theta$09, the law of $m_\theta$10 is Gaussian with explicit filter and smoother formulas. The paper presents CGNS simultaneously as a cheap surrogate model and a fast preconditioner for data assimilation, parameter estimation with uncertainty quantification, PDF reconstruction, trajectory sampling, and linear response prediction.

7. Assumptions, limitations, and recurring misconceptions

Across the supplied literature, NL-SME does not imply unconditional replacement of the original model. In (Beja-Battais et al., 7 Apr 2026), the surrogate is explicitly a warm-start layer and “not a replacement for the true solve”; the final solution remains the original NMPC solution. In (Chen et al., 2021), the surrogate is valuable precisely because it preserves much of the nonlinear physics while making part of the system analytically tractable. In (Olucha et al., 31 Mar 2025), the method is not a reduced-order model and does not simply stitch local linears; it reconstructs a nonlinear field from interpolated Jacobians. In (Xia et al., 26 Sep 2025), NL-SME is not a controller or predictor but an attack model for reconstructing private data from multi-step updates.

Several limitations recur. Surrogate quality depends on data coverage and representational adequacy. If surrogate accuracy is too low, SDO may fail to accelerate or may degrade the warm start (Beja-Battais et al., 7 Apr 2026). The least-squares system for QB identification may be ill-conditioned, rank-deficient, and sensitive to noise (Karachalios et al., 2021). The on-the-fly linearization method requires online numerical integration, tuning of the width parameter $m_\theta$11, and does not yet provide error bounds (Olucha et al., 31 Mar 2025). S2K still faces scalability issues in high-dimensional state spaces and ignores cross-output correlations by modeling components independently (Chenga et al., 2024). CGNS depends on the quality of the chosen closure; if the closure is poor, parameter estimates or response predictions can be biased (Chen et al., 2021). In the federated-learning setting, NL-SME remains vulnerable to defenses that alter update directions and magnitudes; the paper explicitly notes that adaptive optimizers such as AdamW can disrupt gradient-direction-based attacks (Xia et al., 26 Sep 2025).

A recurring misconception is that “nonlinear extension” necessarily means a monolithic black-box regressor. The supplied literature shows otherwise. The nonlinearity may arise from a Bézier trajectory, lifted quadratic-bilinear operators, a line integral of interpolated Jacobians, hard regime partitioning with local reduced models, recursive local spectral expansions, state-space GP closure, or conditional Gaussian decomposition. Another misconception is that greater expressiveness automatically yields formal guarantees. Several papers emphasize the opposite: strong convexity assumptions in the SDO bound are unlikely to hold in practice, noise sensitivity remains open for LoewnerQB, and the on-the-fly method offers convergence intuition rather than formal error estimates (Beja-Battais et al., 7 Apr 2026, Karachalios et al., 2021, Olucha et al., 31 Mar 2025).

Taken together, these works show that NL-SME is most coherently understood as a family of structurally nonlinear surrogate extensions introduced when global linearization, global interpolation, or direct policy imitation is too weak for the target dynamics, optimization path, or regime structure.