MF-GLaMs: Multifidelity Lambda Models

• MF-GLaMs models the full conditional response distribution using a flexible four-parameter generalized lambda distribution.
• It couples low-fidelity and high-fidelity simulators through parameter-level discrepancy functions to correct bias.
• The approach leverages polynomial chaos expansions and robust optimization to efficiently estimate and validate surrogate parameters.

Searching arXiv for the MF-GLaM paper and closely related context. arXiv search query: "MF-GLaM A multifidelity stochastic emulator using generalized lambda models" Multifidelity Generalized Lambda Models (MF-GLaMs) are a non-intrusive multifidelity surrogate modeling approach for stochastic simulators that aims to emulate the full conditional response distribution of a high-fidelity (HF) simulator by exploiting data from a lower-fidelity (LF) stochastic simulator. In this setting, a simulator produces random outputs even at fixed input conditions because of unobservable, uncontrollable, or unmodeled input variables, so the target is not a deterministic response surface but the conditional distribution $p_H(y \mid x)$. MF-GLaMs represent that distribution parametrically through a four-parameter generalized lambda distribution (GLD) whose parameters vary with the input, and then couple LF and HF models through parameter-level discrepancies. The method is designed for black-box access, requires only input/output pairs, and does not require replications at identical inputs or access to internal random seeds (Giannoukou et al., 14 Jul 2025).

1. Problem setting and modeling objective

For inputs $x \in \mathbb{R}^d$, a stochastic simulator $M_s$ produces random outputs $Y$, so that $Y \mid x$ has a conditional distribution $p(y \mid x)$. MF-GLaMs consider two such simulators: a high-fidelity simulator with conditional distribution $p_H(y \mid x)$ and a low-fidelity simulator with conditional distribution $p_L(y \mid x)$. The objective is to build a surrogate that predicts the full HF conditional distribution for any $x$, while leveraging LF information.

A defining feature of the framework is its non-intrusive setting. The emulator requires no access to internal random seeds or to the simulator’s stochastic mechanism. Only black-box input/output pairs $\{(x_i,y_i)\}$ are needed, with one random output per input. The training set consists of single realizations $x \in \mathbb{R}^d$0 with independent $x \in \mathbb{R}^d$1, and neither common-random-numbers nor repeated $x \in \mathbb{R}^d$2 at the same $x \in \mathbb{R}^d$3 are required.

This formulation addresses a gap between deterministic multifidelity surrogate modeling and stochastic simulation emulation. While multifidelity surrogate modeling techniques are well-established for deterministic settings, MF-GLaMs target the full conditional response distribution rather than only a mean, variance, or selected quantiles. A plausible implication is that the method is especially relevant when HF runs are expensive and the LF simulator carries enough structure to regularize the estimation of the HF distribution.

2. Generalized lambda representation of conditional distributions

MF-GLaMs build on the generalized lambda model (GLaM), which represents the conditional distribution at each input by a flexible, four-parameter generalized lambda distribution. The formulation uses the FKML parameterization of the GLD (Freimer–Kollia–Mudholkar–Lin). Its quantile function is

$x \in \mathbb{R}^d$4

where $x \in \mathbb{R}^d$5 are input-dependent parameters. For completeness, the RS parameterization, which is not used, is

$x \in \mathbb{R}^d$6

Validity of the induced distribution is enforced through quantile monotonicity. The FKML derivative is

$x \in \mathbb{R}^d$7

Monotonicity holds for all $x \in \mathbb{R}^d$8 if $x \in \mathbb{R}^d$9 and the fractional powers are well-defined and positive. Accordingly, $M_s$0 is the key positivity constraint. For analytic mean and variance expressions, the framework uses $M_s$1.

The induced density can be written through the quantile function. If $M_s$2 is the cdf implied by $M_s$3, then for $M_s$4 in the support,

$M_s$5

Numerically, evaluating the density typically requires solving $M_s$6.

Under the FKML parameterization, for $M_s$7,

$M_s$8

and

$M_s$9

with

$Y$0

$Y$1

The single-fidelity GLaM represents each GLD parameter as a function of $Y$2 using polynomial chaos expansions (PCEs):

$Y$3

$Y$4

$Y$5

The basis functions $Y$6 are multivariate orthogonal polynomials matched to the input marginals; the truncation sets are selected through q-norm or hyperbolic truncation. The exponential mapping for $Y$7 enforces positivity, and because FKML monotonicity is then automatic, the representation yields a valid quantile function under the stated conditions (Giannoukou et al., 14 Jul 2025).

3. Single-fidelity estimation and identifiability

Given observed single draws $Y$8, $Y$9, single-fidelity GLaM estimation proceeds by maximum likelihood. With parameter functions $Y \mid x$0, where $Y \mid x$1, the fitted model solves

$Y \mid x$2

where the density is induced by the FKML quantile representation through

$Y \mid x$3

Initialization and basis construction are structured. PCEs are first fitted to the conditional mean and variance across the single-observation dataset using hybrid LARS, which provides candidate bases $Y \mid x$4 and initial coefficients for $Y \mid x$5 and $Y \mid x$6. The truncation sets for $Y \mid x$7 and $Y \mid x$8 are then selected via a BIC sweep over increasing degrees, and the corresponding coefficients are initialized accordingly. Sparsity is controlled through BIC to prevent overfitting, and initial truncations for $Y \mid x$9 and $p(y \mid x)$0 are obtained from feasible GLS on mean and variance.

Optimization uses a derivative-based trust-region method to maximize the unconstrained likelihood. If constraints or support violations arise, the procedure switches to constrained $p(y \mid x)$1-CMA-ES to recover feasibility. Because density evaluation requires $p(y \mid x)$2, practical implementation depends on robust one-dimensional root finding with good bracketing and monotone quantile guarantees.

The procedure can be summarized in five steps: construct candidate PCE bases and truncation sets using q-norm hyperbolic truncation; initialize $p(y \mid x)$3 from mean and variance fits and $p(y \mid x)$4 from BIC selection; define the log-likelihood; maximize it with trust-region and then constrained CMA-ES if needed; and return the fitted parameter functions and predictive distribution. The framework explicitly states that identifiability is achieved through variation across $p(y \mid x)$5: the four FKML parameters cover location, scale, and shape, and across diverse input samples there is sufficient information to infer the parameter functions even without replications at the same input, provided there is adequate input-space coverage and sample size.

A common misconception is that full distribution emulation for stochastic simulators necessarily requires replicated runs at identical inputs. In MF-GLaMs, the stated assumption is weaker: single observations at varying inputs can be sufficient because identifiability is induced by the functional dependence of the GLD parameters on $p(y \mid x)$6 rather than by repeated sampling at fixed $p(y \mid x)$7.

4. Multifidelity coupling and joint likelihood

The multifidelity extension operates at the parameter level. For each GLD parameter,

$p(y \mid x)$8

where $p(y \mid x)$9 is learned from LF data and $p_H(y \mid x)$0 is a discrepancy function that captures LF-to-HF bias. For $p_H(y \mid x)$1,

$p_H(y \mid x)$2

so that

$p_H(y \mid x)$3

For the scale parameter,

$p_H(y \mid x)$4

which preserves $p_H(y \mid x)$5.

The framework adopts a default simplification for shape discrepancies. Empirically, HF and LF distribution shapes, represented by $p_H(y \mid x)$6 and $p_H(y \mid x)$7, often match closely, so the default choice is

$p_H(y \mid x)$8

hence

$p_H(y \mid x)$9

This reduces model complexity and improves data efficiency. If known shape mismatch exists, low-order $p_L(y \mid x)$0 and $p_L(y \mid x)$1 can be introduced.

Training is based on a weighted joint log-likelihood over LF and HF data. Let $p_L(y \mid x)$2 collect all LF and discrepancy coefficients. For a general LF/HF weight $p_L(y \mid x)$3,

$p_L(y \mid x)$4

The neutral weighting used is $p_L(y \mid x)$5, yielding

$p_L(y \mid x)$6

and estimation is

$p_L(y \mid x)$7

Optimization again uses trust-region and, when constraints or support are violated, constrained $p_L(y \mid x)$8-CMA-ES. Basis selection is split into stages: the LF-only GLaM is trained first to fix the LF truncation sets, then candidate discrepancy truncations for $p_L(y \mid x)$9 and $x$0 are swept over varying degrees and q-norms, and the selected model minimizes

$x$1

The principal assumptions are that LF and HF are correlated so that LF parameters plus a discrepancy can explain HF parameters, and that the discrepancy is simpler than directly learning the HF parameter functions. This suggests that MF-GLaMs are most advantageous when the LF simulator captures the dominant structure of location, scale, and possibly shape, while the HF simulator can be represented through a relatively low-complexity correction (Giannoukou et al., 14 Jul 2025).

5. Computational characteristics, sample allocation, and reported performance

Let $x$2 denote the number of basis functions in each parameter expansion. Each log-likelihood evaluation requires computing the parameter functions as linear combinations of basis terms, at cost $x$3 per sample, plus numerical inversion $x$4 and evaluation of $x$5 for the density, with per-sample inversion cost denoted $x$6. The total per-iteration cost is

$x$7

Scalability increases with input dimension $x$8 through basis size; hyperbolic truncation controls combinatorial growth, and sparsity through BIC and low-order discrepancies improves tractability. Trust-region optimization benefits from derivatives, while CMA-ES improves robustness when constraint handling becomes difficult.

The reported sample-allocation guidance is to use substantially more LF than HF points when LF is cheaper and reasonably correlated with HF. A practical heuristic is to start with $x$9 in the low hundreds, augment with $\{(x_i,y_i)\}$0 an order of magnitude larger when budget permits, and tune $\{(x_i,y_i)\}$1 if HF data are very scarce or abundant.

Three empirical examples are reported. In a synthetic GLD example with $\{(x_i,y_i)\}$2, MF-GLaM with $\{(x_i,y_i)\}$3 and $\{(x_i,y_i)\}$4 achieves approximately one-tenth the error of HF-only at $\{(x_i,y_i)\}$5. More specifically, the median normalized Wasserstein error $\{(x_i,y_i)\}$6 at $\{(x_i,y_i)\}$7 is reported as $\{(x_i,y_i)\}$8 for MF-GLaM versus $\{(x_i,y_i)\}$9 for HF-only on the stated $x \in \mathbb{R}^d$00 scale, with error reduction of about threefold and lower variability. In the borehole example, with HF dimension $x \in \mathbb{R}^d$01 and LF dimension $x \in \mathbb{R}^d$02, the median $x \in \mathbb{R}^d$03 at $x \in \mathbb{R}^d$04 is reported as $x \in \mathbb{R}^d$05 for MF-GLaM versus $x \in \mathbb{R}^d$06 for HF-only on the same scale, and MF-GLaM at $x \in \mathbb{R}^d$07 matches HF-only at $x \in \mathbb{R}^d$08. In the earthquake example, with dependent inputs in $x \in \mathbb{R}^d$09, HF time step $x \in \mathbb{R}^d$10s, LF time step $x \in \mathbb{R}^d$11s, and cost ratio approximately $x \in \mathbb{R}^d$12 for LF:HF, MF-GLaM with $x \in \mathbb{R}^d$13 and $x \in \mathbb{R}^d$14 has median $x \in \mathbb{R}^d$15 versus $x \in \mathbb{R}^d$16 for HF-only, and the MF-GLaM configuration costs about $x \in \mathbb{R}^d$17 HF-equivalent runs compared with $x \in \mathbb{R}^d$18 for HF-only at $x \in \mathbb{R}^d$19, yielding comparable accuracy at about $x \in \mathbb{R}^d$20 lower cost (Giannoukou et al., 14 Jul 2025).

These examples are evaluated with distributional metrics rather than only summary-statistic errors. The reported metrics include the Wasserstein-2 distance,

$x \in \mathbb{R}^d$21

the normalized global error

$x \in \mathbb{R}^d$22

the Kullback–Leibler divergence when tractable,

$x \in \mathbb{R}^d$23

the continuous ranked probability score,

$x \in \mathbb{R}^d$24

and calibration diagnostics based on the probability integral transform and interval coverage.

6. Practical guidance, limitations, and relation to alternative methods

The reported implementation workflow has seven steps. First, collect LF and HF input/output pairs, one output per input. Second, choose PCE inputs assuming independent marginals for basis construction and handle dependence only in sampling; hyperbolic truncation uses $x \in \mathbb{R}^d$25. Third, fit the LF-only GLaM by estimating $x \in \mathbb{R}^d$26 and $x \in \mathbb{R}^d$27 from mean and variance PCEs, selecting $x \in \mathbb{R}^d$28 and $x \in \mathbb{R}^d$29 via BIC, and solving the MLE for the LF coefficients. Fourth, define the multifidelity discrepancy with $x \in \mathbb{R}^d$30 by default and low-degree candidate sets for $x \in \mathbb{R}^d$31 and $x \in \mathbb{R}^d$32, initialized at zero. Fifth, maximize the weighted joint likelihood with $x \in \mathbb{R}^d$33 unless tuned otherwise. Sixth, sweep discrepancy bases and select the model via MF-BIC. Seventh, for prediction at a new input $x \in \mathbb{R}^d$34, compute $x \in \mathbb{R}^d$35 and then obtain the predictive quantile and density through the FKML formulas.

The numerical guidance is specific. The recommended practice is to expand $x \in \mathbb{R}^d$36 and $x \in \mathbb{R}^d$37 with higher degrees than $x \in \mathbb{R}^d$38 and $x \in \mathbb{R}^d$39, to keep $x \in \mathbb{R}^d$40 and $x \in \mathbb{R}^d$41 low-order unless LF-to-HF bias is strong, to enforce $x \in \mathbb{R}^d$42 positivity through the exponential mapping, to constrain $x \in \mathbb{R}^d$43 and $x \in \mathbb{R}^d$44 for finite variance and stability, and to use robust one-dimensional root finding for $x \in \mathbb{R}^d$45 with monotonicity-based bracketing over $x \in \mathbb{R}^d$46. MF-BIC is used to avoid overfitting of the discrepancy, especially when HF data are scarce.

Several limitations are explicit. If LF and HF are weakly correlated, especially when shapes or locations differ markedly, the default low-order discrepancy may be insufficient and higher-order terms or nonzero $x \in \mathbb{R}^d$47 may be necessary; otherwise multifidelity modeling may underperform. Identifiability can weaken in sparse, high-dimensional input spaces, so adequate HF sample size and LF coverage are needed. The FKML GLD is flexible for unimodal, skewed, and heavy-tailed distributions, but multimodal responses are not well captured; mixture extensions or stochastic PCE are suggested as alternatives. Boundary and support issues can arise when $x \in \mathbb{R}^d$48 or $x \in \mathbb{R}^d$49 approach $x \in \mathbb{R}^d$50, causing numerical instability near $x \in \mathbb{R}^d$51 or $x \in \mathbb{R}^d$52. Input dependence is another caveat: PCE bases assume independent inputs, and although independent bases can still give good empirical fits under dependence, extrapolation should be treated cautiously.

The reported comparison to alternative methods is organized by what those methods target. Heteroscedastic Gaussian processes capture mean and variance rather than the full distribution and require replication or explicit noise modeling. Quantile regression fits selected quantiles across $x \in \mathbb{R}^d$53 and needs many quantile levels to reconstruct the full distribution. Mixture density networks are flexible but described as data-hungry, harder to regularize non-intrusively, and potentially brittle without replications. Distribution regression and kernel methods are nonparametric but may require many samples and can be computationally heavy. Co-kriging of moments or quantiles focuses on summary statistics rather than the full conditional distribution. Against these, MF-GLaMs trade parametric stability, low data requirements, and non-intrusiveness against the representational limitation of the FKML GLD family.

A worked one-dimensional example illustrates the construction. With $x \in \mathbb{R}^d$54, LF parameters are specified as

$x \in \mathbb{R}^d$55

and HF discrepancies are

$x \in \mathbb{R}^d$56

Thus

$x \in \mathbb{R}^d$57

With $x \in \mathbb{R}^d$58 LF samples and $x \in \mathbb{R}^d$59 HF samples, using bases $x \in \mathbb{R}^d$60, $x \in \mathbb{R}^d$61, $x \in \mathbb{R}^d$62, LF truncation sets $x \in \mathbb{R}^d$63, $x \in \mathbb{R}^d$64, $x \in \mathbb{R}^d$65, $x \in \mathbb{R}^d$66, and discrepancy truncations $x \in \mathbb{R}^d$67, $x \in \mathbb{R}^d$68 with $x \in \mathbb{R}^d$69, the MF stage estimates $x \in \mathbb{R}^d$70 by maximizing $x \in \mathbb{R}^d$71 with $x \in \mathbb{R}^d$72. At $x \in \mathbb{R}^d$73, the fitted MF parameters approximately match

$x \in \mathbb{R}^d$74

leading to the predictive quantile

$x \in \mathbb{R}^d$75

and derivative

$x \in \mathbb{R}^d$76

with

$x \in \mathbb{R}^d$77

This example illustrates the central mechanism of MF-GLaM: LF data determine a full conditional distribution model, and a small HF sample corrects location and scale through discrepancy terms while inheriting shape from the LF model (Giannoukou et al., 14 Jul 2025).