1 Introduction↩︎

Estimating probability density functions (PDFs) from data is a central task across the quantitative sciences and the backbone of likelihood-based and simulation-based inference (SBI) [1], where parameters of interest are extracted by comparing observed data to a model of its expected distribution. When the PDF cannot be written in closed form, or is accessible only through an expensive forward simulation, it must be learned from samples.

Traditionally such densities are estimated with histograms or kernel density estimators, which struggle in high dimensions. Normalizing Flows (NFs) [2] have emerged as a powerful alternative: they learn an invertible mapping between a simple base distribution, such as a multivariate Gaussian, and the target density, giving both efficient sampling and an explicit, normalized likelihood in high-dimensional spaces. They have found widespread success across the physical sciences, with applications ranging from astrophysics to lattice field theory [3]–[6].

In many settings, however, a single fixed density is not enough: the object of interest is how the density changes as a function of continuous parameters. These may be physical parameters of interest, instrument calibrations, or systematic uncertainties that reshape the distribution. This dependence is rarely known analytically, so it too must be inferred from simulations produced at different parameter values. Doing so naively does not scale: capturing the joint effect of \(K\) parameters would require sampling their full \(K\)-dimensional configuration space, whose cost grows exponentially with \(K\).

This challenge is acute in high energy physics. At the LHC, parameters of interest are extracted by fitting observed data against simulated templates while profiling hundreds of nuisance parameters that encode systematic uncertainties [7], [8]. The effect of each uncertainty is propagated through a pair of \(\pm 1\sigma\) template variations and interpolated with a low-order polynomial: a robust and widely used procedure, but one tied to binned, low-dimensional summaries of the data. Unbinned fits over high-dimensional observables would extract substantially more information per event, provided the response of the full density to each systematic can be modeled accurately and propagated to the likelihood at comparable scale.

In this work we introduce Factorizable Normalizing Flows (FNFs), a scalable approach to this parametric density-morphing problem. An FNF composes a fixed, high-fidelity flow for a reference configuration with a learnable transformation that is polynomial in the parameters and factorized over them, so that the effect of each parameter is learned in isolation and the joint response is recovered by summation at inference. This directly generalizes the \(\pm 1\sigma\) template-variation procedure to a continuous, differentiable, high-dimensional setting, while requiring only the per-systematic variations that analyses already produce: the combined nuisance space is never sampled. Here we focus on validating the parametric density modeling itself: we show that the FNF reproduces the imposed deformations and closely matches the optimal likelihood, with optional interaction terms capturing the residual correlations between parameters. Its deployment in a profiled, unbinned likelihood fit to extract parameters of interest is the subject of dedicated forthcoming work.

This paper is organized as follows: in Section 2 we discuss the challenges of modeling systematic uncertainties in unbinned likelihood fits and motivate the need for a new approach. In Section 3 we introduce the Factorizable Normalizing Flow architecture and describe its key components. In Sections 4 and 5 we present experimental results demonstrating the effectiveness of our approach on a toy dataset, and in Section 6 we discuss related work in the field. Finally, we conclude with a summary of our findings and potential future directions in Section 7.

2 Systematic uncertainties in unbinned likelihood fits↩︎

We consider a typical high energy physics (HEP) measurement, where the goal is to infer parameters of interest from observed data while accounting for systematic uncertainties. Each source of uncertainty is encoded by a nuisance parameter, collected in a vector \(\nu\), and the central object of the analysis is the conditional density \(p(y \mid x, \nu)\), where \(y\) are the observed features whose distribution we model, \(x\) are conditioning variables, and \(\nu\) are the nuisance parameters. Traditionally the likelihood is built by binning the data in one or two features and comparing the observed and expected counts per bin; an unbinned fit in higher dimensions instead requires the full density \(p(y \mid x, \nu)\) itself. In either case the dependence on \(\nu\) is not known analytically and must be learned from simulation, which is the central modeling challenge addressed in this work.

In HEP, the effect of systematic uncertainties is traditionally modeled through template variations [7], [8]. The impact of each uncertainty source is captured independently: the corresponding nuisance is shifted to the edges of its \(68\%\) confidence interval (\(\nu_k = \pm 1\), the “\(\pm 1\sigma\)” variations), the simulation is re-run, and the resulting “up” and “down” histograms bracket the nominal one. The per-bin yield is then interpolated as a smooth polynomial in each nuisance, of at least quadratic order [9], so that the nuisances can be profiled continuously during the likelihood fit [10]. This construction is the workhorse of LHC analyses: it is robust, and because each source is treated separately it scales to hundreds of nuisance parameters. Its limitations are equally structural: it operates on binned, low-dimensional summaries of the data, and it assumes that the effects of different sources factorize, neglecting their correlations, an approximation that simplifies the fit but can bias it.

Carrying this idea to unbinned, high-dimensional fits places several demands on the modeling method:

• Fidelity: it must accurately capture the effect of each systematic on the input density, so that the learned deformation reflects the true response to each source.

• Profiling: it must expose a smooth, differentiable dependence on \(\nu\), so that the nuisances can be profiled during the fit.

• Scalability: it must scale gracefully with the number of nuisances, both in the simulation required to train it (no sampling of the joint nuisance space) and in the cost of evaluating the density and its Jacobian during the fit.

• Interpretability: the learned transformation should remain interpretable, exposing how each nuisance reshapes the distribution.

A direct way to bring this picture to generative models is to model the effect of systematic uncertainties as parametric deformations of a learned density. Two routes are available:

1. make the generative model directly conditional on the nuisances, learning \(p(y \mid x, \nu)\) end-to-end from samples generated at different values of \(\nu\);

2. keep a fixed nominal density and compose it with an additional transformation that captures the deformation induced by the nuisances, \(p(y \mid x, \nu) = p_{\text{nom}}(y \mid x) \circ T(y \mid x, \nu)\).

The first route is conceptually simple but does not scale: to capture the joint effect of \(K\) nuisances the model must be trained on samples spanning their full \(K\)-dimensional space, whose size grows exponentially with \(K\). The second route is more promising, as it isolates the deformation in a dedicated transformation \(T\) and, in the spirit of the template-variation approach, opens the door to factorizing it over the individual nuisances. This is the route we take. Realizing it, however, requires a transformation that respects this factorization, stays polynomial in \(\nu\) for smooth profiling, and keeps the likelihood tractable, a combination that calls for a dedicated architecture, which we introduce next.

3 Factorizable Normalizing Flows↩︎

To address the challenge of scalable density estimation under systematic uncertainties, we introduce the Factorizable Normalizing Flow (FNF). Its central design choice is to structurally decouple the modeling of the complex, multidimensional distribution from the modeling of systematic variations: the density is a composition of a fixed, high-fidelity nominal model \(p_{\text{nom}}(y \mid x)\), trained once on high-statistics nominal simulation, and a learnable systematic transformation \(T_\nu\) that captures how probability mass migrates under systematic perturbations.

3.1 The idea↩︎

The FNF rests on a simple separation of concerns. Rather than learning, in one shot, both the shape of the data and how every systematic distorts it, we split the problem in two. First, we learn an accurate model of the data in its reference configuration, the nominal simulation at \(\nu = 0\), once and for all. Then, for each systematic, we learn a smooth deformation of the feature space, a warping that morphs the reference distribution into the one observed when that systematic is switched on (Figure 1).

To judge how likely an observed event is for a given setting of the systematics, we apply the deformation in reverse: we carry the event back to the reference frame and ask how probable it was there, correcting for how much the warp stretched or compressed the space around it. This last correction is the familiar change-of-variables (Jacobian) factor that makes the construction a proper, normalized probability density.

The key property is that each systematic contributes its own deformation and, to a good approximation, these deformations add up. We can therefore learn each one in isolation, from just the \(\pm 1\sigma\) samples already produced for that systematic, and sum their contributions to recover the joint response, with no need for simulations in which several systematics are varied together. This additivity is exact when the systematics act independently and otherwise an approximation, whose neglected correlations, and the optional interaction terms that restore them, we examine in Section 4.4. In this sense the FNF is a direct, continuous generalization of the template-variation recipe of Section 2: the discrete \(\pm 1\sigma\) histograms become a smooth, learnable warp of the full density, and the per-bin polynomial interpolation becomes a polynomial dependence of that warp on \(\nu\).

3.2 Morphing as a pullback↩︎

Mathematically, we treat the systematic-aware density as the pullback of the nominal density by the diffeomorphism \(T_\nu\). If \(y\) is an observable in the “distorted” space (for example, data affected by detector effects), the transformation \(T_\nu(y \mid x, \nu)\) maps it back to the “nominal” reference space \(y_{\text{nom}}\). The likelihood is evaluated by correcting the event and scoring it under the nominal model, penalized by the volume change: \[y_{\text{nom}} = T_\nu(y \mid x, \nu) \quad \Rightarrow \quad p(y \mid x, \nu) = p_{\text{nom}}\!\left(T_\nu(y \mid x, \nu) \mid x\right) \left|\det \nabla_y T_\nu(y \mid x, \nu)\right|.\] For this construction to serve an unbinned fit, \(T_\nu\) must meet the two requirements: its Jacobian determinant must be cheap to evaluate, and its dependence on \(\nu\) must be polynomial and factorized over the individual nuisances. We satisfy both with the construction below.

3.3 A tractable, factorized transformation↩︎

We realize \(T_\nu\) as an affine [11], autoregressive [2], [12], [13] transformation. The autoregressive structure makes the Jacobian triangular, so its determinant is simply the product of the diagonal and costs only \(O(D)\) for \(D\) feature dimensions. For each dimension \(j\), \[y_{\text{nom}, j} = y_j \cdot \exp\!\left(s_j(y_{<j}, x, \nu)\right) + t_j(y_{<j}, x, \nu),\] where the scale \(s_j\) and shift \(t_j\) depend only on the preceding features \(y_{<j}\), the conditioning \(x\), and the nuisances \(\nu\).

We interpret the systematic deformation as a Taylor expansion of the transformation around the nominal point \(\nu = 0\), and accordingly decompose the scale and shift as polynomials in \(\nu\), factorized over the individual nuisances and including optional pairwise interaction terms (Figure 2): \[\begin{align} \label{eq:factorizable95flows} s_j(y_{<j}, x, \nu) &= \sum_{k=1}^K \left( \nu_k\,\alpha_{j}^{k}(y_{<j}, x) + \nu_k^2\,\beta_{j}^{k}(y_{<j}, x) \right) + \sum_{1 \le k < \ell \le K} \nu_k\,\nu_{\ell}\,\phi_{j}^{k\ell}(y_{<j}, x) \\ t_j(y_{<j}, x, \nu) &= \sum_{k=1}^K \left( \nu_k\,\gamma_{j}^{k}(y_{<j}, x) + \nu_k^2\,\delta_{j}^{k}(y_{<j}, x) \right) + \sum_{1 \le k < \ell \le K} \nu_k\,\nu_{\ell}\,\psi_{j}^{k\ell}(y_{<j}, x) \end{align}\tag{1}\] The dependence on \(\nu\) is thus additive across nuisances, each acting through its own set of coefficients. This additivity in the affine parameters is what averts the curse of dimensionality of learning the joint conditional \(p(y \mid x, \nu)\) directly, and what lets each nuisance be trained in isolation. The coefficients \(\alpha, \beta, \gamma, \delta, \phi, \psi\) themselves carry the full kinematic complexity: they are the outputs of deep neural networks \(\Psi^{k}\) conditioned on \(y\) and \(x\). We implement \(\Psi^{k}\) as masked multi-layer perceptrons (masked MLPs) [2], [14], so a single network produces the autoregressive coefficients for all dimensions \(j\) in one forward pass while respecting the causal constraint (dependence only on \(y_{<j}\)). The result is expressive enough to capture deep feature correlations yet structured enough to scale linearly with the number of systematics.

The cross-terms \(\phi\) and \(\psi\) are optional: they capture pairwise interactions between nuisances, but can be dropped whenever the systematics are treated as independent (factorizable), recovering the standard assumption of the template approach. Finally, stacking several such layers, interleaved with permutations [15] that mix the feature dimensions, yields a highly expressive transformation that retains precise, polynomial control over the dependence on \(\nu\).

Together, these choices realize the transformation anticipated in Section 2: factorized over the nuisances, polynomial in \(\nu\), and equipped with a tractable likelihood. The FNF thereby generalizes the template-variation method to a continuous, high-dimensional setting and, because the response is additive, each systematic can be trained independently from only its \(\pm 1\sigma\) samples, with the joint effect recovered by summation. This is what makes the approach feasible even for analyses with many systematic uncertainties.

4 Experiments↩︎

We validate the FNF on a controlled toy dataset that reproduces, in a low-dimensional setting, the essential features of a systematic-aware density estimation problem: a nontrivial nominal density and interpretable nuisance deformations.

4.1 Toy dataset↩︎

Each event is described by a binary class label \(c \in \{A, B\}\) (drawn with equal probability), a two-dimensional vector \(x=(x_1,x_2)\), and a two-dimensional feature vector \(y=(y_1,y_2)\), also called the “score”. The pair \((c, x)\) provides the conditioning information, while \(y\) is the observable whose density we model. The kinematics are drawn from two class-dependent Gaussian clusters, \[x \mid c \sim \mathcal{N}\!\big(\mu_c,\;\mathrm{diag}(\sigma_c^2)\big), \qquad \begin{align} \mu_A &= (-0.5,\, 0), & \sigma_A &= (0.9,\, 0.6),\\ \mu_B &= (+0.5,\, 0), & \sigma_B &= (0.6,\, 0.4). \end{align}\] The nominal feature density \(p_{\text{nom}}(y \mid x, c)\) is a bivariate Gaussian whose mean, per-axis spread and correlation depend nontrivially on the kinematics, \[\begin{align} \mu(x) &= \big(\,\sin(1.5\,x_1) + 0.3\,x_2,\;\;0.3\,x_1^2 - 1.2 + 0.5\sin(x_2)\,\big),\\ \sigma_1(x) &= \mathrm{softplus}(0.4\,x_1 + 0.1), \qquad \sigma_2(x) = \mathrm{softplus}(-0.2\,x_1 + 0.4),\\ \rho(x) &= 0.8\,\tanh\!\big(0.5\,(x_1 + x_2)\big). \end{align}\]

Figure 3 illustrates the \(p(x)\) density and the marginal \(p(y \mid c)\).

Two nuisance parameters drive the systematic variations, each acting with opposite sign on the two classes (anti-correlated response) and simultaneously on the kinematics and on the conditional feature density:

• \(\nu_{\text{shift}}\) translates the cluster centroids along \(x_1\), \(x \to x \mp \nu_{\text{shift}}\, s_{\text{shift}}\, \hat{d}\) (with \(\hat{d}=(1,0)\); \(-\) for class \(A\) and \(+\) for class \(B\)), and adds a linear, \(x\)-dependent shift to the score mean;

• \(\nu_{\text{squeeze}}\) applies a volume-preserving (\(\det = 1\)) axis-anti-correlated squeeze of the kinematics about each centroid, \(x \to \mu_c + \mathrm{diag}(e^{+\alpha}, e^{-\alpha})\,(x - \mu_c)\) with \(\alpha = \nu_{\text{squeeze}}\, s_{\text{squeeze}}\) (the sign mirrored between classes), together with an exponential rescaling of the score spread.

Because both nuisances perturb the conditioning kinematics and the score itself, they induce a genuine, interpretable dependence \(p(y \mid x, c, \nu)\) that the FNF must capture. Figure 4 shows their effect on the kinematic density \(p(x \mid c)\) (top) and on the marginal score density \(p(y \mid c)\) (bottom), comparing the nominal distribution with the \(+1\sigma\) variation of the two nuisances.

The generator deformations are set by the kinematic scales \(s_{\text{shift}}=0.3\), \(s_{\text{squeeze}}=0.2\) and the feature-space response scales \(0.3\) and \(0.2\), so that both nuisances induce a sizeable, \(x\)-dependent deformation of the density. All quantitative results below are evaluated on a dataset of \(10^5\) samples.

4.2 Model architecture and training↩︎

Since both the kinematics \(x\) and the score \(y\) are affected by the nuisances, we factorize the joint density as \(p(x, y \mid c, \nu) = p(x \mid c, \nu)\, p(y \mid x, c, \nu)\) and model each factor with a dedicated FNF. This split is a modeling choice rather than a necessity: a single FNF could model the full four-dimensional joint density \(p(x, y \mid c, \nu)\) directly. We adopt it because it mirrors the common structure of a HEP analysis, where the conditioning kinematics \(x\) are kept separate from the discriminating score \(y\), so that each component can be modeled, trained, and validated on its own. Each factor is in turn the composition of a fixed nominal base flow, trained once on the nominal (\(\nu = 0\)) simulation, and a learnable factorizable residual transformation \(T_\nu\) that pulls the observed sample back to the nominal reference, exactly as in Section 3.

4.2.0.1 Nominal base flows.

The nominal densities \(p_{\text{nom}}(x \mid c)\) and \(p_{\text{nom}}(y \mid x, c)\) are modeled with autoregressive Neural Spline Flows (NSF) [16], using monotonic rational-quadratic spline transforms with \(20\) bins. The kinematic flow uses \(2\) transforms and is conditioned on the one-hot class label \(c\); the score flow uses \(3\) transforms and is conditioned on \((c, x)\). Both use masked MLPs with hidden layers \((128, 128, 128)\). They are trained by maximum likelihood (minimizing the negative log-likelihood) on nominal samples drawn on the fly from the generator at \(\nu = 0\), and are then frozen.

4.2.0.2 Factorizable residual.

The residual transformation \(T_\nu\) is implemented as described in Section 3. Each internal network’s output layer is initialized to zero so that \(T_\nu\) starts as the identity. For this experiment we use two independent nuisances and disable the pairwise cross-terms. We stack \(2\) residual layers for the score flow, interleaved with fixed reverse permutations to mix the two feature dimensions, and \(1\) layer for the kinematic flow; the per-nuisance masked MLPs use hidden layers \((128,128,128)\) for the score and \((128,128)\) for the kinematics. The architecture and hyperparameters are summarized in Table 1. All models are implemented in PyTorch [17] on top of the zuko normalizing-flow library [18], and our implementation, together with the toy dataset, is publicly available [19].

4.2.0.3 Training the residual.

With the base flows frozen, only the residual layers are trained, again by maximum likelihood. Crucially, the residual is trained only on the per-nuisance \(\pm 1\sigma\) templates: for each nuisance we generate samples with that nuisance set to \(\nu_k = \pm 1\) and all others at zero, yielding four template generators for the two nuisances (\(\nu_{\text{shift}} = \pm 1\), \(\nu_{\text{squeeze}} = \pm 1\)). Each training batch mixes equal-size chunks from the four templates, labeled by their nuisance vector \(\nu\), and the model maximizes their likelihood under \(p(\cdot \mid x, c, \nu) = p_{\text{nom}} \circ T_\nu\). Because each nuisance is only ever seen in isolation, the factorized additive structure is what allows the joint, multi-nuisance response to be recovered at inference by summation, never requiring training samples from the combined \(K\)-dimensional variation space. We assess the validity of this factorization assumption, and the effect of enabling the pairwise cross-terms, in Section 4.4.

4.2.0.4 Optimization.

All flows are optimized with AdamW [20] at a learning rate of \(10^{-4}\) and batch size \(1024\), using a linear warm-up over the first epoch followed by cosine decay. The nominal base flows are trained for \(100\) epochs and the residual transformations for \(60\) epochs, each epoch comprising \(1000\) steps.

4.3 Results↩︎

We first assess the FNF applied to the kinematic density \(p(x \mid c, \nu)\). Figure 5 (a) shows the learned affine response of the residual transformation, the scale factor \(\exp(s)\) and the shift \(t\), evaluated at several representative points in \(x\) for both nuisances. For the \(\nu_{\text{shift}}\) deformation the scale remains at unity (i.e.the log-scale \(s \approx 0\)), while the shift grows linearly with \(\nu_{\text{shift}}\), exactly as expected for a pure translation. For the \(\nu_{\text{squeeze}}\) deformation, instead, both the scale and the shift are nontrivial and depend on the specific point in \(x\), as expected for the volume-preserving squeeze. Figure 5 (b) shows the corresponding displacement field \(\Delta x = T_\nu(x \mid c, \nu) - x\) encoded by the residual: the transformation maps the deformed distribution back onto the nominal reference (dashed gray contours), demonstrating that the model learns a smooth deformation that recovers the nominal density.

We then validate the residual on the conditional score density \(p(y \mid x, c, \nu)\). This is the more stringent test: unlike the kinematics, the score is correlated and its nominal shape varies with the conditioning kinematics, so the residual must reproduce an \(x\)-dependent deformation. Figure 6 (a) shows the learned affine response of \(y_1\) at the conditioning point \(x=(1,1)\), evaluated at four representative points in \(y\) and for both classes.

The two nuisances leave qualitatively distinct signatures, and the residual assigns each to the correct affine component. Under \(\nu_{\text{shift}}\) the scale stays close to unity while the shift grows linearly, reproducing the pure translation of the score mean that defines this nuisance. Under \(\nu_{\text{squeeze}}\) the response is instead carried by the scale, whose logarithm grows linearly with the nuisance (an exponential rescaling of the score spread), accompanied by a smaller, \(y\)-dependent shift. The two classes respond with opposite sign throughout (solid versus dashed), recovering the anti-correlated class dependence built into the generator, while the spread of the curves across the four \(y\) points shows that the learned response is local rather than a single global affine, as the autoregressive networks permit.

Figure 6 (b) confirms this at the level of the full distribution. At \(+1\sigma\) the displacement field \(\Delta y = T_\nu(y \mid x, \nu) - y\) is a near-uniform translation for \(\nu_{\text{shift}}\) and a convergent or divergent squeeze for \(\nu_{\text{squeeze}}\), with opposite orientation for the two classes. In every case the deformed model density (solid) overlays the generator truth (black), and the residual transports it back onto the nominal reference (dashed gray), demonstrating closure not only in the learned affine parameters but in the morphed score distribution itself.

4.4 Additivity and correlated effects↩︎

The results above validate each nuisance in isolation. We now turn to the combined response, with both nuisances active simultaneously. By construction, the factorized FNF reconstructs the joint deformation as the sum of the per-nuisance scale and shift contributions, so that the scale factors compose multiplicatively across nuisances. Each nuisance’s own response is retained up to second order through its linear and quadratic terms; the additive structure across nuisances, with the cross-terms omitted, encodes the assumption that distinct systematic sources act independently on the density. It is exact whenever that assumption holds and deviates only through the neglected pairwise interactions. This is precisely the factorization underlying the standard template-variation approach, where each uncertainty is described by its own independent variations and correlated shape effects between sources are not modeled at the template level.

To quantify the impact of this approximation, we compare two configurations: (i) the factorized model, with cross-terms disabled and each nuisance trained on its own \(\pm 1\sigma\) samples; and (ii) the correlated model, with the pairwise cross-terms \(\phi\) and \(\psi\) enabled and trained on additional samples in which both nuisances are varied simultaneously. For both we quantify the goodness of fit across the full nuisance plane through the per-event excess negative log-likelihood with respect to the optimal (truth) likelihood. Averaged over events generated by the truth at a given \(\nu\), this excess is the Kullback-Leibler (KL) divergence \[\label{eq:kl} \Delta\mathrm{NLL}(\nu) \equiv \mathrm{KL}\!\left[\,p_{\text{truth}}\,\big\|\,p_{\text{model}}\,\right] = \mathbb{E}_{p_{\text{truth}}}\!\left[\,\log p_{\text{truth}}(x, y \mid c, \nu) - \log p_{\text{model}}(x, y \mid c, \nu)\,\right] \ge 0,\tag{2}\] where the expectation runs over events \((c, x, y)\) drawn from the truth at \(\nu\). It vanishes if and only if the model matches the truth, giving an absolute, per-event measure of fit quality in nats.

Figure 7 compares the two configurations across the nuisance plane. The factorized model reproduces the optimal likelihood to within a few times \(10^{-3}\) nats per event near the training axes and around the nominal point. Away from the axes, however, the missing interaction grows steadily, exceeding \(10^{-2}\) in the corners where both nuisances are large and act simultaneously. Enabling the cross-terms, constrained by the additional combined-variation samples, suppresses this residual across the whole plane, which stays below \(3.3\times10^{-3}\) nats per event, an improvement of more than a factor of four in the most correlated region. The factorized approximation is thus adequate when the nuisances are moderate or well constrained, whereas the cross-terms become necessary to preserve fidelity under large, simultaneous variations.

Figure 8 makes this interaction explicit. Decomposing the correlated model’s affine response into its additive part and the bilinear remainder, the cross-term is the small correction that the factorized model omits: it vanishes on the axes, where a single nuisance is active, follows the \(\nu_{\text{shift}}\,\nu_{\text{squeeze}}\) product across the plane, and is largest in the corners. Its magnitude stays one to two orders of magnitude below the full response, which is why the factorized approximation holds for moderate or single-nuisance variations and breaks down only where both nuisances are simultaneously large, exactly the region where the goodness of fit in Figure 7 degrades.

These two configurations span the available precision-versus-cost trade-off. The factorized model is the cheapest: it requires only the per-nuisance variations and scales linearly with the number of uncertainties, exactly as in the established template approach. When residual correlations between sources are non-negligible, the FNF can account for them natively by enabling the cross-terms, at the cost of generating combined-variation samples and a number of interaction terms that grows with the number of nuisance pairs. The same architecture therefore interpolates between the standard factorized approximation and a fully correlated model, with the level of approximation chosen explicitly.

5 Discussion↩︎

This formulation has several key features:

• Factorized, modular structure: Each nuisance \(\nu_k\) is handled by a separate network \(\Psi^{k}\), so its effect is modeled independently (when excluding cross-terms). This lets us train each variation from its own \(\pm 1\sigma\) samples alone, avoids the combinatorial explosion of sampling the full \(K\)-dimensional space, and allows new uncertainties to be added without retraining the existing model.

• Interpretability: The learned coefficients carry a direct meaning: \(\alpha^{k}\) and \(\gamma^{k}\) give the linear response of the feature distribution to \(\nu_k\), while \(\beta^{k}\) and \(\delta^{k}\) capture the quadratic effects.

• Efficiency: The joint transformation is evaluated in a single forward pass by summing the per-nuisance outputs, and its triangular Jacobian makes the density fast to evaluate during the fit.

• Expressivity: Despite the factorized structure, the transformation remains highly expressive: the neural-network coefficients and the autoregressive structure capture complex, non-linear dependencies on the features \(y\) and kinematics \(x\).

5.0.0.1 Interplay with Optimal Transport

The morphing transformation encoded by the Factorizable Normalizing Flow is not guaranteed to be “optimal” in the Optimal Transport (OT) sense, i.e., the unique map that transports the probability mass between the two distributions at minimal cost [21]. This is not a limitation of our approach. On the contrary, by letting the solution live in a larger space, without imposing uniqueness or minimal-cost constraints, the FNF can be trained directly by maximum likelihood, whereas constructing an explicit OT map, for example through convex potentials or continuous flows, is typically far more computationally expensive [22]–[25]. Optimal transport has nonetheless found fruitful application in high energy physics, for instance to define a metric on the space of collider events [26] and to map between distributions for calibration and inference [27], [28].

6 Related work↩︎

Encoding systematic uncertainties in unbinned likelihood fits is a crucial topic that has been addressed from several directions. A broad family of methods falls under simulation-based inference [1], [29], where the intractable likelihood, or the likelihood ratio, is learned directly from simulated samples rather than modeled analytically. In the effective field theory (EFT) context, the known parametric dependence of the matrix element can be exploited to learn per-event likelihood ratios and scores [30]. Closely related to our approach are methods that parameterize the density or the likelihood ratio with generative models trained on simulation: Normalizing Flows have been used to unify simulation and inference [31], generative models enable data-driven high-dimensional inference [32], and machine-learned systematic uncertainties have been incorporated into unbinned cross-section measurements [33] and global SMEFT analyses [34].

The ATLAS Collaboration has recently deployed a neural SBI technique for parameter estimation [35], in which neural-network classifiers are trained to discriminate between samples generated under different parameter hypotheses. The resulting parameterized classifier yields the per-event likelihood ratio as a continuous function of the parameters of interest and the nuisance parameters, enabling fully unbinned inference on high-dimensional observables; the technique was applied to the measurement of off-shell Higgs boson production in the \(H \to ZZ \to 4\ell\) channel [36]. A complementary, tree-based strategy is the refinable modeling of Ref. [37], where boosted information trees learn the parametric dependence of the likelihood ratio as a refinement of a base model for unbinned SMEFT analyses, replacing the neural-network estimator with a gradient-boosted tree ensemble.

These approaches share the same underlying idea as our work: the effect of the parameters is encoded as a parametric deformation of the model, learned from simulation and interpolated smoothly across the parameter space. The key distinction is what is deformed and how it is learned. Both the ATLAS and the tree-based methods parameterize the likelihood ratio and learn it discriminatively with classifiers, neural networks in one case and boosted decision trees in the other. In contrast, the FNF parameterizes the density itself through an invertible normalizing-flow transformation, providing an explicit and normalized model of \(p(y \mid x, \nu)\) whose deformation is interpretable and whose Jacobian makes the likelihood directly tractable.

6.1 Future directions↩︎

Factorizable Normalizing Flows are expected to find many direct applications in future unbinned likelihood-fit workflows. They are particularly well suited to the methodology of HEP analyses, but they are a general density-estimation tool that could also prove useful in other fields that already rely heavily on Normalizing Flows for inference, such as cosmology [38] and gravitational-wave astrophysics [39]. Within HEP, we foresee direct applications to detector calibration [40], unbinned parameter estimation, and unfolding [41]–[45]. In particular, the use of the FNF likelihood to profile nuisance parameters and extract parameters of interest in a full unbinned fit is the subject of a dedicated forthcoming publication.

7 Conclusions↩︎

We have introduced the Factorizable Normalizing Flow (FNF), a generative framework for modeling how a probability density deforms under continuous parameters such as systematic uncertainties. By decoupling a fixed, high-fidelity nominal flow from a learnable transformation that is polynomial in the parameters and factorized over them, the FNF generalizes the classical template-variation approach to a continuous, high-dimensional, and fully differentiable setting. Its central practical advantage is that each variation is learned independently from only its own \(\pm 1\sigma\) samples, avoiding the combinatorial cost of sampling the joint parameter space, while the joint response is recovered by summing the per-parameter contributions at inference.

On a controlled problem with two interpretable nuisances acting jointly on the kinematics and on the conditional score, the learned scale and shift responses reproduce the imposed deformations and transport the distorted distributions back onto the nominal reference. Quantitatively, the factorized model attains closure to the optimal likelihood across the trained region, and enabling the optional cross-terms recovers the residual correlations where parameters vary strongly together, so that a single architecture spans the full range from the standard factorized approximation to a fully correlated model, with the level of approximation chosen explicitly. The resulting model is interpretable, scales linearly with the number of parameters, and is efficient to evaluate, with a tractable likelihood suited to profiling. By making systematic-aware density estimation practical at scale, the FNF lowers a key barrier to the adoption of unbinned likelihood fits in high energy physics, and offers a general, reusable tool for any inference workflow in which probability densities must be morphed as a function of continuous parameters. Its deployment in a profiled, unbinned likelihood fit to extract parameters of interest is the subject of dedicated forthcoming work.

References↩︎