Multifidelity Trajectory Simulation Estimators

• Multifidelity trajectory simulation estimators are methods that combine a high-fidelity simulator with lower-cost models to efficiently estimate trajectory-dependent quantities under budget constraints.
• They leverage techniques such as regression, Bayesian conditioning, control variates, and multilevel telescoping to minimize variance and mean-squared error.
• Applied to problems in parametric ODEs, uncertainty propagation, and risk quantification, these estimators deliver significant computational speed-ups while ensuring accuracy.

Multifidelity trajectory simulation estimators are statistical constructions that combine one trusted high-fidelity simulator with cheaper lower-fidelity models to estimate trajectory-dependent quantities of interest under explicit computational budgets. In the literature, the target may be a scalar expectation such as $\mu=\mathbb{E}[f(Y)]$, the full probabilistic response $p(y_H)$, trajectory statistics such as means, standard deviations, and peak magnitudes, stacked multi-output means, variances, and Sobol main effect variances, or path-integral control and safety probabilities (Xu et al., 2021, Keshavarzzadeh et al., 2018, Nitzler et al., 2020, Dixon et al., 2023, Wang et al., 9 Oct 2025, Warner et al., 19 Sep 2025). The common premise is that low-fidelity trajectory information, if exploited through regression, Bayesian conditioning, control variates, multilevel telescoping, or hyperparameter tuning, can reduce estimator variance or mean-squared error for a fixed cost.

1. Problem classes and mathematical setting

A recurring setup uses one high-fidelity simulator $Y$ with cost $c_0$ and $n$ lower-fidelity surrogates $X^{(i)}$, $i=1,\dots,n$, with costs $c_i$, all sampled subject to a total budget

$\sum_{i=0}^n c_i\,n_i \le B.$

The objective is to estimate $\mu:=\mathbb{E}[f(Y)]$ under that cost constraint (Xu et al., 2021). In the trajectory-simulation setting, the simulator output may itself be a time series, so the estimator can target either the entire trajectory or a scalar functional extracted from it.

For parameterized ordinary differential equations,

$p(y_H)$0

a natural hierarchy is induced by time-integration levels with timesteps $p(y_H)$1 and $p(y_H)$2, giving discrete trajectories $p(y_H)$3 at progressively finer resolutions (Keshavarzzadeh et al., 2018). In this setting, the fidelity notion is numerical rather than statistical: the lower-fidelity models are trajectories computed with larger timesteps.

In uncertainty propagation, the high-fidelity output $p(y_H)$4 and low-fidelity output $p(y_H)$5 are both driven by a random input $p(y_H)$6. The target density satisfies

$p(y_H)$7

so the multifidelity task becomes learning the conditional relationship between high- and low-fidelity responses from a small paired dataset while cheaply sampling the low-fidelity model many times (Nitzler et al., 2020).

For multi-output estimators, each model $p(y_H)$8 may map $p(y_H)$9 to $Y$0, and the target estimator can stack means, variances, and Sobol main effect variances into a single vector-valued quantity (Dixon et al., 2023). For stochastic control and risk quantification, the quantity of interest may be a level-dependent approximation $Y$1 to either a path-integral control law or a safety probability, with trajectory samples generated across a hierarchy of time or state discretizations (Wang et al., 9 Oct 2025).

A plausible implication is that “multifidelity” in trajectory simulation is not tied to a single construction. The literature treats at least three distinct sources of fidelity variation: timestep or discretization level, model-form simplification, and surrogate hyperparameters such as $Y$2 in a coarse simulator (Keshavarzzadeh et al., 2018, Wang et al., 9 Oct 2025, Warner et al., 19 Sep 2025).

2. Linear-surrogate and bandit estimators under a budget

One line of work formulates multifidelity approximation as a linear regression problem over subsets of available surrogates. For any subset $Y$3,

$Y$4

where $Y$5, $Y$6, and $Y$7 is zero-mean noise with variance $Y$8, independent of $Y$9 (Xu et al., 2021). The estimator is built in two phases: an exploration stage with $c_0$0 joint samples $c_0$1 and an exploitation stage in which a chosen surrogate model is sampled alone.

With exploration estimate $c_0$2, the linear-regression Monte Carlo estimator is

$c_0$3

Its conditional mean-squared error is

$c_0$4

and the average conditional MSE has an asymptotic form

$c_0$5

This decomposition separates exploitation error, which scales with the remaining budget, from exploration error, which decays with $c_0$6 (Xu et al., 2021).

The same work casts each subset $c_0$7 as an arm in a modified stochastic bandit. Policies that explore each arm uniformly and then exploit the empirically best arm are termed uniform explore-then-exploit policies. The proposed Adaptive Explore-Then-Commit (AETC) algorithm computes empirical upper bounds

$c_0$8

with

$c_0$9

and commits when further exploration is no longer preferable (Xu et al., 2021).

The theoretical guarantees are explicit. Theorem 3.1 states that under mild sub-Gaussian design assumptions, for fixed $n$0 and large $n$1, with probability at least $n$2,

$n$3

hence

$n$4

Theorem 4.1 states that if $n$5 denote AETC’s exploration length and exploitation model, then almost surely as $n$6,

$n$7

where $n$8 are the oracle choices minimizing the leading-order MSE. The paper therefore describes AETC as “trajectory-wise optimal” (Xu et al., 2021).

The vector-valued extension replaces scalar noise by $n$9 in

$X^{(i)}$0

and measures error through a quadratic $X^{(i)}$1-risk,

$X^{(i)}$2

The resulting bounds depend only on scalar traces such as $X^{(i)}$3 and $X^{(i)}$4, so the approach avoids estimating high-dimensional parameters directly (Xu et al., 2021).

A common misconception is that multifidelity estimation requires an ordered fidelity hierarchy or prior covariance knowledge. This framework states the opposite: it does not assume that fidelity $X^{(i)}$5 is strictly cheaper or less accurate than $X^{(i)}$6, does not require any nested hierarchy, and does not require prior knowledge of $X^{(i)}$7 or $X^{(i)}$8; these quantities are learned during exploration (Xu et al., 2021).

3. Time-integration hierarchies and convergence acceleration

For parametric ODE trajectories, one three-level construction begins by defining low-, medium-, and high-fidelity models through timesteps

$X^{(i)}$9

with discrete solutions computed by any global-order-$i=1,\dots,n$0 one-step or multistep scheme (Keshavarzzadeh et al., 2018). The exact trajectory is viewed in the Hilbert space $i=1,\dots,n$1, while each discrete trajectory lies in a level-dependent space $i=1,\dots,n$2 equipped with a weighted inner product based on composite Newton–Cotes weights.

The offline stage constructs a reduced basis of size $i=1,\dots,n$3 by a greedy sweep on inexpensive level-1 solutions. Using a large training set $i=1,\dots,n$4, one computes $i=1,\dots,n$5 and selects “important” parameters $i=1,\dots,n$6 via pivoted Cholesky, QR, or the greedy rule

$i=1,\dots,n$7

For each selected $i=1,\dots,n$8, the expensive trajectories $i=1,\dots,n$9 and $c_i$0 are then computed (Keshavarzzadeh et al., 2018).

For a new parameter $c_i$1, coefficients $c_i$2 are obtained from the level-1 least-squares problem

$c_i$3

or equivalently by the normal equations

$c_i$4

The multifidelity surrogates are

$c_i$5

Under standard $c_i$6 global-truncation and quadrature-accuracy assumptions, Theorem 3.1 gives

$c_i$7

so the $c_i$8-dependent term decays like $c_i$9 (Keshavarzzadeh et al., 2018).

The same paper then adds a three-step acceleration procedure. First, the discrete surrogates $\sum_{i=0}^n c_i\,n_i \le B.$0 are built at all three levels. Second, each level is lifted to continuous time by a global $\sum_{i=0}^n c_i\,n_i \le B.$1th-order B-spline interpolant $\sum_{i=0}^n c_i\,n_i \le B.$2 defined on a normalized knot vector and fitted by interpolation at the discrete time nodes. Third, Richardson extrapolation is applied pointwise in level index $\sum_{i=0}^n c_i\,n_i \le B.$3. For fixed $\sum_{i=0}^n c_i\,n_i \le B.$4 one sets $\sum_{i=0}^n c_i\,n_i \le B.$5, estimates the observed order

$\sum_{i=0}^n c_i\,n_i \le B.$6

and computes

$\sum_{i=0}^n c_i\,n_i \le B.$7

The resulting accelerated continuous surrogate $\sum_{i=0}^n c_i\,n_i \le B.$8 is reported to be practically $\sum_{i=0}^n c_i\,n_i \le B.$9-accurate in many cases (Keshavarzzadeh et al., 2018).

The numerical examples establish the trajectory-specific character of the method. For the damped harmonic oscillator with $\mu:=\mathbb{E}[f(Y)]$0, $\mu:=\mathbb{E}[f(Y)]$1, $\mu:=\mathbb{E}[f(Y)]$2, $\mu:=\mathbb{E}[f(Y)]$3, $\mu:=\mathbb{E}[f(Y)]$4, and methods RK2, RK3, RK4 and AB2, AB3, AB4, the estimated orders $\mu:=\mathbb{E}[f(Y)]$5 recover the nominal method orders, the spline-accelerated errors $\mu:=\mathbb{E}[f(Y)]$6 attain their minima at the predicted Richardson weights, and log-error versus $\mu:=\mathbb{E}[f(Y)]$7 shows an approximate $\mu:=\mathbb{E}[f(Y)]$8 slope (Keshavarzzadeh et al., 2018). When $\mu:=\mathbb{E}[f(Y)]$9 and 1,000 Monte Carlo samples are used to estimate statistical moments, direct high-fidelity moments have cost $p(y_H)$00, whereas trifidelity moments via $p(y_H)$01 have cost approximately $p(y_H)$02 and achieve an order-of-magnitude lower $p(y_H)$03 error in mean and standard deviation. In the predator–prey example, the accelerated surrogate again shows roughly $p(y_H)$04 decay in $p(y_H)$05 for AB2, AB3, and RK4 (Keshavarzzadeh et al., 2018).

4. Bayesian multifidelity Monte Carlo for trajectory uncertainty propagation

A different construction begins from the identity

$p(y_H)$06

and learns the conditional density from a small paired dataset

$p(y_H)$07

In the Bayesian multi-fidelity Monte-Carlo framework, the conditional relation is modeled as

$p(y_H)$08

with a Gaussian-process prior

$p(y_H)$09

At a new low-fidelity value $p(y_H)$10, the posterior predictive distribution is

$p(y_H)$11

(Nitzler et al., 2020).

For any functional $p(y_H)$12,

$p(y_H)$13

Replacing the unknown conditional by the GP posterior estimate yields the BMFMC estimator

$p(y_H)$14

For the mean, $p(y_H)$15, this reduces to

$p(y_H)$16

Its bias is

$p(y_H)$17

and the variance is approximated by

$p(y_H)$18

with mean-squared error $p(y_H)$19 (Nitzler et al., 2020).

The practical algorithm separates cheap and expensive work. First, one draws $p(y_H)$20 low-fidelity samples. Second, one adds one or two low-dimensional features $p(y_H)$21 at no extra cost and forms

$p(y_H)$22

Third, one selects $p(y_H)$23 training indices by a space-filling design in $p(y_H)$24-space, evaluates the high-fidelity model there, and fits the GP $p(y_H)$25. Fourth, one predicts posterior means and variances at all low-fidelity samples and forms the Monte Carlo average (Nitzler et al., 2020).

The error analysis distinguishes model error, due to a misspecified conditional family, from epistemic uncertainty, due to finite $p(y_H)$26. A key inequality is

$p(y_H)$27

showing that informative additional features reduce the average conditional variance. The computational cost of GP fitting is $p(y_H)$28 and prediction at $p(y_H)$29 points costs $p(y_H)$30, with total work dominated by cheap low-fidelity evaluations plus $p(y_H)$31 high-fidelity runs (Nitzler et al., 2020).

For trajectory-type outputs, the framework allows at least three representations: applying BMFMC to scalar trajectory statistics such as maxima or time averages, concatenating observations at a small set of time-points into a multi-output GP, or embedding time through delays via feature vectors of the form $p(y_H)$32 (Nitzler et al., 2020). The paper also states that credible intervals from the GP posterior predictive variance can be used to certify reliability, and that overly wide intervals motivate additional high-fidelity training points or richer features.

The numerical demonstrations are directly trajectory-oriented. In stochastic flow past a cylinder, the high-fidelity model is high-order DG with $p(y_H)$33, the low-fidelity model is $p(y_H)$34 and is approximately $p(y_H)$35 cheaper, the random inputs include cylinder radius, position, viscosity, and an uncertain-inflow Gaussian field of dimension approximately $p(y_H)$36 reduced to $p(y_H)$37 via Karhunen–Loève, and $p(y_H)$38 low-fidelity runs with $p(y_H)$39 high-fidelity runs are sufficient for BMFMC to recover the full PDF of the maximum lift coefficient to within negligible Kullback–Leibler divergence of order $p(y_H)$40, with overall speed-up approximately $p(y_H)$41 (Nitzler et al., 2020). In the bending-wall fluid–structure interaction problem, the quantity of interest is the tip deflection trajectory $p(y_H)$42 or its peak magnitude, $p(y_H)$43, $p(y_H)$44, and the cheapest low-fidelity model is $p(y_H)$45 cheaper; the reported speed-up is greater than $p(y_H)$46 while credible intervals reliably bracket the reference (Nitzler et al., 2020).

5. Control variates, multi-output covariance structure, and multilevel–multifidelity fusion

Approximate control variates form another major class of multifidelity trajectory estimators. In a general linear estimator,

$p(y_H)$47

$p(y_H)$48 is the high-fidelity quantity of interest, $p(y_H)$49 are parametrized low-fidelity models, and the total cost is

$p(y_H)$50

(Warner et al., 19 Sep 2025). For fixed sample counts, the optimal weights are

$p(y_H)$51

where $p(y_H)$52 and $p(y_H)$53, giving

$p(y_H)$54

An equivalent explicit form expresses the variance through $p(y_H)$55, correlation vectors $p(y_H)$56, correlation matrices $p(y_H)$57, and sample-allocation matrices $p(y_H)$58 and $p(y_H)$59 (Warner et al., 19 Sep 2025).

For multi-output and multi-statistic estimation, the covariance structure becomes the central object. If each model output is $p(y_H)$60, then the literature derives closed-form covariance expressions for mean–mean, variance–variance, mean–variance, main-effect Sobol, and variance–main-effect estimators. Denoting the high-fidelity stacked estimator by $p(y_H)$61 and the stacked low-fidelity differences by

$p(y_H)$62

the multi-output ACV estimator is

$p(y_H)$63

with variance

$p(y_H)$64

Minimizing the determinant or trace yields

$p(y_H)$65

and therefore

$p(y_H)$66

(Dixon et al., 2023).

The sample-allocation problem is then posed under costs $p(y_H)$67 as

$p(y_H)$68

with the objective of minimizing $p(y_H)$69 over integer sample counts (Dixon et al., 2023). In simple scalar cases, the ratio

$p(y_H)$70

shows explicitly that cheap and highly correlated low-fidelity models should be sampled more heavily (Dixon et al., 2023).

A closely related but distinct construction fuses multilevel Monte Carlo and multi-fidelity Monte Carlo. With level-dependent approximations $p(y_H)$71 to a quantity of interest $p(y_H)$72, MLMC uses the telescoping identity

$p(y_H)$73

whereas MFMC combines model outputs through control variates. The integrated MLMF estimator is

$p(y_H)$74

It can also be rewritten as

$p(y_H)$75

(Wang et al., 9 Oct 2025).

The theory for MLMF is explicit. The unbiasedness theorem states that $p(y_H)$76. Under the optimal coefficients

$p(y_H)$77

the variance becomes

$p(y_H)$78

with $p(y_H)$79. Under the stated assumptions of descending $p(y_H)$80 and increasing costs $p(y_H)$81, the estimator is unbiased, consistent, and achieves RMS error $p(y_H)$82 at cost $p(y_H)$83 (Wang et al., 9 Oct 2025).

These two strands are complementary rather than competing. One emphasizes closed-form covariance blocks for complex stacked trajectory statistics; the other emphasizes unbiased telescoping across levels and models with provable $p(y_H)$84 cost scaling (Dixon et al., 2023, Wang et al., 9 Oct 2025).

6. Automated tuning, empirical domains, and recurrent limitations

The effectiveness of multifidelity trajectory estimators depends on the relative correlations and computational costs of the available models, and one recent direction treats low-fidelity hyperparameters themselves as optimization variables (Warner et al., 19 Sep 2025). In that formulation, the target is to choose both model hyperparameters $p(y_H)$85 and sample allocations $p(y_H)$86 to minimize $p(y_H)$87 subject to a cost constraint: $p(y_H)$88 The equivalent penalized form uses a Lagrangian

$p(y_H)$89

with KKT conditions

$p(y_H)$90

(Warner et al., 19 Sep 2025).

The practical algorithm is bi-level. The inner loop solves the multifidelity sample-allocation problem for fixed $p(y_H)$91, while the outer loop solves for $p(y_H)$92 using Efficient Global Optimization, that is, Bayesian optimization with a Gaussian-process surrogate and Expected Improvement (Warner et al., 19 Sep 2025). The online workflow evaluates pilot samples to estimate $p(y_H)$93, $p(y_H)$94, and $p(y_H)$95, then solves an approximate control variates allocation for the remaining budget. Three pilot modes are discussed—offline, projection, and online—and the paper uses the online mode to hedge pilot cost against accuracy (Warner et al., 19 Sep 2025).

The reported trajectory application is entry, descent, and landing for the ADEPT SR-1 flight test using NASA’s POST2 code. The uncertain inputs comprise approximately $p(y_H)$96 variables. Four models are used: full POST2 with $p(y_H)$97 and cost $p(y_H)$98; reduced-physics POST2 at approximately $p(y_H)$99; coarse-timestep POST2 at approximately $Y$00; and a machine-learning surrogate at approximately $Y$01 (Warner et al., 19 Sep 2025). The tunable models vary over $Y$02, and the correlation-versus-cost sweep shows that as $Y$03 increases, cost decreases monotonically while correlations can be nonmonotonic, including a sharp drop around $Y$04 (Warner et al., 19 Sep 2025). This directly supports the paper’s recommendation that hand-tuning $Y$05 is unreliable and that models that look “very inaccurate” may still be useful if correlations remain high.

The budget dependence of tuning overhead is also explicit. For single-parameter tuning of the time-of-flight quantity of interest, with budgets $Y$06, pilot sizes $Y$07, and EGO iterations $Y$08 across $Y$09 random trials, tuned ACV lies between hand-tuned and oracle ACV for all budgets and strictly beats hand-tuned ACV at $Y$10 (Warner et al., 19 Sep 2025). For two-parameter tuning up to $Y$11, the tuning overhead falls below $Y$12 once $Y$13, and for $Y$14 tuned estimators consistently outperform hand-tuned ones and approach the oracle solution; by contrast, at low $Y$15, heavy tuning with large $Y$16 or $Y$17 can degrade performance (Warner et al., 19 Sep 2025). A parametric ACV enhancement based on a search over $Y$18 admissible recursion trees in the Generalized Multifidelity framework yields up to $Y$19 more variance reduction for the landing-latitude quantity of interest at $Y$20, $Y$21 (Warner et al., 19 Sep 2025).

The diversity of empirical domains is summarized below.

Setting	Estimator mechanism	Reported outcome
Parametric ODE moments	Three-level surrogate + spline lifting + Richardson extrapolation (Keshavarzzadeh et al., 2018)	Cost $Y$22 vs $Y$23 with order-of-magnitude lower $Y$24 error
Flow past a cylinder	BMFMC with DG $Y$25 HF and $Y$26 LF (Nitzler et al., 2020)	KL divergence of order $Y$27 with $Y$28 HF runs; speed-up $Y$29
Bending-wall FSI	BMFMC on tip deflection trajectory or peak (Nitzler et al., 2020)	$Y$30 HF runs; speed-up $Y$31
5-dimensional Ornstein–Uhlenbeck risk estimation	MLMF with feature-level hierarchy (Wang et al., 9 Oct 2025)	$Y$32 lower cost at the same error
Path-integral control	MLMF with surrogate dynamics (Wang et al., 9 Oct 2025)	Optimal cost recovered at one-third the sampling cost
EDL trajectory UQ	MOACV and tuned ACV/GMF (Dixon et al., 2023, Warner et al., 19 Sep 2025)	$Y$33 variance reduction for touchdown latitude; up to $Y$34 more variance reduction with GMF tuning

Several limitations recur across the literature. Heavy pilot or tuning cost can be detrimental when the total budget is small (Warner et al., 19 Sep 2025). Multi-output ACV requires $Y$35 to be invertible, so redundant statistics may need to be removed or regularized (Dixon et al., 2023). Higher-dimensional hyperparameter spaces are identified as challenging and may require more advanced derivative-free or cost-aware optimizers (Warner et al., 19 Sep 2025). More broadly, a plausible implication is that practical success depends less on any single fidelity hierarchy than on the joint management of correlation structure, estimator covariance, and budget partitioning across exploration, training, and exploitation.