Conditional Flow Models Convert Epistemic Uncertainty into Apparent Aleatoric Variance

Introduction

The role of generative models in machine learning has expanded rapidly. Rather than returning a single point prediction, these models produce a conditional distribution over possible outputs. Perhaps surprisingly, they have gained traction even in settings where predicting a single value, rather than distribution over values, suffices. These setting, which we will call point-prediction problems, encompass predicting high-performing actions in robotics ([8]), accurate estimates of a policy value function in reinforcement learning ([1]), and numerous other scientific estimation and time-series problems ([17]).

In this work, we focus on conditional flow-matching models (CFMs), which produce samples by integrating a learned velocity field from base noise to the output space, and suggest that these models confer a complementary benefit: they can naturally express estimation uncertainty, even in point-prediction problem in which observations are nearly noiseless.

In the idealized setting, the optimal flow $v^{\star}$ transports the base distribution to the true conditional law ([21]; [3]). In practice, however, the flow field is parametrized by a neural network $v_{\theta}$ that is learned imperfectly, due to a combination of limited training data, network capacity, and incomplete optimization (e.g. early stopping, popular in the robot learning community). Thus the learned conditional distribution is shaped not only by the target distribution, but also by the inductive biases caused by fitting the ground-truth vector-field learning with a neural network. In particular, neural networks tend to oversmooth high-frequency features during learning [30, 41, 37], causing the learned vector-field $v_\theta$ to underfit the ground-truth $v^*$.

We demonstrate, both across domains and in natural theoretical setups, that under this common underfitting behavior, conditional flow models convert epistemic uncertainty into apparent aleatoric variance. That is, the approximation error or uncertainty about the dependence of the generated variable on the conditioning variable is expressed as a larger spread in the learned conditional distribution.

Beyond scientific interest, this has actionable consequences. One can train conditional flow-matching models to fit potentially deterministic target functions and, as long as the model does not overfit, use the spread of their resulting samples to define confidence sets containing the ground-truth function, with no conformal correction.

Across the application domains we study, including RL policy evaluation, cell perturbation prediction, and robot policy learning, predictive sets built this way are approximately calibrated empirically in the underfitting regime, a property the variance-inflation phenomenon enables. A practitioner receives an uncertainty estimate “for free” alongside the point prediction.

More concretely, we make the following contributions:

• We identify a surprisingly general relationship between learning error in the conditional flow field and the resulting generated samples: whenever a learned CFM underfits the ground-truth velocity field as a function of the conditioning variable, its generated samples have conditional variance strictly larger than that of the true distribution, in expectation over the conditioning variable. This stands in sharp contrast to deterministic regression, where underfitting manifests only as oversmoothing of the predictive mean.

• We trace this variance inflation to a single mechanism: smoothing bias in vector-field learning. Any inductive bias that smooths the learned conditional vector field along the conditioning variable, whether early stopping, or explicit regularization, produces the same effect. We demonstrate this empirically on synthetic regression problems, controlled smoothing estimators, and high-dimensional offline RL benchmarks.

• We propose an idealized setting of regularized smoothing to theoretically study the variance inflation phenomenon. Under an equivalent-kernel approximation to regularized vector-field learning, the learned conditional distribution is a local mixture of true conditionals around the query point. This yields a clean bias–variance relation: on average, the learned conditional variance is at least the true conditional variance plus the squared bias of the learned conditional mean. This simplified setting serves as a useful explanatory model to understand the variance inflation phenomenon empirically observed in neural CFMs. We also discuss the validity of this approximation to neural CFMs and the regimes where the approximation fails.

• Motivated by our theoretical study, we identify another empirical consequence of the variance inflation phenomenon – the distribution learned by CFMs is approximately calibrated in the underfitting regime. We give a practical method for leveraging variance inflation for predictive coverage. Across three application domains, policy evaluation, cell perturbation prediction, and behavior cloning for robotics, predictive sets built directly from CFM samples are approximately

calibrated in the underfitting regime, and are competitive with or tighter than conformal-prediction baselines.

Related Work

Flow matching and conditional transport. Flow matching and related stochastic-interpolant and rectified-flow methods learn continuous-time transports by regressing a time-dependent vector field along a path between base noise and data [21, 22, 4]. Conditional flow matching indexes this transport by side information, yielding a model of the conditional law of an output given a conditioning variable. Conditional Flow Models (CFMs) have been applied broadly in robotics [9, 31, 44] and scientific modeling [17, 23].

Inductive biases and generalization in diffusion and flow models. Several recent works study why diffusion or score-based models generalize rather than simply memorizing the training data. Recent work has attributed this behavior to several sources of bias: variance in the denoising score-matching target [38], architectural bias [16], and selective underfitting of the learned score [35]. These papers are closest in spirit to ours, however, they focus on unconditional generation, whereas we study conditional generation and ask how underfitting with respect to the conditioning variable changes the conditional law.

Implicit regularization, spectral bias, and smoothing. Our explanation is based on smoothing bias in the learned conditional vector field. Neural networks are known to preferentially fit low-frequency or smoother structure early in training, a phenomenon often called spectral bias [30, 41]. In wide-network regimes, training dynamics can be approximated by kernel gradient flow, connecting early stopping and finite-time training to classical regularization in function space [15]. Related ideas appear throughout nonparametric regression, smoothing splines, kernel ridge regression, and inverse problems, where regularized estimators can be described by equivalent kernels or linear smoothing operators. We build on this viewpoint, but study a different object: not the bias of a point estimator, but the distribution obtained after integrating a smoothed conditional vector field.

Preliminaries

We study conditional generative modeling. Given samples from a joint distribution $p^{\star}$ over random variables $(X,C) \in\mathbb{R}^{d_x} \times\mathbb{R}^{d_c}$, where $X$ is the generated variable and $C$ is the conditioning variable, conditional generative models learn to generate new samples from the conditional distribution of $X$ given $C = c$, denoted $p^{\star}(\cdot\mid c)$. Flow models represent such distributions as continuous normalizing flows (CNF) and ODEs transporting a simple initial distribution to a terminal distribution approximating $p^{\star}$. Here, we consider the conditional variant.

Flow-matching [21, 3, 22] provides a generic recipe for learning (C)CNFs from data. Let $p_0$ and differentiable coefficients $\alpha_t,\beta_t : [0,1] \to\mathbb{R}$ satisfying $\alpha_0 = \beta_1 = 1$ and $\beta_0 = \alpha_1 = 0$ both be given. The induced (linear) stochastic interpolant is the random variable $X_t = \alpha_t X_0 + \beta_t X_1$, where $X_0 \sim p_0(\cdot\mid c)$ and $X_1 \sim p^{\star}(\cdot\mid c)$. Flow-matching is based on the following observation:

Motivated by Proposition 3.1, flow matching trains a velocity field $v_{\theta} = (v_{\theta;t}(x,c))$ by (approximately) minimizing an empirical approximation to the flow matching loss $\mathcal{L}_{\mathrm{FM}}$.

Variance Inflation in Underfit Conditional Flows

In this section, we uncover two surprising phenomena about underfit conditional flow models: First, the resulting samples have larger conditional variance than the ground-truth distribution (in expectation over the conditioning variable), even in nearly deterministic regression problems. Second, this variance inflation is a consequence of smoothing bias in the learned vector field: any mechanism that smooths the conditional vector field along the conditioning variable produces the same effect. We demonstrate the phenomenon on a synthetic 1D problem (Figure 1), establish that it is driven by smoothing (Figure 2), and verify that it transfers to a high-dimensional benchmark used in offline RL (Figure 3).

Image not extracted

Image not extracted

Image not extracted

Figure 1, given in the introduction, illustrates a one-dimensional synthetic regression problem in which the target function contains both low and high-frequency structure. When a neural-network regressor is trained with early stopping, underfitting leads to oversmoothing of the predictive mean; the learned function captures the coarse trend but fails to recover the finer oscillations.

For a conditional flow model trained under the same regime, the learned conditional mean is similarly oversmoothed, but unlike the deterministic regressor, the flow also learns substantial additional variance, visible in the wide 95% band around the mean. Even in a nearly deterministic problem, underfitting in the conditional vector field appears not only as bias in the learned mean, but also as extra variance in the generated conditional distribution. Panel (d) shows the full learned conditional density $\hat{p}(x \mid c)$, which is broad and roughly bimodal in regions where $f^{\star}$ (the conditional mean) has high-frequency oscillations within a smoothing neighborhood; the flow’s conditional distribution effectively averages over the local oscillations of $f^{\star}$. This finding is not specific to this example, and Appendix C shows that similar phenomena arise across a number of synthetic target functions.

Why do we call this uncertainty epistemic? The excess variance arises from unresolved dependence of the target on the conditioning variable: the learned model smooths over distinctions between nearby conditions that the true conditional mean resolves. This uncertainty is therefore reducible and model-dependent, rather than irreducible randomness in the data-generating process. As the model resolves more of this structure, for example with additional data or training, both the smoothing bias and the excess variance decrease. However, in a CFM this unresolved variation appears as spread within the learned conditional distribution, and hence as apparent aleatoric variance. From the learned distribution alone, this model-induced spread is not distinguishable from genuine aleatoric variation.

Variance inflation tracks smoothing bias

Figure 2 reveals that the phenomenon in Figure 1 is not tied to a single estimator or training recipe. We compare three mechanisms that induce smoothing in the learned conditional vector field: early stopping of a neural network which biases training towards low-frequency components of the target [42, 30, 37], an explicit roughness penalty $\lambda\left(\partial_c v\right)^2$ in the conditioning variable as used in classic non-parametric regression [39, 13], and kernel ridge regression with an RKHS-norm regularization penalty $\lambda$ [6]. All these methods induce smoothing of the conditional mean function, as can be observed from the plotted learned mean at different regularization strengths (Figure 2). In all three cases, this smoothing is accompanied by another consistent effect: the learned conditional variance increases as the learned mean becomes smoother. The average learned conditional variance (reported in each panel title) decreases monotonically as smoothing decreases.

Our key finding, corroborated across additional functions in Appendix C, is that the key object is the smoothing bias in vector-field learning, not flow matching alone. Flow matching simply provides a mechanism through which this smoothing bias is expressed as conditional spread. Also note in Figure 2 that the induced variance is often largest precisely when the learned mean is most biased.

Beyond one-dimensional synthetic regression: Monte-Carlo Policy Evaluation

We next test whether the variance-inflation phenomenon persists when the conditioning variable is high-dimensional. For this, we use a Policy Evaluation task, where we train a CFM to predict $Q$-values conditioned on state-action pairs $c = (s, a)$, using Monte-Carlo targets. We evaluate on two OGBench Park et al., [25] environments, Cube-Single and Antmaze-large, where the conditioning variable is high-dimensional while the target remains scalar.

Figure 3 visualizes the learned and empirical distributions against the first two principal components of the conditioning variable. At early training steps, the learned conditional distribution is noticeably more diffuse than the empirical target, but still ’covers’ the target distribution samples; at later checkpoints, the learned samples contract and align more closely with the target scatter. This mirrors the synthetic experiments: underfitting in the conditional vector field leads to excess spread, which decreases as training proceeds.

An Idealized Model for Variance Inflation

Image not extracted

Due to finite data, capacity constraints on the learned network, and optimization choices (e.g. early stopping), $v_\theta$ is an imperfect approximation of $v^*$. Because the precise form of this approximation is difficult to characterize for neural function approximators, we study the following idealized population-level regularized flow-matching problem.

Above, operator $L_c$ is a linear operator inducing smoothing in the conditioning variable. Specifically, for each fixed $(t, x)$, $L_c$ acts on the function $c \mapsto v_t(c, x)$ while treating $t$ and $x$ as fixed. This matches the regime of interest: the learned vector field may depend flexibly on the flow state $x$ and on time $t$, but its dependence on the conditioning variable $c$ is smoothed. Concrete examples include first-derivative penalties $L_c = \nabla_c$, higher-order smoothness penalties $L_c = \nabla_c^k$, and RKHS-type penalties on $c \mapsto v_t(c, x)$ Wahba, [39]. It is well-known that such regularized regression problems have equivalent kernel representations Schölkopf et al., [33]. This objective serves as an idealized model of neural-network smoothing bias; related regularized models have been studied for conditional diffusion models (Frommer et al., 2025). Appendix § 8 provides further motivation, and Appendix § 8.§ 4 numerically evaluates our equivalent-kernel approximation (Equations (5) and (5.4)) for neural CFMs on synthetic tasks.

Since $L_c$ only induces smoothing in the conditioning variable $c$ and treats $t$ and $x$ as fixed, for any fixed $(t, x)$, up to a multiplicative constant independent of $g$, the regularized least-squares problem for the slice $g(c) := v_t(c,x)$ is

For Tikhonov-type penalties, the solution is a linear smoother in the regression target, with some equivalent kernel $K_{\lambda,t,x}$ depending on $(t,x)$ [34]:

For the constant-preserving smoothers considered here, the equivalent kernel is normalized so that $\int K_{\lambda,t,x}(c,c')p_t^*(x\mid c')p_C(c')\,\mathrm{d}c'=1$. The factor $p_t^*(x\mid c')$ is the design density of the slice problem, not a separate modeling choice. The starting point of our analysis is a structural assumption on the kernel $K_{\lambda,t,x}$.

Appendix § 8 motivates this assumption and discusses its regime of validity (small $t$, large path variance, or locally slowly varying conditional mean) in detail. Under Assumption 5.1, together with mild moment and regularity conditions (Assumptions A.1, A.2, and A.3 in Appendix § 8), the conditional distribution generated by $v_{\lambda}$ satisfies the following bias–variance relation.

The proof proceeds by showing that Assumption 5.1 reduces the regularized flow-matching solution to a density-weighted local average of the true velocity field, which in turn implies that the endpoint distribution generated by the learned velocity field is a local mixture of the true conditional distributions. Taking moments of this mixture and using standard assumptions on the smoothing kernel (Assumptions A.2, A.3) gives Theorem 1. The full proof is deferred to Appendix § 8.

The comparison above is in the Loewner order. Equivalently, in every direction $u$, the average learned conditional variance of $u^\top\hat{X}_1$ is at least the corresponding true conditional variance plus the average squared directional bias of the learned conditional mean. In the scalar case, Theorem 1 says that, on average over the conditioning variable, the learned conditional variance is at least the true conditional variance plus the squared bias of the learned conditional mean. Thus, any smoothing bias incurred in fitting the conditional mean is reflected as additional spread in the generated samples. Figure 5 plots squared bias with excess variance for vector fields trained using different smoothing estimators on multiple one-dimensional synthetic distributions, and the figure reflects the relation predicted in Eq. (7).

Image not extracted

Image not extracted

Image not extracted

Coverage Properties of Conditional Flow Models

We now evaluate whether the variance induced by conditional flow models yields meaningful predictive coverage on real-world tasks. The heatmap in panel (d) of Figure 1 suggests that the excess variance of an underfit CFM has structure: it assigns more mass to regions of higher density under the true distribution. Motivated by this, we aim to test a practical consequence of the variance-inflation phenomenon: whether the learned conditional distributions assign appropriate mass to the true outcomes, even in regimes where the conditional mean may be biased. We consider three application domains: Monte-Carlo policy evaluation, cell perturbation prediction, and behavior cloning. We focus on the underfitting regime, the regime where the smoothing bias of neural networks holds.

Coverage and interval definitions. For a CFM $\hat{p}_\theta(x \mid c)$ and nominal level $\alpha\in(0, 1)$, let $C_{\theta,\alpha}(c)$ denote a predictive set containing approximately $\alpha$ probability mass. We measure empirical coverage on $\mathcal{D} = \{(x_i, c_i)\}_{i=1}^{|\mathcal{D}|}$ as $\widehat{\operatorname{Cov}}_\alpha(\mathcal{D}) = \frac{1}{|\mathcal{D}|} \sum_{i=1}^{|\mathcal{D}|} \mathbf{1}\{x_i \in C_{\theta,\alpha}(c_i)\}$. We report coverage at $\alpha\in\{0.75, 0.90, 0.95\}$. For scalar outputs (policy evaluation), we construct $C_{\theta,\alpha}(c)$ as the central $\alpha$-mass interval of the learned conditional distribution; for high-dimensional outputs (cell-perturbation, behavior cloning), we use sample-based KNN predictive sets (Wang et al., [40]). Full constructions are given in Appendix § 9.

The question investigated through these experiments is: Do predictive sets constructed from CFMs as above provide meaningful coverage? We experiment on the following tasks:

i) Monte-Carlo Policy Evaluation. We use three OGBench environments (Park et al., [25]). The CFM models the distribution of scalar Monte-Carlo returns $x$ conditioned on the state-action pair $c = (s, a)$. Figure 6 reports results for Antmaze-large; Cube-Single and Humanoidmaze-medium are deferred to Appendix § 9.§ 4.

Image not extracted

ii) Single-cell perturbation prediction. We evaluate on Norman, PBMC cytokines, and SciPlex3 from scPerturb (Peidli et al., [27]) datasets, using the CellFlow architecture (Klein et al., [17]). The CFM models the distribution of perturbed-cell states $x$, with control cells as the source distribution and a pretrained perturbation embedding as the condition $c$. Figure 7 reports Norman results; PBMC and SciPlex3 are deferred to Appendix § 9.§ 4.

Image not extracted

iii) Behavior cloning. We evaluate a flow-matching policy on PushT (Chi et al., [8]). The CFM models the distribution of short action chunks $x$ conditioned on the robot’s observation history $c$.

Figure 6–Figure 8 and 12–15 report our results. We plot CFM training and validation losses, empirical coverage, predictive-set sizes, and utility metrics. We include the loss curves to identify when the trained model begins to overfit. In particular, they indicate the training regime over which we expect coverage to remain reliable: from the beginning of training until the train and validation loss curves diverge (or validation improvement saturates). This is shaded in the figures. In practice, models are generally not trained far into the overfitting regime; techniques such as early stopping, regularization, and large-data training are used to prevent this. Here, however, we deliberately train for longer to demonstrate that coverage remains reliable in the underfitting regime and changes thereafter. For scalar-output tasks we additionally plot pointwise predictive variance against squared prediction error. Set sizes are reported as interval widths for scalar outputs and KNN radii for high-dimensional outputs, while validation log-likelihood (or task success rate for behavior cloning) is tracked to assert that high coverage is not obtained from a degenerate predictive distribution.

Image not extracted

Baselines. We compare against six baselines. For scalar outputs, det + SC applies split conformal prediction to a deterministic regressor with the same conditioning backbone as the CFM, while det + CQR uses conformalized quantile regression to obtain locally adaptive intervals. Conf CFM conformalizes the raw CFM predictive sets, whereas CFM mean + SC applies split conformal prediction to the CFM conditional mean. We additionally compare against Deep Evidential Regression (DER) [5] for scalar outputs and Heteroscedastic Gaussian Regression (HGR), using a diagonal Gaussian model for high-dimensional outputs. DER and HGR are reported for fully trained models. Complete experimental details are given in Appendix § 9. The takeaways are as follows:

1. Coverage in the underfitting regime is approximately calibrated. There are two takeaways here: i) Up till the point CFM train and val loss curves diverge (underfitting regime), coverage of CFM on val data stays approx. $\geq$ nominal level ($\operatorname{Cov}_{\alpha}(\mathcal{D}) \gtrsim\alpha$), and dips afterwards. ii) Interestingly,

coverage of CFMs on the training data stays $\geq$ nominal throughout training. This is because the training data is always ’in-distribution’ for the model. Once the model starts to overfit, the training data stays in-distribution w.r.t. to the model’s learned distribution, whereas the validation does not.

2. Variance tracks pointwise prediction error. Examples with larger flow variance also tend to have larger squared prediction error for the flow mean, with positive Spearman correlation across all three environments with one-dimensional outputs (panel (e) in Figure 6, Figure 12, Figure 13). The variance is larger precisely where the model’s mean is less accurate, consistent with the interpretation that underfitting is expressed as conditional sample variance.

Image not extracted

Image not extracted

3. CFM based coverage sets are tight. Throughout the regime of interest (before overfitting sets in): i) CFM-based coverage widths / KNN-radii are tighter than split-conformals. This is because split-conformals are not locally adaptive, and also shows that CFM-based coverage sets are. ii) CFM coverage set sizes are competitive with CQR and other baselines that produce close-to nominal coverage. iii) CFM and conformalized CFM sets have very similar sizes, showing that additional conformalization changes the sets only minimally.

Conclusion and Limitations

We identify a surprising inductive bias of conditional flow models: smoothing bias in vector-field learning causes unresolved dependence on the conditioning variable to appear as excess predictive variance, converting epistemic uncertainty into apparent aleatoric variance. We observe this phenomenon across synthetic and real-world tasks, and provide an idealized model to study it. We empirically find that this excess variance has approximate coverage properties, and, at any point in training before overfitting, predictive sets constructed directly from CFM samples achieve approximately nominal coverage without explicit calibration training.

Limitations. Our theory analyzes an idealized population-level regularized smoothing model and does not establish the variance-inflation result for arbitrary neural-network training. Although the theoretical result characterizes the variance inflation phenomenon, it does not imply formal predictive-coverage guarantees. The coverage properties studied here are empirical, and developing stronger theoretical coverage guarantees is left for future work.

Theoretical Model of Variance Inflation

This section provides the formal setup and proof of Theorem 1 from § 5. The argument proceeds in three steps.

• First, an equivalent-kernel approximation (Section A.1) reduces the regularized flow-matching solution to a density-weighted local average of the true velocity field.

• Second, this representation implies that the endpoint distribution generated by the regularized flow is a local mixture of true conditionals (Proposition A.2).

• Third, taking moments of this mixture gives an explicit decomposition of the learned covariance into true conditional covariance and unresolved variation of the conditional mean. A short operator argument then shows that, on average, this unresolved variation dominates the squared smoothing bias, yielding Theorem 1.

Setup. We restate the regularized flow-matching problem from § 5: we study the population-level optimum

as a proxy for the inductive bias of a learned conditional flow. Here $\mathcal{V}$ is the chosen vector-field class, $\lambda\ge0$ controls the strength of regularization, and $L_c$ is a linear operator acting only on the conditioning variable: for each fixed $(t,x)$, $L_c$ acts on the function $c \mapsto v_t(c,x)$ while treating $t$ and $x$ as fixed. Examples include derivative penalties such as $L_c=\nabla_c$, higher-order smoothness penalties, and RKHS-type penalties. This captures the regime of interest: the learned vector field may depend flexibly on the flow state $x$ and time $t$, but its dependence on the conditioning variable $c$ is smoothed.

Justification for studying this setup: Prior work has established that neural networks exhibit a spectral or smoothness bias during training, preferentially fitting smoother and lower-frequency components before finer-scale structure ([11, 15, 30, 41, 19]). Equation (3) captures this behavior in an idealized population-level model. Closely related, Pfrommer et al. [28] study a special case with $L_c=\nabla_c^2$ as a model of smoothness bias in conditional diffusion models. More broadly, underfitting of learned velocity and score fields has been studied extensively in the diffusion and flow-matching literature. Prior work has linked smoothness and selective underfitting of these fields to generalization, interpolation, and hallucination in generative models ([43, 32, 2, 35]). Finally, in Appendix § 8.§ 4 we numerically evaluate how well the factorized equivalent-kernel representation in Equations (5) and (5.4) approximates trained neural CFMs on synthetic one-dimensional regression tasks.

For brevity, throughout this appendix we write $\hat{v}:=v_\lambda$ for the regularized population solution.

For any vector field $v$, the endpoint distribution generated at condition $c$ is the law of $\Phi_1^v(c,Z)$ with $Z\sim\rho_0$. We write its conditional mean and covariance as

When the generated variable is scalar, we write $s_v^2(c):=\Sigma_v(c)$. In the unregularized realizable population limit, $\Sigma_v(c)$ matches the true conditional covariance of $X_1\mid C=c$. The central question of this section is what happens to this endpoint covariance when the vector field is regularized or underfit in the $c$-direction. We show that the learned endpoint distribution can have substantially larger spread than the true conditional distribution, with the excess spread reflecting unresolved variation of the true conditional mean over the model’s effective smoothing neighborhood.

The above assumption simply means that $L_2$ regression is well-specified and second moments are finite. We make no parametric assumption about $\xi$ beyond its first two conditional moments.

Assumptions: Equivalent-kernel model of regularized flow matching, and regularity conditions

The regularized objective in (A.1) penalizes variation of $c \mapsto v_t(c,x)$. For each fixed $(t,x)$, this is a supervised regression problem in $c$ whose population design distribution is the conditional law of $C \mid X_t=x$, i.e., proportional to $p_t^*(x\mid c')P_C(dc')$. Classical equivalent-kernel calculations for regularized least squares show that the fitted value at $c$ is a local average of the regression target Silverman [34], which gives Eq. (5). We now simplify this equivalent-kernel representation into the form used in our analysis.

Simplifying the Kernel. We first restate Assumption 5.1 for convenience.

We now show that Assumption 5.1 implies Lemma A.1, and then motivate Assumption 5.1 in the light of prior literature.

Proof. Substituting (5.4) into the slice solution Eq. (5) and absorbing $a_{t,x}(c)$ into the normalization, we obtain

The form in (A.5) is exactly the representation in Lemma A.1, with $W_\lambda$ obtained concretely by normalizing $K'_\lambda(c,c')p_C(c')$ in $c'$. □

Motivating Assumption 5.1. The factorization (5.4) is an equivalent-kernel approximation of the kind used for smoothing splines and other roughness-penalized estimators [34, 39, 13]. The exact kernel $K_{\lambda,t,x}$ depends on $(t,x)$ only through the design density $\rho_{t,x}(c) := p_t^*(x\mid c)p_C(c)$, so replacing it by a common kernel is justified when $\rho_{t,x}$ is locally stable over the smoothing length scale of the penalty. Equivalently, using the normal-equation resolvent $R_\rho:= (M_\rho+\lambda L_c^*L_c)^{-1}$, the second-resolvent identity

shows that the exact kernels are stable whenever the perturbation $M_{\rho_{t,x}}-M_{\bar{\rho}}$ is small in the operator norm induced by the regularized problem. This is the standard Fredholm/resolvent picture behind Tikhonov-regularized least squares and kernel ridge regression [10, 6, 18, 12]. For the Gaussian-like paths we consider, the approximation is expected to be best for small $t$, large path variance, or locally slowly varying conditional mean $\mu^*$. It is weakest near $t=1$ in low-noise regions, where $p_t^*(x\mid c)$ sharply localizes around level sets of the conditional mean.

Also, in Section A.4 we numerically test how well Eq. (A.4) approximates the velocity field learned by a neural network, on a synthetic task. This is equivalent to verifying how well Eq. (5) and Eq. (5.4) together approximate a neural CFM’s learned velocity field.

We conclude this subsection with two assumptions on the resulting kernel $W_\lambda$. First, we require mild regularity conditions ensuring that the mixed path above is the path generated by $\hat{v}$.

Second, we assume reversibility and positive semidefiniteness of the kernel — two standard properties of self-adjoint smoothing operators on $L^2(P_C)$. In particular, for example, they are satisfied by the Markov kernels arising from Tikhonov regularization with symmetric quadratic penalties [12].

Formal Main Results

We are now ready to state our formal main results, deferring proofs to Section A.3. First, we show that under our assumptions, smoothing over the flows directly translates to smoothing over the conditionals. This is the key bridge between regularization of the vector field and variance inflation in the generated distribution.

At this point, the flow-dynamics part of the argument is complete: the regularized flow generates a local mixture of the true conditional distributions. The remainder of the variance analysis follows from the moments of this mixture.

Let

By Proposition A.2,

Therefore, by the laws of total expectation and total covariance,

The decomposition in (A.15) already exposes the variance-inflation mechanism. The first term is the true conditional covariance averaged over the smoothing neighborhood, while the second is the variation of nearby conditional means that the smoothed model has failed to resolve.

Variance in the Scalar Case. For simplicity, first consider scalar outputs. Write $\mu^{\star} = f^{\star}$ and $\Sigma^{\star} = \sigma^{2}$, and define

together with

Specializing (A.14)–(A.15) to the scalar case gives

Equation (A.21) makes the source of the additional variance explicit. The term $s_{\lambda}^{2}(c)$ is the variance of $f^{\star}$ over the smoothing neighborhood $W_{\lambda}(\cdot\mid c)$, i.e., the unresolved variation of the regression function over the model’s effective receptive field around $c$. The term $\bar{\sigma}_{\lambda}^{2}(c)$ is the local average of the irreducible noise variance $\sigma^{2}$ over the same neighborhood. In the unregularized limit $(\lambda\to0, W_{\lambda} \to\delta_{c})$, the first term vanishes and the second collapses to $\sigma^{2}(c)$, recovering the correct conditional variance.

Thus, the learned conditional variance exceeds the local average noise level by exactly $s_{\lambda}^{2}(c)$. The remaining question is whether this unresolved variation is large enough to account for the squared bias $\delta_{\lambda}^{2}(c)$ induced by the same smoothing operator.

Average scalar excess variance. Proposition A.3 can be equivalently phrased directly in terms of the conditional variance learned by the flow. Let

and

Combining (A.21) with Proposition A.3 and using $\mathbb{E}[\bar{\sigma}_{\lambda}^{2}(C)] = \mathbb{E}[\sigma^{2}(C)]$ by stationarity,

Equivalently,

Thus, on average over the conditioning variable, the learned conditional variance is at least the true conditional variance plus the squared smoothing bias of the learned conditional mean.

This gives a precise sense in which the conditional flow converts epistemic uncertainty into apparent aleatoric variance. The smoothing operator $S_{\lambda}$ captures unresolved variation of the regression function $f^{\star}$ in the conditioning variable, which the model cannot represent given its inductive bias. This unresolved variation manifests as bias $\delta_{\lambda}(c)$ in the learned conditional mean; (A.27) shows that the same smoothing also induces enough additional conditional spread to cover the squared bias on average. Underfit conditional flows therefore do not merely smooth their predictions toward a less informative estimate; they also produce stochastic outputs whose extra variance carries the unresolved dependence on $C$.

Vector Case. The vector case follows by applying the same scalar argument along every direction in the output space, and admits the same interpretation.

Deferred Proofs

We first prove the local-mixture representation that connects smoothing of the velocity field to smoothing of the generated conditional distribution.

Proof of Proposition A.2. For every fixed $c'$, the exact conditional path satisfies the multivariate continuity equation

Integrating this identity against $W_{\lambda}(\mathrm{d}c' \mid c)$ and exchanging the integral with the derivatives, by Assumption A.2, gives

By Lemma A.1,

The initial law is $\rho_{0}$ since $p_{0}^{\star}(\cdot\mid c') = \rho_{0}$ does not depend on $c'$. Uniqueness in Assumption A.2 then identifies $\hat{p}_{t}$ as the law generated by $\hat{v}$. Evaluating at $t = 1$ gives (A.11).

Proof of Proposition A.3 (scalar case). The proof compares two effects of the same smoothing operator: the unresolved variation $s_{\lambda}^{2}$ created by averaging nearby values of $f^{\star}$, and the squared bias $\delta_{\lambda}^{2}$ created by moving $f^{\star}$ toward $S_{\lambda}f^{\star}$. Their difference reduces to a quadratic form involving $S_{\lambda}(I - S_{\lambda})$.

By definition,

Taking expectation over $C \sim P_{C}$ and using stationarity of $W_{\lambda}$,

Similarly, since

we have

Subtracting (A.34) from (A.33) and writing inner products in $L^2(P_C)$ gives

It remains to show that this quadratic form is nonnegative. By Assumption A.3, $S_\lambda$ is self-adjoint on $L^2(P_C)$ by reversibility and positive semidefinite by assumption. As a Markov operator preserving $P_C$, it is also a contraction. Hence its spectrum lies in $[0,1]$, and so $S_\lambda(I-S_\lambda)$ is positive semidefinite. Therefore,

$\square$

Proof of Proposition A.4 (vector case). The vector result follows by applying the scalar argument to every one-dimensional projection of the conditional mean.

By (A.15),

Since $W_\lambda$ preserves $P_C$,

It therefore remains to show

Fix $u\in\mathbb{R}^d$ and define $g_u(c):=u^\top\mu^\star(c)$. Then

Applying the scalar calculation from Proposition A.3 to $g_u$ gives

where the final inequality follows from Assumption A.3. Since this holds for every $u$, the Loewner inequality follows:

Combining this with the covariance decomposition in (A.15) gives (A.29). $\square$

We emphasize that both (A.27) and (A.29) are average statements over $C$. They do not assert that the scalar excess variance dominates squared bias pointwise, nor that the learned covariance pointwise dominates the true covariance in every heteroscedastic problem. The empirical observations in § 6 go somewhat beyond this, showing that pointwise predicted variance correlates strongly with pointwise squared error of the flow’s mean.

Validity of the Equivariant-Kernel Assumption (Assumption 5.1)

In this section, we numerically test how well Eq. (A.4) approximates the velocity field learned by a neural network, on a synthetic task. This is equivalent to verifying how well Eq. (5) and Eq. (5.4) together approximate a neural CFM’s learned velocity field. Let $\hat{v}_\theta(c,x_t,t)$ denote a neural CFM, and $v_t^\star,p_t^\star$ denote the true velocity field and distribution respectively at time $t$. Specifically, we want to test if there exists a kernel $W$ such that

closely approximates $\hat{v}_\theta$.

To assess the validity of the equivalent-kernel assumption, we use a 1D Gaussian target:

with linear path $X_t=(1-t)X_0+tX_1$, giving

The factorized kernel approximation is most accurate when the slice-design density $\rho_{t,x}(c^\prime)=p^\star(x\mid c^\prime)p_C(c^\prime)$ varies little over the smoothing bandwidth $h_\lambda$. With $p_C$ locally constant, the path mean shifts by $\approx th_\lambda\lvert f^{\star\prime}(c)\rvert$ over this neighborhood while the path std is $s_t$, giving the heuristic condition

For our experiments, we use $f^\star(c)=\sin(2\pi c)+\sin(10\pi c)$, $\sigma=0.05$, and train a CFM $\hat{v}_\theta(c,x_t,t)$, analyzing checkpoints $N_{\mathrm{steps}}\in\{1000,3000,5000\}$ (heavy/moderate/mild underfitting).

To test the equivalent-kernel (EK) approximation, we discretize $c$ on a uniform grid $\{c_j\}_{j=1}^N$ and define

where $v_t^\star$ and $p_t^\star$ are available in closed form. We parameterize $W_\phi$ as a learned-embedding Gibbs kernel,

where $h_\phi:[0,1]\to\mathbb{R}^2$ is a small MLP and $\sigma_K$ is learned. The positive diagonal matrix $D_\phi$ is obtained by symmetric Sinkhorn normalization so that $W_\phi\mathbf{1}=\mathbf{1}$. This construction ensures

Denoting $\hat{v}_{\theta,t}=\hat{v}_\theta(\cdot,\cdot,t)$, we fit $\phi$ by minimizing

where $\hat{p}_t(\cdot\mid c)$ is the marginal induced by integrating the learned field $\hat{v}_\theta$ from the Gaussian base distribution up to time $t$. Thus, we test whether the EK approximation explains the trained field on the region of state space actually visited by the learned model.

Let $\hat{p}_1(\cdot\mid c),\widetilde{p}_1(\cdot\mid c)$ be the endpoint distributions obtained by integrating $\hat{v}_\theta,\widetilde{v}_\phi$ from a shared $\mathcal{N}(0,1)$ base, with means/variances $\hat{f},\hat{\sigma}^2$ and $\widetilde{f},\widetilde{\sigma}^2$ respectively. Writing

we report the field error

as a function of $t$ in Figure 9, Column 1. We also plot $\widetilde{v}_{\phi,t}(c,x_t)$ and $\hat{v}_{\theta,t}(c,x_t)$ as a function of $c$ at $x_t=0$ for $t=0.25,0.5,0.7,0.9$ in Figure 9, Columns 2–5. As expected, Figure 9 shows that $\widetilde{v}_\phi$ is a good approximation of $\hat{v}_\theta$ over much of the learned trajectory, with the approximation becoming weaker in some regions as $t\to1$. In Figure 11, we also plot the endpoint distributions obtained by integrating $\hat{v}_\theta$ and $\widetilde{v}_\phi$. Interestingly, the endpoint distributions look very similar, and the vector-field approximation error close to $t=1$ does not affect the integral significantly. Finally, in Figure 10 we visualize how both $\hat{v}_\theta$ and $\widetilde{v}_\phi$ look like density-weighted smoothed functions of the true velocity field $v^\star$.

Image not extracted

Image not extracted

Image not extracted

Coverage Experiments

This section provides full constructions of the predictive sets and conformal baselines used in § 6, along with experimental details.

Coverage Set Creation

Flow predictive intervals (scalar outputs) In the Monte-Carlo policy-evaluation experiments, the target variable is one-dimensional. We write each example as $(x_i,c_i)$, where $c_i$ is the conditioning input and $x_i \in\mathbb{R}$ is the target value. Given a trained conditional flow with flow map $\Phi_{1}^{\theta}(\cdot\mid c)$ and base distribution $\rho_0$, the learned conditional distribution at condition $c$ is the law of

Because the output is one-dimensional, flow trajectories do not cross. The map $z \mapsto\Phi^\theta_1(z \mid c)$ is therefore order-preserving, and central predictive intervals can be obtained by transporting quantiles of the base distribution. For a nominal coverage level $\alpha\in(0,1)$, let

We define the raw flow predictive interval by

Equivalently, this is the central $\alpha$-mass interval of the learned one-dimensional conditional distribution. The corresponding interval width is

For a dataset $\mathcal{D}$, we report the empirical coverage $\widehat{\operatorname{Cov}}_\alpha(\mathcal{D})$ and the average interval width

KNN predictive sets (high-dimensional outputs) For high-dimensional outputs, central quantile intervals are not well-defined, and the standard conformal baselines (split conformal, CQR) do not generalize naturally to multivariate outputs. We instead use a sample-based KNN predictive set, similar to what is proposed in Wang et al. [40].

For a condition $c$, draw $\widetilde{x}_1,\ldots,\widetilde{x}_M \sim\widehat{p}_\theta(\cdot\mid c)$. For any candidate point $x$, define $d_k(x;c)$ as the distance from $x$ to its $k$-th nearest neighbor among the generated samples. We compute the same score for the generated samples themselves using leave-one-out nearest neighbors, and define

The KNN predictive set is

A held-out target $x_i$ is counted as covered if $d_k(x_i;c_i) \le\tau_{\theta,\alpha}(c_i)$. Throughout, we use $k=1$ for computing coverage metrics.

Baselines

All conformal methods use a held-out calibration set $\mathcal{D}_{\mathrm{cal}}$ disjoint from the training set. The deterministic baselines use the same backbone architecture as the flow’s conditioning encoder and are trained on the same training set with matched optimizer, learning rate, and number of gradient steps. Split conformal prediction is applied either to a deterministic regressor (det + SC) or to the conditional mean of the trained CFM (CFM mean + SC). Conf CFM instead applies conformal calibration directly to the raw predictive set produced by the CFM.

Split conformal prediction. Let $\hat{f}(c)$ denote a point predictor. For det + SC, $\hat{f}(c)$ is the output of the deterministic regressor, whereas for CFM mean + SC, $\hat{f}(c)$ is the conditional mean of the trained CFM.

On $\mathcal{D}_{\mathrm{cal}}$, we compute absolute residual scores

Let $q_{\alpha}^{\mathrm{split}}$ be the $\alpha$-quantile of $\{r_i\}$, with the standard finite-sample correction. The split-conformal interval is

with width $W_{\alpha}^{\mathrm{split}}(c) = 2q_{\alpha}^{\mathrm{split}}$. Every input receives the same expansion, so split conformal achieves marginal coverage but is not locally adaptive.

For high-dimensional outputs, CFM mean + SC uses the analogous construction with Euclidean residual scores

The resulting predictive set is the conformal ball

Conformalized quantile regression. For det + CQR, at each coverage level $\alpha$, we train two predictors $\hat{\ell}_{\alpha}(c)$ and $\hat{u}_{\alpha}(c)$ using quantile losses for the lower and upper tails, with target quantile levels $(1-\alpha)/2$ and $(1+\alpha)/2$, respectively. Before conformal calibration, these define the interval $[\hat{\ell}_{\alpha}(c),\hat{u}_{\alpha}(c)]$.

On $\mathcal{D}_{\mathrm{cal}}$, we compute CQR scores

Let $q_{\alpha}^{\mathrm{CQR}}$ be the $\alpha$-quantile of $\{s_i\}$, with the standard finite-sample correction. The conformalized CQR interval is

with width

Unlike split conformal prediction, CQR can produce different interval widths at different conditioning inputs.

Conformalized CFM. For scalar outputs, let

denote the raw CFM predictive interval. On $\mathcal{D}_{\mathrm{cal}}$, we compute the conformity scores

Let $q_{\alpha}^{\mathrm{CFM}}$ be their $\alpha$-quantile, with the standard finite-sample correction. The conformalized CFM interval is

For high-dimensional outputs, we conformalize the KNN predictive set. Let $\tau_{\theta,\alpha}(c)$ denote the raw CFM KNN threshold defined above. We compute calibration scores

and let $q_{\alpha}^{\mathrm{CFM}}$ be their conformal quantile. The resulting predictive set is

Thus, Conf CFM applies a conformal correction to the boundary of the raw CFM predictive set without changing the underlying learned distribution.

Deep Evidential Regression. We use the standard Deep Evidential Regression method of Amini et al. [5] as an additional probabilistic uncertainty baseline. DER is evaluated only for scalar outputs, and its predictive intervals are obtained from its learned predictive distribution. We report DER at its best validation checkpoint.

Heteroscedastic Gaussian Regression. Heteroscedastic Gaussian Regression models the conditional distribution as

For scalar outputs, the network predicts a conditional mean $\mu_{\psi}(c)$ and variance $\sigma_{\psi}^{2}(c) > 0$, and is trained by minimizing the Gaussian negative log-likelihood

up to an additive constant. The central $\alpha$-mass predictive interval is

where $z_q$ denotes the $q$-quantile of the standard normal distribution.

For $d$-dimensional outputs, HGR predicts a diagonal conditional covariance,

and minimizes

Computing $k$-NN radii across methods. For high-dimensional outputs, we express all methods’ predictive-region sizes on the same $k$-NN radius scale used for the CFM. At each condition $c$, the sample-based radius is the $\alpha$-quantile of the leave-one-out $k$-nearest-neighbor distances among $K$ samples from the predictive distribution; we report its average over conditions. We hold $K$ and $k$ fixed across methods within each task to ensure a common sampling resolution. This preserves the CFM’s original radius statistic while avoiding an ellipsoidal approximation to its predictive region.

For the CFM, we sample from the learned conditional flow; for HGR, from the predicted diagonal Gaussian. HGR radii are computed separately at each condition because Euclidean $k$-NN distances depend on the full shape of the diagonal covariance, not only on the geometric mean of its marginal standard deviations.

For the point-predictor split-conformal baselines (det + SC and CFM mean + SC), no per-condition sampling is required. At each nominal level $\alpha$, we interpret the conformal $L_2$ ball of radius $R_{\alpha}$ as the central $\alpha$-mass region of an isotropic Gaussian $\mathcal{N}\left(\mu(c), s_{\alpha}^{2}I_D\right)$, where $\mu(c)$ is the corresponding point prediction and $D$ is the output dimension. This fixes

where $\chi_{\alpha}^{2}(D)$ denotes the $\alpha$-quantile of the chi-squared distribution with $D$ degrees of freedom. Since pairwise distances are translation-invariant and scale linearly with $s_{\alpha}$, the corresponding radius is

where $q_{\alpha,k}(K,D)$ is the reference radius obtained by applying the same leave-one-out $k$-NN quantile procedure to a cloud of $K$ standard Gaussian samples in $D$ dimensions. We estimate this reference quantity once by simulation for each $(\alpha,k,K,D)$ and reuse it for both ball-based baselines. The distance rescaling is exact under the isotropic-Gaussian representation; only the reference quantity is estimated numerically. This representation is used solely to express the ball radius on the common $k$-NN scale.

For Conf CFM, we report the sample-based CFM threshold after adding the conformal correction defined above, again averaged over conditions.

Datasets and Training Setup

All the plots are created using 3 random seeds. Dataset specific details and training setup are given below.

Biology datasets

We evaluate conditional flow matching on three single-cell perturbation benchmarks from the harmonized scPerturb collection [27]: Norman, SciPlex3, and PBMC cytokines. In each dataset, a cell is represented by a PCA-reduced gene-expression vector $x \in\mathbb{R}^{d}$, and each sample is annotated with a perturbation condition $c$ together with a control condition. The perturbation prediction problem is cast as conditional generation: the source distribution is the population of control cells, the target distribution is the population of cells under the perturbation condition, and the conditioning variable $c$ is a pretrained embedding of the perturbation. With the source distribution of control cells fixed, the flow learns the distribution of perturbed-cell states conditioned on $c$.

Norman (Perturb-seq). The Norman et al. [24] Perturb-seq dataset consists of single-cell RNA-seq profiles of K562 cells subjected to CRISPRa-induced overexpression of one or two genes. We use the preprocessed AnnData release distributed by the authors of biolord [29], which contains $\sim$89k cells across 229 single and double-gene perturbation conditions plus unperturbed controls; this is the same release used by the CellFlow benchmark [17]. Of biolord’s predefined train/test partitions, we use the “train” perturbation conditions ($\sim$42k cells) for our flow-matching experiments and the in-distribution cells are then randomly split into a training set and a validation set (80-20 split). The 7353 unperturbed control cells form the pool from which $x_0$ samples are drawn at training time. The condition $c \in\mathbb{R}^{2560}$ is the concatenation of two ESM2-650M [20] protein embeddings, one per perturbed gene. The target $x \in\mathbb{R}^{50}$ is the 50-component PCA projection of the log-normalized expression vector, with PCA fit on all cells.

The $x_0$ distribution is control-pool + Gaussian noise: $x_0 = z + \sigma\epsilon$ with $z \sim\operatorname{Uniform}(\text{control pool})$, $\epsilon\sim\mathcal{N}(0,I)$, and $\sigma= 1.0$. The velocity field is a CellFlow-style [17] three-branch MLP: a condition encoder $([1024,1024] \to256)$, a cell-state encoder $([512,512,512] \to512)$, and a sinusoidal time embedding (128d) followed by $[512,512,512] \to512$, fused by a $[1024,1024,1024] \to50$ decoder. Activations are SiLU throughout. We train with Adam at learning rate $10^{-3}$, batch size 1024. Sampling at evaluation uses 100 Euler steps.

PBMC (cytokine stimulation). The PBMC cytokine atlas released by Parse Biosciences [26] contains single-cell RNA-seq profiles of peripheral blood mononuclear cells from 12 donors stimulated with 90 cytokines plus PBS controls ($\sim$9.7M cells over 2,000 highly-variable genes). We access this dataset through the CellFlow data loader [17] and restrict to a single donor (Donor 1), which contains $\sim$1.3m cells: $\sim$90k PBS-stimulated controls and $\sim$1.2m cytokine-stimulated cells across the 90 cytokines. To keep the per-condition cell count balanced, we subsample uniformly at random to at most 1,000 cells per cytokine condition (and the same cap for controls), yielding 91,000 cells which are split into a training set and a validation set (80–20 split). The 1000 subsampled PBS controls form the pool from which $x_0$ samples are drawn at training time. The condition $c \in\mathbb{R}^{2560}$ is the ESM2 [20] protein embedding of the cytokine. The target $x \in\mathbb{R}^{100}$ is the 100-component PCA projection of the log-normalized expression vector, with PCA fit on all cells.

The $x_0$ distribution is again control-pool + Gaussian noise with standard deviation $\sigma= 10.0$. The velocity field used a condition encoder $([1024,1024] \to256)$, a cell-state encoder $([512,512,512] \to$ 512), and a sinusoidal time embedding ($1024$d) followed by [1024, 1024, 1024]$\to$1024, fused by a [1024, 1024, 1024]$\to$100 decoder. We use SiLU activations, train with Adam at learning rate $10^{-3}$, batch size 1024. Sampling at evaluation uses 1,000 Euler steps.

SciPlex (small-molecule perturbations). The SciPlex3 dataset of Srivatsan et al. [36] profiles single-cell transcriptomes of three cancer cell lines (A549, MCF7, K562) exposed to $\sim$188 small-molecule drugs at four doses each, with DMSO as the vehicle control. We restrict to the A549 lung-adenocarcinoma line, retrieve the dataset through the pertpy library (Heumos et al., [14]) via pertpy.data.sciplex3_raw() and follow the preprocessing of the CellFlow benchmark (Klein et al., [17]). After filtering, the dataset contains 703 in-distribution conditions covering 187 drugs, with $\sim$136k perturbed cells and $\sim$3k DMSO control cells, which are split into train-val datasets (80-20 split). The condition $c \in\mathbb{R}^{1025}$ is [Morgan-1024-fingerprint(drug), $\log_{10}$(dose in nM)]. The target $x \in\mathbb{R}^{100}$ is the 100-component PCA projection of the log-normalized expression vector, with PCA fit on all cells.

Similar to the previous two datasets, the $x_0$ distribution is constructed by adding gaussian noise with standard-deviation = 10.0 on the control pool. The velocity field architecture is the same as that used in the Norman experiments. Activations are SiLU throughout. We train with Adam at learning rate $10^{-3}$, batch size 1024, for 1,000 epochs. Sampling at evaluation uses 1000 Euler steps.

Policy Evaluation

We construct three policy-evaluation datasets on top of the OGBench (Park et al., [25]) offline RL benchmark, one per environment. For each environment, we take the trained floq agent of Agrawalla et al. [1] as the policy $\pi$ to be evaluated, and distill its true return $Q^\pi(s,a)$ into a flow-matching predictor by Monte-Carlo: for each $(s,a)$ pair we average the discounted return over $n_{\mathrm{MC}}$ independent on-policy rollouts of $\pi$ from $(s,a)$, giving a stochastic target $\widehat{Q}(s,a) \in\mathbb{R}$. The condition is $c = [s,a]$ and the target is $x = \widehat{Q}(s,a)$.

The per-environment details are:

Dataset	$n_{\mathrm{MC}}$	$H$	$\gamma$	dim $s$	dim $a$	dim $c$	$\lvert D\rvert$
Cube-single	20	200	0.99	28	5	33	200,000
AntMaze-large	20	1000	0.998	29	8	37	200,000
Humanoid-Maze-medium	20	2000	0.999	69	21	90	200,000

$x_0$ distribution. All three FloQ runs use an isotropic Gaussian prior $x_0 \sim\mathcal{N}(0,1)$ in the 1-D target space.

Architecture and training. For all three FloQ datasets we use the same MLP velocity field with hidden widths [512, 512, 512, 512] and ReLU activations, taking the full input $[c, x_t, \operatorname{TimeEmb}(t)]$, where TimeEmb is a 32-dimensional cosine-geometric time embedding. The output is a scalar velocity. We optimize with Adam (lr $3 \times10^{-4}$, no weight decay) using batch size 1024.

Behavior Cloning

We train a flow-matching behavior-cloning policy on the state-based variant of the PushT benchmark introduced by Chi et al. [8], where a 2-D end-effector pushes a T-shaped block toward a fixed target pose. The observation is the 5-dimensional state [agent$_x$, agent$_y$, block$_x$, block$_y$, block$_\theta$] and the action is the 2-dimensional target end-effector position. Following the standard receding-horizon BC setup, each example predicts an 8-step action chunk from 2 steps of observation history. Both observations and actions are MinMax-normalized over the training set. We use an 80 − 20 train-val split. The condition $c \in\mathbb{R}^{10}$ is the flattened observation history, the target $x \in\mathbb{R}^{16}$ is the flattened action chunk and a 128-dimensional sinusoidal time embedding is used. The $x_0$ distribution is a standard gaussian.

The velocity field is an 8-block residual MLP with embedding width 512, feed-forward expansion factor 4 (i.e., a 2,048-dim FFN), GELU activations, LayerNorm, and dropout 0.1. Time is encoded with a 128-dim sinusoidal positional embedding with frequencies linear in $[0, 100]$ followed by an MLP. No image encoder is used (state-based variant). We train using Adam with learning rate $10^{-4}$, weight decay $10^{-5}$, gradient clipping at norm 10 and batch size 1024.

Additional Experiments

Policy Evaluation § 6 reports policy-evaluation results for Antmaze-Large in the main paper (Figure 6). The same analysis on the remaining two OGBench environments, Cube-Single and Humanoidmaze-medium, is shown in Figure 12 and Figure 13 below. Both figures use the layout described in § 6: the top row shows training loss and flow coverage at three nominal levels, and the bottom row shows the variance-vs-error scatter along with mean interval widths for the flow, split conformal, and CQR. The qualitative findings discussed in § 6, all hold consistently across these two additional environments.

Cell perturbation § 6 reports coverage results for Norman in the main paper (Figure 7). The same analysis on the remaining PBMC cytokines and scPerturb benchmark SciPlex3 (drug-dose perturbations), is shown in Figure 14 and Figure 15 respectively. The figures use the same layout as Figure 7: training and validation flow-matching loss alongside KNN-based coverage at the three nominal levels.

Image not extracted

Image not extracted

The qualitative findings discussed in § 6 hold here as well: training coverage remains at or above nominal throughout training, and validation coverage is close to nominal during the underfitting regime. The flow’s predictive sets are approximately calibrated on SciPlex3 and PBMC, consistent with the pattern observed on Norman.

Variance Inflation: Additional Synthetic Experiments

In this section, we present evidence of variance inflation trends similar to Figure 2 across several more synthetic datasets. Similar to Figure 2, for all the experiments below, the learned conditional variance increases as the learned mean becomes smoother. The average learned conditional variance decreases monotonically as smoothing decreases across all three estimators.

Image not extracted