Functional History Matching (FHM)

• Functional History Matching (FHM) is a calibration approach that treats simulator outputs as functions, using function-space emulation to rule out implausible regions of the input space.
• It employs projection-based comparisons and specialized emulators, such as the Outer Product Emulator and random projections, to manage high-dimensional functional data efficiently.
• Applications of FHM span hydrology, tsunami warning, and building design, with evolving methods that incorporate uncertainty quantification and fuzzy hierarchical calibration.

Functional History Matching (FHM) denotes a family of history-matching methodologies for computer models whose outputs are functions—most commonly time series, but also spatial fields or other high-dimensional functional objects—rather than only scalar outputs. In this usage, FHM preserves the central history-matching objective of ruling out implausible regions of input space by comparing emulator predictions with observations under uncertainty, while replacing scalar comparison by function-space emulation, functional summaries, or projection-based comparisons (Kanai et al., 4 Sep 2025). A distinct usage of the same acronym appears in work on the “Fuzzy Hierarchical Multiplex,” where “Functional History Matching” is embedded in a hierarchical fuzzy optimization framework calibrated to historical targets through an $L_2$ fit (Kafantaris, 10 Dec 2025).

1. Origins, scope, and terminological variation

Classical history matching was developed as a calibration and screening methodology for complex deterministic or stochastic models. Its core mechanism is iterative elimination: an emulator is fitted on a limited set of simulator evaluations, an implausibility measure is defined from emulator prediction and uncertainty, and successive “waves” restrict attention to the Not-Ruled-Out-Yet (NROY) region rather than searching for a single optimum (Drovandi et al., 2017). In this sense, history matching is intentionally weaker than full Bayesian calibration: it identifies inputs that are not yet contradicted by data, rather than a posterior mode or a local best fit.

The transition from scalar history matching to FHM arose because many simulators produce natural outputs as functions of time or space. Hydrological calibration is an early explicit instance of this pressure: rainfall–runoff models return hydrographs, and the calibration target is itself a time series rather than a single number (Bhattacharjee et al., 2017). Related developments extended the same logic to functionals of stochastic outputs, such as overheating risk probabilities or expected annual energy use derived from yearly temperature trajectories in building design (Baker et al., 2020). The term “Functional History Matching” is made explicit in work that treats both simulator outputs and observations as elements of a function space and builds a dedicated functional history-matching framework around that premise (Kanai et al., 4 Sep 2025).

A second terminological strand uses “FHM” differently. In “Fuzzy Hierarchical Multiplex,” the term is essentially synonymous with the multiplex framework itself and denotes the calibration process by which historical input–output patterns are matched in a hierarchical fuzzy network (Kafantaris, 10 Dec 2025). This suggests that the acronym is not fully stabilized across subfields, even though the dominant statistical meaning concerns history matching for functional outputs.

Variant	Functional object	Representative role
Modified HM for hydrological models	Time-series hydrographs	Match selected time points, then minimize full-series distance
History matching for level sets	Functionals of stochastic outputs	Rule out designs failing threshold criteria
Functional History Matching	Time series as functions in $\mathcal{F}$	Functional emulator + projection-based implausibility
Fuzzy Hierarchical Multiplex	Hierarchical fuzzy mappings	Calibrate inner-to-outer metric mappings

2. Statistical formulation of functional implausibility

In the standard history-matching formulation, the simulator output is compared with observations through an implausibility measure that standardizes mismatch by the combined uncertainty from emulator error, observation error, model discrepancy, and, when relevant, intrinsic simulator stochasticity. One canonical scalar form is

$$\mathcal{I}(\theta) = \frac{|y_\theta-y_{\mathrm{obs}}|} {\sqrt{s_{m,\theta}^2+s_{e,\theta}^2+s_d^2+s_{\mathrm{obs}}^2}},$$

with wavewise non-implausible sets obtained by thresholding $\mathcal{I}$ or its multi-output extensions (Drovandi et al., 2017). In Bayes Linear history matching, the same logic appears as

$$I(x)^2=\frac{(E_D[g(x)]-z)^2}{\operatorname{Var}_D[g(x)]+\operatorname{Var}[e]+\operatorname{Var}[\epsilon(x)]},$$

so that emulator uncertainty, observational error, and model discrepancy enter additively in the denominator (Iskauskas et al., 25 Feb 2026).

FHM generalizes this structure from scalar outputs to functional outputs. In the explicit functional formulation, the simulator is

$$f:\Theta\times T\to\mathbb{R}, \qquad (\theta,t)\mapsto f(\theta,t),$$

or equivalently $f:\Theta\to\mathcal{F}$, where $\mathcal{F}$ is a function space. Observations satisfy

$$z(t)=y(t)+\phi(t), \qquad y(t)=f(\theta)(t)+\delta(t),$$

and with an emulator $f_{\text{emu}}$ the working decomposition becomes

$\mathcal{F}$0

The functional problem is therefore still a problem of screening inputs $\mathcal{F}$1, but the object being matched is now a curve rather than a vector of a few summary statistics (Kanai et al., 4 Sep 2025).

A distinctive feature of the functional framework is that implausibility is defined after projecting curves into scalar spaces. Let $\mathcal{F}$2 be a Hilbert space with inner product

$\mathcal{F}$3

For a random functional direction $\mathcal{F}$4, the observed and emulated curves are projected as $\mathcal{F}$5 and $\mathcal{F}$6, and the projection-wise functional implausibility is

$\mathcal{F}$7

where $\mathcal{F}$8 and $\mathcal{F}$9 contains projected emulator uncertainty together with observation and discrepancy terms (Kanai et al., 4 Sep 2025). In that work, because projected functional distributions may be non-unimodal, the implausibility threshold is taken as $$\mathcal{I}(\theta) = \frac{|y_\theta-y_{\mathrm{obs}}|} {\sqrt{s_{m,\theta}^2+s_{e,\theta}^2+s_d^2+s_{\mathrm{obs}}^2}},$$0 by invoking Chebyshev’s inequality rather than the classical $$\mathcal{I}(\theta) = \frac{|y_\theta-y_{\mathrm{obs}}|} {\sqrt{s_{m,\theta}^2+s_{e,\theta}^2+s_d^2+s_{\mathrm{obs}}^2}},$$1 rule.

Other FHM-style formulations achieve the same transition to functional outputs by reduction rather than direct function-space matching. Hydrological calibration selects a small discretization-point set $$\mathcal{I}(\theta) = \frac{|y_\theta-y_{\mathrm{obs}}|} {\sqrt{s_{m,\theta}^2+s_{e,\theta}^2+s_d^2+s_{\mathrm{obs}}^2}},$$2, emulates each scalar mapping $$\mathcal{I}(\theta) = \frac{|y_\theta-y_{\mathrm{obs}}|} {\sqrt{s_{m,\theta}^2+s_{e,\theta}^2+s_d^2+s_{\mathrm{obs}}^2}},$$3, and combines the resulting pointwise implausibilities through

$$\mathcal{I}(\theta) = \frac{|y_\theta-y_{\mathrm{obs}}|} {\sqrt{s_{m,\theta}^2+s_{e,\theta}^2+s_d^2+s_{\mathrm{obs}}^2}},$$4

before choosing the best parameter vector by minimizing the full time-series distance $$\mathcal{I}(\theta) = \frac{|y_\theta-y_{\mathrm{obs}}|} {\sqrt{s_{m,\theta}^2+s_{e,\theta}^2+s_d^2+s_{\mathrm{obs}}^2}},$$5 (Bhattacharjee et al., 2017). Building-design history matching similarly treats simulator outputs as functional in time and stochastic weather, but history matches functionals of those outputs, notably the overheating probability and expected annual heating energy, through threshold-based implausibility rather than direct curve matching (Baker et al., 2020).

3. Emulation strategies for functional outputs

The principal technical difficulty in FHM is emulation: a naive GP treatment of a finely discretized time series as a very high-dimensional vector is statistically and computationally inefficient. The most explicit solution is the Outer Product Emulator (OPE), a GP-based emulator for deterministic functional codes that models $$\mathcal{I}(\theta) = \frac{|y_\theta-y_{\mathrm{obs}}|} {\sqrt{s_{m,\theta}^2+s_{e,\theta}^2+s_d^2+s_{\mathrm{obs}}^2}},$$6 jointly over input and time through a regression-plus-residual decomposition

$$\mathcal{I}(\theta) = \frac{|y_\theta-y_{\mathrm{obs}}|} {\sqrt{s_{m,\theta}^2+s_{e,\theta}^2+s_d^2+s_{\mathrm{obs}}^2}},$$7

where the regressors are outer products

$$\mathcal{I}(\theta) = \frac{|y_\theta-y_{\mathrm{obs}}|} {\sqrt{s_{m,\theta}^2+s_{e,\theta}^2+s_d^2+s_{\mathrm{obs}}^2}},$$8

In the tsunami application, the input basis functions are Legendre polynomials up to degree 2 in each input dimension, the time basis functions are low-order Fourier terms up to order 2, and the residual covariance is separable in input and time (Kanai et al., 4 Sep 2025).

The residual covariance in OPE is written as

$$\mathcal{I}(\theta) = \frac{|y_\theta-y_{\mathrm{obs}}|} {\sqrt{s_{m,\theta}^2+s_{e,\theta}^2+s_d^2+s_{\mathrm{obs}}^2}},$$9

with componentwise kernels

$\mathcal{I}$0

$\mathcal{I}$1

A conjugate NIG prior is then used for regression coefficients, residual GP, and variance scaling, permitting analytical posterior prediction of mean curves and covariance functions (Kanai et al., 4 Sep 2025).

An alternative emulation strategy is scalarization. In hydrological calibration, one GP emulator is built for each selected time point $\mathcal{I}$2, with predictive mean

$\mathcal{I}$3

and predictive variance

$\mathcal{I}$4

This avoids functional GP construction by reducing the hydrograph to a small number of informative landmarks (Bhattacharjee et al., 2017).

Further reductions appear in applications where the functional object is represented by selected bins or scalar functionals. In non-perturbative event-generator calibration, the outputs are distributions and spectra represented by 432 histogram bins; principal variable analysis is used to select 50–70 bins per wave such that at least 95% of the variability per observable is captured, and Bayes Linear emulators are fitted only to those selected scalar components (Iskauskas et al., 25 Feb 2026). In building-design screening, the temperature time series is reduced to a CIBSE overheating classification and annual heating energy usage, and the emulators are correspondingly a GP classifier for the latent logit overheating probability and a heteroscedastic GP for mean energy usage (Baker et al., 2020).

Random projection supplies a second layer of reduction in explicit FHM. Independent directions $\mathcal{I}$5 are sampled from a Wiener process on the time grid, and each curve is mapped to

$\mathcal{I}$6

The stated motivation is that random projections preserve enough information about distances and distributions to compare curves in a low-dimensional space while also accommodating derivatives, so that both waveform shape and phase information contribute to screening (Kanai et al., 4 Sep 2025).

4. Wave-based algorithms, NROY geometry, and uncertainty quantification

FHM inherits the multi-wave structure of classical history matching. At each wave, a design is generated in the current admissible region, the simulator is evaluated there, emulators are fitted or updated, implausibility is computed over a much larger candidate set, and the NROY region is redefined by thresholding. In hydrological calibration, the plausible set at iteration $\mathcal{I}$7 is

$\mathcal{I}$8

where $\mathcal{I}$9 is a space-filling test set of size $$I(x)^2=\frac{(E_D[g(x)]-z)^2}{\operatorname{Var}_D[g(x)]+\operatorname{Var}[e]+\operatorname{Var}[\epsilon(x)]},$$0, and the algorithm stops when the new plausible set is empty; among all evaluated points, the estimated inverse solution is

$$I(x)^2=\frac{(E_D[g(x)]-z)^2}{\operatorname{Var}_D[g(x)]+\operatorname{Var}[e]+\operatorname{Var}[\epsilon(x)]},$$1

over the full time series (Bhattacharjee et al., 2017). In the tsunami framework, wave 2 restricts the parameter space to the hyperrectangle defined by the min/max of each parameter in $$I(x)^2=\frac{(E_D[g(x)]-z)^2}{\operatorname{Var}_D[g(x)]+\operatorname{Var}[e]+\operatorname{Var}[\epsilon(x)]},$$2, a new Latin hypercube is sampled within that restricted region, and the emulator is retrained before recomputing functional implausibility (Kanai et al., 4 Sep 2025).

The geometry of the non-implausible region is often highly irregular. The SMC formulation makes this explicit by defining a target sequence

$$I(x)^2=\frac{(E_D[g(x)]-z)^2}{\operatorname{Var}_D[g(x)]+\operatorname{Var}[e]+\operatorname{Var}[\epsilon(x)]},$$3

which is the prior restricted to the current non-implausible region. This work shows that non-implausible regions can be multi-modal, highly irregular, and very difficult to sample uniformly, especially once outputs become rich or functional (Drovandi et al., 2017). In that semi-automatic formulation, the next threshold $$I(x)^2=\frac{(E_D[g(x)]-z)^2}{\operatorname{Var}_D[g(x)]+\operatorname{Var}[e]+\operatorname{Var}[\epsilon(x)]},$$4 is chosen as the empirical $$I(x)^2=\frac{(E_D[g(x)]-z)^2}{\operatorname{Var}_D[g(x)]+\operatorname{Var}[e]+\operatorname{Var}[\epsilon(x)]},$$5-quantile of current particle implausibilities so that approximately $$I(x)^2=\frac{(E_D[g(x)]-z)^2}{\operatorname{Var}_D[g(x)]+\operatorname{Var}[e]+\operatorname{Var}[\epsilon(x)]},$$6 particles survive, thereby automating a major tuning choice in practical history matching.

Uncertainty decomposition remains central throughout. In scalar HM and its functional extensions, the denominator of the implausibility measure aggregates emulator uncertainty, observation error, and model discrepancy, with stochastic simulator variance added when required (Drovandi et al., 2017). In explicit FHM, the projected variance is

$$I(x)^2=\frac{(E_D[g(x)]-z)^2}{\operatorname{Var}_D[g(x)]+\operatorname{Var}[e]+\operatorname{Var}[\epsilon(x)]},$$7

where the projected emulator variance is computed from the OPE predictive covariance and the observation-plus-discrepancy term is estimated through a cross-validation procedure that inflates a fixed additional uncertainty until a specified fraction of curves is non-implausible under a worst-case scenario (Kanai et al., 4 Sep 2025). In building-design level-set estimation, model discrepancy is not used in the case study, but the paper states directly that standard history matching would incorporate it by replacing $$I(x)^2=\frac{(E_D[g(x)]-z)^2}{\operatorname{Var}_D[g(x)]+\operatorname{Var}[e]+\operatorname{Var}[\epsilon(x)]},$$8 with $$I(x)^2=\frac{(E_D[g(x)]-z)^2}{\operatorname{Var}_D[g(x)]+\operatorname{Var}[e]+\operatorname{Var}[\epsilon(x)]},$$9 in the implausibility denominator (Baker et al., 2020).

For high-dimensional functional outputs, uncertainty handling may also be componentwise. Event-generator calibration treats Monte Carlo variability as an $$f:\Theta\times T\to\mathbb{R}, \qquad (\theta,t)\mapsto f(\theta,t),$$0-independent model discrepancy term because the simulator is stochastic and, in about 8% of bins, this variability exceeds observational error, sometimes by a factor of 5; additional structural discrepancy proportional to the magnitude of the observation is also included (Iskauskas et al., 25 Feb 2026). A plausible implication is that functional history matching becomes practically reliable only when discrepancy and stochastic noise are represented at the same resolution as the functional summaries being matched.

5. Applications and reported empirical performance

Applications of FHM and FHM-like methodology span hydrology, building design, tsunami warning, and high-energy event-generator calibration. In hydrology, the modified history-matching procedure was applied to a Matlab–Simulink runoff model with a 5445-point hydrograph and a SWAT watershed model with 84 monthly observations (Bhattacharjee et al., 2017). For the Matlab–Simulink case, the reported metrics changed from RMSE $$f:\Theta\times T\to\mathbb{R}, \qquad (\theta,t)\mapsto f(\theta,t),$$1 to $$f:\Theta\times T\to\mathbb{R}, \qquad (\theta,t)\mapsto f(\theta,t),$$2, from $$f:\Theta\times T\to\mathbb{R}, \qquad (\theta,t)\mapsto f(\theta,t),$$3 to $$f:\Theta\times T\to\mathbb{R}, \qquad (\theta,t)\mapsto f(\theta,t),$$4, from NSE $$f:\Theta\times T\to\mathbb{R}, \qquad (\theta,t)\mapsto f(\theta,t),$$5 to $$f:\Theta\times T\to\mathbb{R}, \qquad (\theta,t)\mapsto f(\theta,t),$$6, from PPTS$$f:\Theta\times T\to\mathbb{R}, \qquad (\theta,t)\mapsto f(\theta,t),$$7 to $$f:\Theta\times T\to\mathbb{R}, \qquad (\theta,t)\mapsto f(\theta,t),$$8, and from PPTS$$f:\Theta\times T\to\mathbb{R}, \qquad (\theta,t)\mapsto f(\theta,t),$$9 to $f:\Theta\to\mathcal{F}$0. For SWAT, the reported change was from RMSE $f:\Theta\to\mathcal{F}$1 to $f:\Theta\to\mathcal{F}$2, from $f:\Theta\to\mathcal{F}$3 to $f:\Theta\to\mathcal{F}$4, from NSE $f:\Theta\to\mathcal{F}$5 to $f:\Theta\to\mathcal{F}$6, from PPTS$f:\Theta\to\mathcal{F}$7 to $f:\Theta\to\mathcal{F}$8, and from PPTS$f:\Theta\to\mathcal{F}$9 to $\mathcal{F}$0.

In building design under climate change, the level-set history-matching formulation targets the set

$\mathcal{F}$1

with overheating probability below $\mathcal{F}$2 and mean heating energy usage below $\mathcal{F}$3 kWh/m$\mathcal{F}$4. Starting from 1,000,000 candidate designs, the reported proportions were: wave 1 NROY $\mathcal{F}$5 and tenable ($\mathcal{F}$6) $\mathcal{F}$7; wave 2 NROY $\mathcal{F}$8 and tenable $\mathcal{F}$9; wave 3 NROY $$z(t)=y(t)+\phi(t), \qquad y(t)=f(\theta)(t)+\delta(t),$$0, tenable $$z(t)=y(t)+\phi(t), \qquad y(t)=f(\theta)(t)+\delta(t),$$1, and ruled-in ($$z(t)=y(t)+\phi(t), \qquad y(t)=f(\theta)(t)+\delta(t),$$2) $$z(t)=y(t)+\phi(t), \qquad y(t)=f(\theta)(t)+\delta(t),$$3 (Baker et al., 2020). The methodology is history matching “in spirit,” but here the matched objects are functionals of stochastic yearly trajectories rather than direct time-series curves.

The tsunami-warning application is the most direct operational realization of FHM. Using wave elevation time series from DART buoy data to infer tsunami source parameters, the method applies OPE plus 1000 random Wiener projections, works jointly with original waveforms and their derivatives, and uses an implausibility cutoff of $$z(t)=y(t)+\phi(t), \qquad y(t)=f(\theta)(t)+\delta(t),$$4 (Kanai et al., 4 Sep 2025). In the comparison at Mumbai 001, wave 2 functional HM yielded an NROY ensemble of $$z(t)=y(t)+\phi(t), \qquad y(t)=f(\theta)(t)+\delta(t),$$5 emulated curves, whereas landmark-based scalar HM yielded $$z(t)=y(t)+\phi(t), \qquad y(t)=f(\theta)(t)+\delta(t),$$6 curves. The reported envelope amplitude was $$z(t)=y(t)+\phi(t), \qquad y(t)=f(\theta)(t)+\delta(t),$$7 for FHM and $$z(t)=y(t)+\phi(t), \qquad y(t)=f(\theta)(t)+\delta(t),$$8 for landmark HM, and the landmark-based method was reported to include outlier curves with qualitatively incorrect shapes. The same study reports a JAGURS runtime of about $$z(t)=y(t)+\phi(t), \qquad y(t)=f(\theta)(t)+\delta(t),$$9 s per simulation, OPE training time of about $f_{\text{emu}}$0 s for 100 runs and 5 parameters, emulation time of about $f_{\text{emu}}$1 s per curve with derivative, and random-projection time of about $f_{\text{emu}}$2 s.

In event-generator calibration, the functional object is not a time series but a set of distributions and spectra. The calibration uses 29 observables represented by 432 bins and Bayes Linear emulation across 3 waves for the cluster model and 5 waves for the string model (Iskauskas et al., 25 Feb 2026). Reported non-implausible volume reductions are extreme: for Ahadic, from $f_{\text{emu}}$3 initially to about $f_{\text{emu}}$4 after wave 1, $f_{\text{emu}}$5 after wave 2, and $f_{\text{emu}}$6 after wave 3; for Lund string, from $f_{\text{emu}}$7 to $f_{\text{emu}}$8, $f_{\text{emu}}$9, $\mathcal{F}$00, $\mathcal{F}$01, and $\mathcal{F}$02 across waves 1–5, with a final test at $\mathcal{F}$03. These results are presented as evidence that history matching can retain all non-implausible regions, including disjoint or curved ones, rather than collapsing the problem to a single local tune.

6. The Fuzzy Hierarchical Multiplex interpretation

In “Fuzzy Hierarchical Multiplex,” the term “Functional History Matching” is used in a markedly different way. The paper states that the term is essentially synonymous with the Fuzzy Hierarchical Multiplex itself and is “not defined as a separate module but is embedded in how the multiplex is trained and used” (Kafantaris, 10 Dec 2025). The framework extends Fuzzy Cognitive Maps (FCMs) by replacing a flat graph of concepts $\mathcal{F}$04 and causal weights $\mathcal{F}$05 with a hierarchical multiplex whose inner layer consists of subsystem nodes $\mathcal{F}$06 and whose outer layer consists of system-level metrics $\mathcal{F}$07.

The basic update structure is

$\mathcal{F}$08

where $\mathcal{F}$09 is the inner multiplex update function and $\mathcal{F}$10 is the outer multiplex objective function. The history-matching aspect consists in fitting historical data so that the multiplex outputs match target metrics $\mathcal{F}$11 as closely as possible, explicitly through minimizing an $\mathcal{F}$12 norm. The paper gives the example residual

$\mathcal{F}$13

and states the corresponding loss functional as

$\mathcal{F}$14

or more generally

$\mathcal{F}$15

This framework is explicitly fuzzy: input data, internal activations, and output metrics are treated as fuzzy quantities with membership in $\mathcal{F}$16, and the paper states that metric predictions “would result in an Interval valued Fuzzy Set (IVFS) for the metrics.” The intended application is service process optimization, with outer metrics including Wait, Throughput, Utilization, Patience, and a composite Contribution measure. The illustrative table reports, for five nodes, Wait values from $\mathcal{F}$17 to $\mathcal{F}$18, Throughput from $\mathcal{F}$19 to $\mathcal{F}$20, Utilization from $\mathcal{F}$21 to $\mathcal{F}$22, Patience from $\mathcal{F}$23 to $\mathcal{F}$24, and Contribution values $\mathcal{F}$25, $\mathcal{F}$26, $\mathcal{F}$27, $\mathcal{F}$28, and $\mathcal{F}$29 (Kafantaris, 10 Dec 2025).

This use of “Functional History Matching” differs from the statistical meaning of FHM in several respects. Here, “history” refers to observed sequences of service-process inputs and metric outcomes, “matching” refers to gradient-based calibration of a hierarchical fuzzy mapping, and “functional” refers to the learned mapping from inner configurations to outer metrics rather than to function-valued simulator outputs. A plausible implication is that the two traditions share an emphasis on calibrating mappings to historical behavior, but not the same mathematical object.

7. Limitations, methodological tensions, and future directions

The current FHM literature identifies several recurring limitations. In explicit functional HM, model discrepancy is handled pragmatically as a fixed additional variance term estimated by cross-validation, and the implausibility threshold of $\mathcal{F}$30 is justified by Chebyshev’s inequality; the authors note that the separable OPE covariance and low-order basis choices may be restrictive for highly nonstationary signals or more complex outputs, and that sparse or noisy functional observations may require further smoothing or regularization (Kanai et al., 4 Sep 2025). In hydrological calibration, the discretization-point set is critical but subjective, the algorithmic parameters $\mathcal{F}$31, $\mathcal{F}$32, $\mathcal{F}$33, and $\mathcal{F}$34 require tuning, and no explicit model structural error is included; the procedure returns the closest available approximation rather than a quantified posterior over parameters (Bhattacharjee et al., 2017).

For general history matching, the geometry of the non-implausible region remains a central difficulty. The SMC analysis emphasizes that such regions can be multi-modal and highly irregular, that uniform sampling is hard, and that richer functional information can intensify this difficulty (Drovandi et al., 2017). High-dimensional functional outputs also require aggressive reduction: event-generator calibration uses principal variable analysis and selected bins, precisely because emulating all 432 outputs jointly is not straightforward (Iskauskas et al., 25 Feb 2026). A plausible implication is that practical FHM is often defined as much by its reduction strategy as by its implausibility formula.

The “Fuzzy Hierarchical Multiplex” strand adds a different set of limitations. Its ODE-style derivation is described as heuristic, the paper does not provide formal convergence, existence, or stability proofs in a rigorous way, and the framework is characterized as producing outputs “aligned with the data and not so much the dynamics,” which indicates a predominantly data-fitting orientation (Kafantaris, 10 Dec 2025). This is not a criticism of the framework’s intended optimization role, but it marks a different epistemic stance from the uncertainty-aware screening tradition of statistical history matching.

Future directions already appear in the literature. Explicit FHM proposes broader application to climate, environmental, oceanic, and engineering models with functional outputs, as well as possible use in sensor-network design (Kanai et al., 4 Sep 2025). Level-set history matching suggests extension from scalar functionals to more explicit functional thresholds over trajectories (Baker et al., 2020). Bayes Linear HM for event generators suggests replacing raw histograms by lower-dimensional hyperparameters such as peak location, width, and tails (Iskauskas et al., 25 Feb 2026). The fuzzy multiplex framework suggests further formalization of implication logic, more diverse optimization applications, integration with neural networks, and exploration of prediction and generation beyond history matching (Kafantaris, 10 Dec 2025). Together these directions indicate that FHM is developing along two connected but distinct trajectories: one centered on uncertainty-aware screening for function-valued simulators, and another centered on hierarchical fuzzy calibration of metric-producing mappings.