On the Role of Strain and Vorticity in Numerical Integration Error for Flow Matching

# On the Role of Strain and Vorticity in Numerical Integration Error for Flow Matching
Source: [https://arxiv.org/html/2605.06680](https://arxiv.org/html/2605.06680)
Seung\-Kyum Choi Georgia Institute of Technology schoi@me\.gatech\.edu

###### Abstract

Flow matching generates data by integrating a learned velocity field, where the number of integration steps \(NFE\) directly determines inference cost\. Yet a precise understanding of*which properties of the velocity field govern integration error*has been lacking\. We provide such an understanding by decomposing the velocity Jacobian∇xv\\nabla\_\{x\}vinto its symmetric partSS\(strain rate\) and antisymmetric partΩ\\Omega\(vorticity\), and proving that they play fundamentally different roles: strain controls*exponential*error amplification via the logarithmic normμ2=λmax​\(S\)\\mu\_\{2\}=\\lambda\_\{\\max\}\(S\), while vorticity contributes only*linearly*to the local truncation error\. This asymmetry has three implications\. First, we derive a separated error bound showing that suppressing strain alone eliminates exponential error growth, while suppressing vorticity alone does not\. Second, we prove that the optimal transport velocity field is automatically irrotational \(Ω=0\\Omega=0\) and has zero material derivative, which upgrades Euler integration from first\-order to second\-order accuracy\. For exact OT displacement interpolation, the corresponding Lagrangian particle dynamics are in fact integrated exactly by Euler; we verify this on both Gaussian and nonlinear OT flows, where errors reach machine precision \(∼10−14\\sim\\\!10^\{\-14\}\)\. Third, we show that a weighted Jacobian regularizer with strain weightα\\alphaexceeding vorticity weightβ\\betais theoretically favored, a prediction we validate on synthetic benchmarks and probe on CIFAR\-10\. Experiments on 2D distributions confirm the main theoretical predictions, demonstrating up to2\.7×2\.7\\timesreduction in integration error at NFE==5\. Preliminary CIFAR\-10 experiments show consistent trends, with a lightweight fine\-tuning procedure yielding 14% FID improvement at NFE==10 while preserving high\-NFE quality\. A matched fine\-tuning control experiment \(same training, no regularization\) shows no comparable improvement, indicating that the gains are associated with Jacobian regularization rather than additional training alone\. Ablations further illustrate the predicted bias\-complexity tradeoff and support strain\-dominant weighting in the low\-dimensional setting\.

## 1Introduction

Flow Matching\(Lipman et al\.,[2023](https://arxiv.org/html/2605.06680#bib.bib1); Liu et al\.,[2023](https://arxiv.org/html/2605.06680#bib.bib2)\)has emerged as a powerful paradigm for generative modeling, training a velocity fieldvθ​\(t,x\)v\_\{\\theta\}\(t,x\)whose ODE integration transports noise to data\. A central practical challenge is that accurate integration requires many function evaluations \(high NFE\), making inference slow\. A variety of methods address this problem: Rectified Flow\(Liu et al\.,[2023](https://arxiv.org/html/2605.06680#bib.bib2)\)straightens trajectories via reflow, Consistency Models\(Song et al\.,[2023](https://arxiv.org/html/2605.06680#bib.bib3)\)enforce self\-consistency, and MeanFlow\(Geng et al\.,[2025](https://arxiv.org/html/2605.06680#bib.bib4)\)learns an average velocity enabling one\-step generation\.

Despite these practical advances, a fundamental question remains:*what properties of the learned velocity field determine how many integration steps are needed?*Standard numerical analysis bounds the Euler error using the Lipschitz constantL=sup‖∇xv‖L=\\sup\\\|\\nabla\_\{x\}v\\\|, yieldingO​\(h⋅eL​T\)O\(h\\cdot e^\{LT\}\)\. But this bound treats all components of the Jacobian equally, potentially missing structure that could yield tighter analysis and better\-targeted regularization\.

In this paper, we provide a finer\-grained analysis by decomposing the velocity Jacobian into itssymmetric partSS\(strain rate tensor\) andantisymmetric partΩ\\Omega\(vorticity tensor\)\. We prove that these two components affect integration error in fundamentally different ways:

- •Strain controls exponential error amplification\.The error propagation factor is governed by the logarithmic normμ2​\(∇v\)=λmax​\(S\)\\mu\_\{2\}\(\\nabla v\)=\\lambda\_\{\\max\}\(S\), which depends*only*onSS\. Large strain causes errors to grow aseμ\+​Te^\{\\mu\_\{\+\}T\}\.
- •Vorticity contributes only linearly\.Vorticity affects the local truncation error through the termΩ​v\\Omega v, but this contribution is purely additive — it does not enter the exponential amplification factor\.
- •Suppressing both yields the tightest Euler bound\.WhenS→0S\\to 0andΩ→0\\Omega\\to 0\(the vanishing\-strain\-and\-vorticity regime\), the error collapses fromO​\(h⋅eL​T\)O\(h\\cdot e^\{LT\}\)toO​\(h​T⋅Mt\)O\(hT\\cdot M\_\{t\}\), whereMt=sup‖∂tv‖M\_\{t\}=\\sup\\\|\\partial\_\{t\}v\\\|\.

We further connect this analysis to optimal transport theory, proving that the OT velocity field from Brenier’s theorem is irrotational \(Ω=0\\Omega=0\) and has zero material derivative \(D​v/D​t=0Dv/Dt=0\)\. The latter implies that Euler integration is automatically*second\-order*accurate on OT flows\.

Contributions\.

1. 1\.Aseparated error bound\(Theorem[1](https://arxiv.org/html/2605.06680#Thmtheorem1)\) proving the asymmetric roles of strain and vorticity in ODE integration error, together with a three\-regime analysis \(Corollary[2](https://arxiv.org/html/2605.06680#Thmtheorem2)\)\.
2. 2\.Proof that theOT velocity field is irrotational\(Theorem[4](https://arxiv.org/html/2605.06680#Thmtheorem4)\) and haszero material derivative\(Theorem[5](https://arxiv.org/html/2605.06680#Thmtheorem5)\), which yields second\-order Euler accuracy in the Eulerian error analysis\. For exact displacement interpolation, we further observe exact Lagrangian Euler integration on both Gaussian and nonlinear OT flows, a stronger phenomenon discussed in Remark[7](https://arxiv.org/html/2605.06680#Thmtheorem7)\.
3. 3\.Abias\-complexity tradeoff analysisshowing that strain regularization is theoretically more valuable than vorticity regularization for controlling Euler discretization error \(Proposition[8](https://arxiv.org/html/2605.06680#Thmtheorem8)\)\.
4. 4\.Experimental validationon 2D benchmarks confirming the main theoretical predictions, together with supporting CIFAR\-10 experiments that show consistent low\-NFE improvements and ablations overα\\alpha,β\\beta, and fine\-tuning duration\.

Figure[1](https://arxiv.org/html/2605.06680#S1.F1)illustrates the core idea: as both strain and vorticity are suppressed, particle trajectories become progressively straighter and non\-crossing, enabling accurate integration with fewer steps\.

![Refer to caption](https://arxiv.org/html/2605.06680v1/figures/hero_figure.png)Figure 1:Overview of the three flow regimes\. \(a\) Standard FM: chaotic, crossing trajectories with high strain and vorticity, requiring many integration steps\. \(b\) Vorticity suppressed \(Ω→0\\Omega\\to 0\): smoother but still curved trajectories; exponential error amplification persists due to strain\. \(c\) Double\-null limit \(S,Ω→0S,\\Omega\\to 0\): nearly straight, parallel trajectories; error growth is linear, enabling few\-step generation\. Error bounds from Theorem[1](https://arxiv.org/html/2605.06680#Thmtheorem1)and Corollary[2](https://arxiv.org/html/2605.06680#Thmtheorem2)are shown below each panel\.
## 2Background

Flow Matching\.Given source distributionp0=𝒩​\(0,I\)p\_\{0\}=\\mathcal\{N\}\(0,I\)and data distributionp1p\_\{1\}, flow matching\(Lipman et al\.,[2023](https://arxiv.org/html/2605.06680#bib.bib1)\)learnsvθ:\[0,1\]×ℝd→ℝdv\_\{\\theta\}:\[0,1\]\\times\\mathbb\{R\}^\{d\}\\to\\mathbb\{R\}^\{d\}by minimizing:ℒFM\(θ\)=𝔼t,x0,x1∥vθ\(t,xt\)−ut\(xt\|x1\)∥2\\mathcal\{L\}\_\{\\text\{FM\}\}\(\\theta\)=\\mathbb\{E\}\_\{t,x\_\{0\},x\_\{1\}\}\\\|v\_\{\\theta\}\(t,x\_\{t\}\)\-u\_\{t\}\(x\_\{t\}\|x\_\{1\}\)\\\|^\{2\}, wherext=\(1−\(1−σ\)​t\)​x0\+t​x1x\_\{t\}=\(1\-\(1\-\\sigma\)t\)x\_\{0\}\+tx\_\{1\}is the conditional OT interpolation\.

Euler Integration Error\.At inference, we solved​x/d​t=vθ​\(t,x\)\\mathrm\{d\}x/\\mathrm\{d\}t=v\_\{\\theta\}\(t,x\)with step sizeh=1/Nh=1/N:x^n\+1=x^n\+h⋅v​\(tn,x^n\)\\hat\{x\}\_\{n\+1\}=\\hat\{x\}\_\{n\}\+h\\cdot v\(t\_\{n\},\\hat\{x\}\_\{n\}\)\. The standard error bound is‖eN‖=O​\(h⋅\(eL​T−1\)/L\)\\\|e\_\{N\}\\\|=O\(h\\cdot\(e^\{LT\}\-1\)/L\), whereL=supt‖∇xv​\(t,⋅\)‖L=\\sup\_\{t\}\\\|\\nabla\_\{x\}v\(t,\\cdot\)\\\|is the Lipschitz constant\. This bound is tight but*pessimistic*: it treats all components of∇xv\\nabla\_\{x\}vuniformly\.

Jacobian Decomposition\.Any matrixAAdecomposes uniquely asA=S\+ΩA=S\+\\OmegawhereS=\(A\+A⊤\)/2S=\(A\+A^\{\\top\}\)/2is symmetric andΩ=\(A−A⊤\)/2\\Omega=\(A\-A^\{\\top\}\)/2is antisymmetric\. These are Frobenius\-orthogonal:‖A‖F2=‖S‖F2\+‖Ω‖F2\\\|A\\\|\_\{F\}^\{2\}=\\\|S\\\|\_\{F\}^\{2\}\+\\\|\\Omega\\\|\_\{F\}^\{2\}\.

Logarithmic Norm\.The logarithmic norm \(matrix measure\)μ2​\(A\)=λmax​\(SA\)\\mu\_\{2\}\(A\)=\\lambda\_\{\\max\}\(S\_\{A\}\)satisfies‖et​A‖≤et​μ2​\(A\)\\\|e^\{tA\}\\\|\\leq e^\{t\\mu\_\{2\}\(A\)\}\. Crucially,μ2\\mu\_\{2\}depends*only on the symmetric part*ofAA, entirely ignoring the antisymmetric part\. This classical result\(Söderlind,[2006](https://arxiv.org/html/2605.06680#bib.bib8)\)is the foundation of our analysis\.

## 3Main Results: Separated Error Bound

### 3\.1Asymmetric Roles of Strain and Vorticity

We define the key quantities along the flow: the supremal logarithmic normμ\+=suptλmax​\(St\)\\mu\_\{\+\}=\\sup\_\{t\}\\lambda\_\{\\max\}\(S\_\{t\}\); the temporal variationMt=sup‖∂tv‖M\_\{t\}=\\sup\\\|\\partial\_\{t\}v\\\|; the strain\-induced accelerationMS=sup‖S​v‖M\_\{S\}=\\sup\\\|Sv\\\|; and the vorticity\-induced accelerationMΩ=sup‖Ω​v‖M\_\{\\Omega\}=\\sup\\\|\\Omega v\\\|\.

###### Theorem 1\(Global Error with Jacobian Decomposition\)\.

Under standard regularity assumptions \(vvisC2C^\{2\}, uniformly Lipschitz\), the Euler global error satisfies:

‖eN‖≤h​\(Mt\+MS\+MΩ\)2​μ\+​\(eμ\+​T−1\)\+O​\(h2\)\.\\\|e\_\{N\}\\\|\\leq\\frac\{h\(M\_\{t\}\+M\_\{S\}\+M\_\{\\Omega\}\)\}\{2\\mu\_\{\+\}\}\\left\(e^\{\\mu\_\{\+\}T\}\-1\\right\)\+O\(h^\{2\}\)\.\(1\)

The proof proceeds in three steps\. First, the error recursionen\+1=\(I\+h​∇v\)​en−τne\_\{n\+1\}=\(I\+h\\nabla v\)e\_\{n\}\-\\tau\_\{n\}is analyzed using the logarithmic norm to bound‖I\+h​∇v‖≤eh​μ\+\\\|I\+h\\nabla v\\\|\\leq e^\{h\\mu\_\{\+\}\}, which depends only onSS\(notΩ\\Omega\)\. Second, the local truncation error is decomposed via the material derivative:τn=h22​\(∂tv\+S​v\+Ω​v\)\+O​\(h3\)\\tau\_\{n\}=\\frac\{h^\{2\}\}\{2\}\(\\partial\_\{t\}v\+Sv\+\\Omega v\)\+O\(h^\{3\}\), separating contributions fromSSandΩ\\Omega\. Third, the discrete Grönwall lemma combines these to yield \([1](https://arxiv.org/html/2605.06680#S3.E1)\)\. Full proof in Appendix[A](https://arxiv.org/html/2605.06680#A1)\.

The key insightis the asymmetry:SSappears in*both*the exponential amplification factor \(eμ\+​Te^\{\\mu\_\{\+\}T\}\) and the truncation error \(MSM\_\{S\}\), whileΩ\\Omegaappears*only*in the truncation error \(MΩM\_\{\\Omega\}\)\. This has immediate consequences:

###### Corollary 2\(Three Regularization Regimes\)\.

\(A\) Strain suppression only \(S→0S\\to 0,Ω\\Omegaarbitrary\):‖eN‖≤h​T2​\(Mt\+MΩ\)\\\|e\_\{N\}\\\|\\leq\\frac\{hT\}\{2\}\(M\_\{t\}\+M\_\{\\Omega\}\)\. Exponential growtheliminated, but error still depends onΩ\\Omega\.

\(B\) Vorticity suppression only \(Ω→0\\Omega\\to 0,SSarbitrary\):‖eN‖≤h​\(Mt\+MS\)2​μ\+​\(eμ\+​T−1\)\\\|e\_\{N\}\\\|\\leq\\frac\{h\(M\_\{t\}\+M\_\{S\}\)\}\{2\\mu\_\{\+\}\}\(e^\{\\mu\_\{\+\}T\}\-1\)\. Exponential growthpersists\.

\(C\) Double\-null limit \(S→0S\\to 0andΩ→0\\Omega\\to 0\):

∥eN∥≤h​T2⋅Mt\+O\(h2\)\.\\boxed\{\\\|e\_\{N\}\\\|\\leq\\frac\{hT\}\{2\}\\cdot M\_\{t\}\+O\(h^\{2\}\)\.\}\(2\)

Regime C is the tightest Euler bound within our decomposition: linear inhh, linear inTT, controlled only by the velocity field’s intrinsic time variation\. In terms of NFE complexity, achieving accuracyϵ\\epsilonrequiresO​\(eL​T/ϵ\)O\(e^\{LT\}/\\epsilon\)steps for standard FM, but onlyO​\(Mt​T/ϵ\)O\(M\_\{t\}T/\\epsilon\)in the vanishing\-strain\-and\-vorticity regime — an exponential\-to\-linear reduction\.

## 4Connection to Optimal Transport

We now establish that the optimal transport velocity field naturally satisfies half of the vanishing\-strain\-and\-vorticity condition, and enjoys a surprising additional property\.

###### Theorem 4\(OT Velocity Field is Irrotational\)\.

LetT∗=∇ΨT^\{\*\}=\\nabla\\Psibe the Brenier OT map withΨ∈C3\\Psi\\in C^\{3\}strictly convex\. The Eulerian velocity field of McCann’s displacement interpolation satisfiesΩO​T​\(t,y\)=0\\Omega^\{OT\}\(t,y\)=0for allt∈\[0,1\)t\\in\[0,1\)\.

###### Proof sketch\.\.

The Jacobian ofvO​Tv^\{OT\}is\(∇2Ψ−I\)​\[\(1−t\)​I\+t​∇2Ψ\]−1\(\\nabla^\{2\}\\Psi\-I\)\[\(1\-t\)I\+t\\nabla^\{2\}\\Psi\]^\{\-1\}\. Since∇2Ψ\\nabla^\{2\}\\Psiis symmetric, both factors are polynomials in a symmetric matrix, hence commute, and their product is symmetric\. ThereforeΩ=0\\Omega=0\. Full proof in Appendix[B](https://arxiv.org/html/2605.06680#A2)\. ∎

###### Theorem 5\(OT Flow Has Zero Material Derivative\)\.

The OT velocity field satisfiesD​vO​TD​t=∂tv\+\(∇v\)​v=0\\frac\{Dv^\{OT\}\}\{Dt\}=\\partial\_\{t\}v\+\(\\nabla v\)v=0for allt∈\[0,1\)t\\in\[0,1\)\. Consequently:

*\(i\)*The local truncation error isτn=O​\(h3\)\\tau\_\{n\}=O\(h^\{3\}\)instead ofO​\(h2\)O\(h^\{2\}\)\.

*\(ii\)*Euler integration achievessecond\-orderglobal convergence:‖eN‖≤C​h2⋅\(eμ\+​T−1\)/μ\+\\\|e\_\{N\}\\\|\\leq Ch^\{2\}\\cdot\(e^\{\\mu\_\{\+\}T\}\-1\)/\\mu\_\{\+\}\.

###### Proof sketch\.\.

In Lagrangian coordinates, each particle has constant velocityφ˙t​\(x\)=∇Ψ​\(x\)−x\\dot\{\\varphi\}\_\{t\}\(x\)=\\nabla\\Psi\(x\)\-x, independent oftt\. The accelerationφ¨t=0\\ddot\{\\varphi\}\_\{t\}=0equals the material derivative in Eulerian coordinates\. TheO​\(h2\)O\(h^\{2\}\)truncation error term vanishes, leavingO​\(h3\)O\(h^\{3\}\)\. ∎

Combining these results yields a hierarchy:

∇v=0⏟vanishing strain & vorticity⊊Ω=0⏟OT \(irrotational\)⊊General⏟Standard FM\\underbrace\{\\nabla v=0\}\_\{\\text\{vanishing strain \\& vorticity\}\}\\;\\subsetneq\\;\\underbrace\{\\Omega=0\}\_\{\\text\{OT \(irrotational\)\}\}\\;\\subsetneq\\;\\underbrace\{\\text\{General\}\}\_\{\\text\{Standard FM\}\}\(3\)with corresponding error scalings that improve fromO​\(h⋅eL​T/L\)O\(h\\cdot e^\{LT\}/L\)in the general case, toO​\(h2​eμ\+​T/μ\+\)O\(h^\{2\}e^\{\\mu\_\{\+\}T\}/\\mu\_\{\+\}\)for OT flows, and further toO​\(h​T​Mt\)O\(hTM\_\{t\}\)in the vanishing\-strain\-and\-vorticity regime\.

![Refer to caption](https://arxiv.org/html/2605.06680v1/figures/hierarchy.png)Figure 2:Hierarchy of velocity field regularity\. General FM \(outer\) has both strain and vorticity, with first\-order Euler convergence\. OT flows \(middle\) are irrotational \(Ω=0\\Omega=0\) with zero material derivative, upgrading Euler to second\-order convergence \(O​\(h2\)O\(h^\{2\}\)global error\)\. The double\-null limit \(inner\) additionally eliminates strain \(S=0S=0\), yieldingO​\(h​T\)O\(hT\)linear error growth\. Arrows indicate the effect of each regularization penalty\.
## 5Implications for Regularization

Our analysis provides principled guidance for designing regularizers to reduce NFE\.

### 5\.1Weighted Jacobian Regularization

Consider augmenting the FM loss with a weighted Jacobian penalty:

ℒ=ℒFM\+α​𝔼​\[‖S‖F2\]\+β​𝔼​\[‖Ω‖F2\]\.\\mathcal\{L\}=\\mathcal\{L\}\_\{\\text\{FM\}\}\+\\alpha\\,\\mathbb\{E\}\[\\\|S\\\|\_\{F\}^\{2\}\]\+\\beta\\,\\mathbb\{E\}\[\\\|\\Omega\\\|\_\{F\}^\{2\}\]\.\(4\)
###### Proposition 8\(Design Principle: Strain\-Dominant Weighting\)\.

For fixed regularization budgetα\+β=λ\\alpha\+\\beta=\\lambda, the allocation minimizing total error \(regularization bias \+ discretization error\) satisfiesα∗\>λ/2\>β∗\\alpha^\{\*\}\>\\lambda/2\>\\beta^\{\*\}wheneverμ\+​T\>1\\mu\_\{\+\}T\>1\. This follows from the exponential sensitivity of discretization error toμ\+\\mu\_\{\+\}\(controlled byα\\alpha\) versus the linear sensitivity toMΩM\_\{\\Omega\}\(controlled byβ\\beta\)\.

We state this as a design principle rather than an optimization theorem, since the precise mapping from\(α,β\)\(\\alpha,\\beta\)to post\-training values of\(μ\+,MS,MΩ\)\(\\mu\_\{\+\},M\_\{S\},M\_\{\\Omega\}\)depends on the optimization landscape\. We validate the prediction empirically in Section[7](https://arxiv.org/html/2605.06680#S7)\.

Whenα=β\\alpha=\\beta, the regularizer simplifies toα​‖∇v‖F2\\alpha\\\|\\nabla v\\\|\_\{F\}^\{2\}, computable via a single Hutchinson VJP:‖∇v‖F2=𝔼z​\[‖\(∇v\)⊤​z‖2\]\\\|\\nabla v\\\|\_\{F\}^\{2\}=\\mathbb\{E\}\_\{z\}\[\\\|\(\\nabla v\)^\{\\top\}z\\\|^\{2\}\],z∼𝒩​\(0,I\)z\\sim\\mathcal\{N\}\(0,I\)\.

### 5\.2Frobenius vs\. Spectral Norm

A subtlety is that the error bound depends onλmax​\(S\)\\lambda\_\{\\max\}\(S\)\(spectral\), while the regularizer penalizes‖S‖F\\\|S\\\|\_\{F\}\(Frobenius\)\. These are related byλmax​\(S\)≤‖S‖F≤d​λmax​\(S\)\\lambda\_\{\\max\}\(S\)\\leq\\\|S\\\|\_\{F\}\\leq\\sqrt\{d\}\\,\\lambda\_\{\\max\}\(S\)\. Reducing‖S‖F2\\\|S\\\|\_\{F\}^\{2\}is therefore a*sufficient*condition for reducingμ\+\\mu\_\{\+\}, though not tight: in the worst case \(d=3072d=3072for CIFAR\-10\), the gap isd≈55\\sqrt\{d\}\\approx 55\. However, for velocity fields with approximately isotropic Jacobian spectra \(as observed empirically\), the effective gap is much smaller\. On our 2D experiments, we measured bothλmax​\(S\)\\lambda\_\{\\max\}\(S\)and‖S‖F\\\|S\\\|\_\{F\}directly and found them tightly correlated \(R2\>0\.95R^\{2\}\>0\.95\), validating that the Frobenius penalty effectively controls spectral amplification in practice\. On CIFAR\-10, we measured‖S‖F≈180\\\|S\\\|\_\{F\}\\approx 180along trajectories \(Section[7](https://arxiv.org/html/2605.06680#S7)\); estimatingλmax​\(S\)\\lambda\_\{\\max\}\(S\)at this scale requires spectral methods \(e\.g\., power iteration\) and is left to future work\.

### 5\.3Vorticity Penalty as Soft OT Constraint

Theorem[4](https://arxiv.org/html/2605.06680#Thmtheorem4)reveals a principled interpretation of the vorticity penaltyβ​‖Ω‖F2\\beta\\\|\\Omega\\\|\_\{F\}^\{2\}: since the OT velocity field satisfiesΩO​T=0\\Omega^\{OT\}=0, its enstrophy \(integrated‖Ω‖F2\\\|\\Omega\\\|\_\{F\}^\{2\}\) achieves the global minimum of zero\. Any velocity field with nonzero vorticity has strictly higher enstrophy\. The vorticity penalty therefore acts as a*soft optimal transport constraint*, encouraging the learned flow toward the irrotational structure of the Brenier solution\. We note that this is distinct from the Benamou–Brenier result, which shows OT minimizes*kinetic energy*\(∫‖v‖2\\int\\\|v\\\|^\{2\}\); the zero\-enstrophy property follows separately from Theorem[4](https://arxiv.org/html/2605.06680#Thmtheorem4)\.

### 5\.4Gradient\-Field Parameterization

An alternative to explicit regularization is to parameterizevθ=∇xϕθv\_\{\\theta\}=\\nabla\_\{x\}\\phi\_\{\\theta\}for a scalar potentialϕθ\\phi\_\{\\theta\}\. This enforcesΩ≡0\\Omega\\equiv 0by construction, with the Jacobian∇v=∇2ϕ\\nabla v=\\nabla^\{2\}\\phiautomatically symmetric\. This connects to Brenier’s theorem: the OT mapT∗=∇ΨT^\{\*\}=\\nabla\\Psiimplies the OT velocity is a gradient field\.

### 5\.5Normalized Regularization Weights

In high\-dimensional settings \(d≫1d\\gg 1\),‖∇v‖F2=O​\(d\)\\\|\\nabla v\\\|\_\{F\}^\{2\}=O\(d\)whileℒFM=O​\(1\)\\mathcal\{L\}\_\{\\text\{FM\}\}=O\(1\), so the raw weightα\\alphamust scale asO​\(1/d\)O\(1/d\)\. We recommend the normalized parameterizationα~=α⋅d\\tilde\{\\alpha\}=\\alpha\\cdot d, which is comparable across dimensions\.

## 6Gaussian Case: Exact Computation

Forp0=𝒩​\(0,I\)p\_\{0\}=\\mathcal\{N\}\(0,I\)andp1=𝒩​\(μ1,Σ1\)p\_\{1\}=\\mathcal\{N\}\(\\mu\_\{1\},\\Sigma\_\{1\}\), all quantities admit closed\-form expressions\. The OT velocity Jacobian has eigenvalues\(σi−1\)/\(\(1−t\)\+t​σi\)\(\\sigma\_\{i\}\-1\)/\(\(1\-t\)\+t\\sigma\_\{i\}\)whereσi=λi​\(Σ1\)\\sigma\_\{i\}=\\sqrt\{\\lambda\_\{i\}\(\\Sigma\_\{1\}\)\}\. The strain norm is:

‖SO​T​\(t\)‖F2=∑i=1d\(σi−1\(1−t\)\+t​σi\)2\.\\\|S^\{OT\}\(t\)\\\|\_\{F\}^\{2\}=\\sum\_\{i=1\}^\{d\}\\left\(\\frac\{\\sigma\_\{i\}\-1\}\{\(1\-t\)\+t\\sigma\_\{i\}\}\\right\)^\{2\}\.\(5\)
Key observations: \(i\)ΩO​T≡0\\Omega^\{OT\}\\equiv 0\(confirming Theorem[4](https://arxiv.org/html/2605.06680#Thmtheorem4)\); \(ii\)SO​T=0S^\{OT\}=0iffΣ1=I\\Sigma\_\{1\}=I\(the vanishing\-strain\-and\-vorticity regime is achievable only for translated distributions\); \(iii\) strain is highest neart=0t=0andt=1t=1, motivating time\-dependent regularization schedules\. These exact expressions serve as a sanity check for our general theory and provide concrete intuition for the strain/vorticity decomposition\.

Empirical verification of Eulerian second\-order structure and exact Lagrangian OT integration\.We verify Theorem[5](https://arxiv.org/html/2605.06680#Thmtheorem5)by measuring Euler error at varying step sizeshhon exact OT velocity fields with non\-OT controls\. We test two settings: \(i\) Gaussian targets \(p1=𝒩​\(μ1,Σ1\)p\_\{1\}=\\mathcal\{N\}\(\\mu\_\{1\},\\Sigma\_\{1\}\)\), where the OT velocity is affine; and \(ii\) nonlinear targets with potentialΨ​\(x\)=12​‖x‖2\+ε4​∑ixi4\\Psi\(x\)=\\frac\{1\}\{2\}\\\|x\\\|^\{2\}\+\\frac\{\\varepsilon\}\{4\}\\sum\_\{i\}x\_\{i\}^\{4\}\(ε∈\{0\.3,0\.5\}\\varepsilon\\in\\\{0\.3,0\.5\\\}\), giving a genuinely nonlinear OT mapT∗​\(x\)=x\+ε​x3T^\{\*\}\(x\)=x\+\\varepsilon x^\{3\}with non\-constant Hessian∇2Ψ=I\+3​ε​diag​\(x2\)\\nabla^\{2\}\\Psi=I\+3\\varepsilon\\,\\mathrm\{diag\}\(x^\{2\}\)\. Non\-OT controls add a time\-dependent rotational perturbation to the OT velocity\.

Gaussian case\(Figure[3](https://arxiv.org/html/2605.06680#S6.F3)\)\. The affine OT velocity yields Euler errors at machine precision \(∼10−14\\sim\\\!10^\{\-14\}\) for all step sizes, consistent with exact integration since all truncation error terms vanish\. Non\-OT flows show standard first\-order convergence \(slope≈1\.0\\approx 1\.0\)\. The gap exceeds1012×10^\{12\}\\times\.

Nonlinear case\(Figure[4](https://arxiv.org/html/2605.06680#S6.F4)\)\. Despite the nonlinear OT map, Euler errors again reach machine precision across all tested configurations \(d∈\{2,5\}d\\in\\\{2,5\\\},ε∈\{0\.3,0\.5\}\\varepsilon\\in\\\{0\.3,0\.5\\\}\)\. Non\-OT flows consistently exhibit slope≈1\.0\\approx 1\.0\. This stronger\-than\-predicted result is explained by a property of displacement interpolation that goes beyond Theorem[5](https://arxiv.org/html/2605.06680#Thmtheorem5): under OT, each particle’s Lagrangian velocityφ˙t​\(x\)=∇Ψ​\(x\)−x\\dot\{\\varphi\}\_\{t\}\(x\)=\\nabla\\Psi\(x\)\-xis time\-independent, so the Euler method applied to the particle trajectory is exact at*any*step size — not merely second\-order\. The Eulerian coordinate transformation introduces only floating\-point\-level errors\. TheO​\(h2\)O\(h^\{2\}\)regime predicted by Theorem[5](https://arxiv.org/html/2605.06680#Thmtheorem5)should be interpreted as the generic Eulerian guarantee; the machine\-precision behavior here is a stronger consequence of exact displacement interpolation\. For learned neural velocity fields that only approximately satisfy the zero\-material\-derivative condition, one would generally expect the second\-order regime rather than exact integration\.

![Refer to caption](https://arxiv.org/html/2605.06680v1/figures/euler_convergence.png)Figure 3:Gaussian OT verification\. \(a\) Euler errors reach machine precision \(∼10−14\\sim\\\!10^\{\-14\}\) for all step sizes, since the affine OT velocity yields zero truncation error at all orders\. Non\-OT flows show first\-order convergence \(slope≈1\\approx 1\)\. The gap exceeds1012×10^\{12\}\\times\. \(b\) Same pattern atd=10d\\\!=\\\!10\.![Refer to caption](https://arxiv.org/html/2605.06680v1/figures/nonlinear_ot_convergence.png)Figure 4:Nonlinear OT verification \(Ψ​\(x\)=12​‖x‖2\+ε4​∑xi4\\Psi\(x\)=\\frac\{1\}\{2\}\\\|x\\\|^\{2\}\+\\frac\{\\varepsilon\}\{4\}\\sum x\_\{i\}^\{4\}\)\. Despite the nonlinear OT map, Euler errors reach machine precision in all configurations, while non\-OT flows followO​\(h\)O\(h\)\(slope≈1\\approx 1\)\. This confirms that the zero\-material\-derivative property of displacement interpolation yields exact Lagrangian integration, a result strictly stronger than the second\-order bound of Theorem[5](https://arxiv.org/html/2605.06680#Thmtheorem5)\.
## 7Experiments

Our experiments validate the theoretical predictions of Sections[3](https://arxiv.org/html/2605.06680#S3)–[5](https://arxiv.org/html/2605.06680#S5)\. We focus on three questions: \(Q1\) Does reducing‖S‖F2\\\|S\\\|\_\{F\}^\{2\}via regularization actually reduce Euler integration error? \(Q2\) Is strain suppression more valuable than vorticity suppression? \(Q3\) Do these effects transfer to realistic settings?

### 7\.12D Synthetic Experiments

Setup\.We train flow matching models on a pinwheel distribution \(5 arms with radial twist\), using a 5\-layer MLP \(256 hidden units\), 8000 epochs, batch size 512\. Integration error is measured as the L2 distance between Euler samples at a given NFE and reference samples generated with NFE==500\. All methods are monitored for‖S‖F2\\\|S\\\|\_\{F\}^\{2\}and‖Ω‖F2\\\|\\Omega\\\|\_\{F\}^\{2\}via exact Jacobian computation\. The regularization weights are chosen to probe the theory rather than to optimize a single scalar metric: we first sweepα\\alphawithβ=0\\beta=0to isolate the effect of strain suppression, then compare matched mixed settings such as\(α,β\)=\(0\.1,0\.1\)\(\\alpha,\\beta\)=\(0\.1,0\.1\)and\(0\.3,0\.05\)\(0\.3,0\.05\)to distinguish isotropic Jacobian penalization from strain\-dominant weighting at comparable budget\.

Q1: Strain reduction→\\toerror reduction\.Table[1](https://arxiv.org/html/2605.06680#S7.T1)shows a systematic sweep of strain regularization weightα\\alpha\(withβ=0\\beta=0\)\. Asα\\alphaincreases,‖S‖F2\\\|S\\\|\_\{F\}^\{2\}decreases monotonically from 1\.93 to 0\.45, and the Euler L2 error at NFE==5 drops correspondingly from 0\.63 to 0\.23 — a2\.7×2\.7\\timesimprovement\. FM loss increases modestly \(3\.17→\\to3\.60\), confirming the bias\-complexity tradeoff of Proposition[8](https://arxiv.org/html/2605.06680#Thmtheorem8)\. The optimalα\\alphadepends on the target NFE: largerα\\alphafavors low\-NFE regimes\.

Table 1:Effect of strain regularization weightα\\alphaon 2D pinwheel \(β=0\\beta=0\)\. L2@kk: Euler integration error vs\. NFE==500 reference\.NFE comparison\.Figure[5](https://arxiv.org/html/2605.06680#S7.F5)shows Euler error and Sliced Wasserstein distance across NFE values for five methods: FM baseline, VFM \(α=β=0\.1\\alpha=\\beta=0\.1\), VFM \(α=0\.3,β=0\.05\\alpha=0\.3,\\beta=0\.05\), VFM \(α=0\.5,β=0\\alpha=0\.5,\\beta=0\), and a gradient\-field model \(v=∇ϕv=\\nabla\\phi,α=0\.1\\alpha=0\.1\)\. VFM models consistently achieve lower error, with the gap most pronounced at low NFE\. At NFE==5, VFM \(α=0\.5\\alpha=0\.5\) achieves error comparable to the baseline at NFE≈\\approx20, a4×4\\timesreduction in required steps\.

Q2: Strain\>\>vorticity\.The gradient\-field model \(Ω≡0\\Omega\\equiv 0by construction\) achieves the highest trajectory straightness but*not*the lowest integration error\. VFM \(α=0\.5\\alpha=0\.5,β=0\\beta=0\) performs better, confirming Corollary[2](https://arxiv.org/html/2605.06680#Thmtheorem2): suppressing strain \(Regime A\) is more important than suppressing vorticity \(Regime B\)\. This is visible in the training curves \(Figure[6](https://arxiv.org/html/2605.06680#S7.F6)\), where‖S‖F2\\\|S\\\|\_\{F\}^\{2\}shows clear separation across methods while‖Ω‖F2\\\|\\Omega\\\|\_\{F\}^\{2\}remains uniformly small\.

![Refer to caption](https://arxiv.org/html/2605.06680v1/figures/nfe_error_comparison.png)Figure 5:NFE vs\. integration error on 2D pinwheel\. Left: L2 error vs\. NFE==500 reference\. Middle: Sliced Wasserstein distance\. Right: Trajectory straightness\. VFM models \(warm colors\) consistently achieve lower error than the FM baseline \(blue\), with the gap largest at low NFE\. At NFE==5, VFM \(α=0\.5\\alpha=0\.5\) matches the baseline at NFE≈\\approx20\.![Refer to caption](https://arxiv.org/html/2605.06680v1/figures/training_curves.png)Figure 6:Training dynamics on 2D pinwheel\. Left: FM loss\. Middle:‖S‖F2\\\|S\\\|\_\{F\}^\{2\}\(strain\)\. Right:‖Ω‖F2\\\|\\Omega\\\|\_\{F\}^\{2\}\(vorticity\)\. VFM with largerα\\alphaachieves lower strain throughout training\. Vorticity is naturally small \(∼10−3\\sim\\\!10^\{\-3\}\) for all methods — two orders of magnitude below strain — confirming that learned velocity fields are approximately irrotational\.
### 7\.2CIFAR\-10 Experiments

Setup\.To test whether the theoretical predictions extend to high\-dimensional settings, we train a SimpleUNet \(∼\\sim27M parameters\) on CIFAR\-10 for 200 epochs using standard flow matching \(FM baseline, final FM loss 0\.174\)\. We then fine\-tune with VFM regularization at reduced learning rate \(5×10−55\\times 10^\{\-5\}\)\. FID is computed on 50K generated samples\. We ablate over the regularization weightα\\alpha, the inclusion of vorticity penaltyβ\\beta, and the fine\-tuning duration\. The scale ofα\\alphais chosen using the dimensional normalization discussed in Section[5](https://arxiv.org/html/2605.06680#S5): since‖∇v‖F2=O​\(d\)\\\|\\nabla v\\\|\_\{F\}^\{2\}=O\(d\), weights of order10−610^\{\-6\}on CIFAR\-10 correspond to normalized strengthsα~=O​\(10−3\)\\tilde\{\\alpha\}=O\(10^\{\-3\}\), small enough to avoid overwhelming the FM objective while still producing a measurable change in Jacobian statistics\. We useβ=0\\beta=0to isolate strain regularization,β=α\\beta=\\alphafor equal\-budget comparisons against the standard Jacobian penalty, andα=0,β=10−6\\alpha=0,\\beta=10^\{\-6\}as a diagnosticβ\\beta\-only control\.

Table 2:CIFAR\-10 unconditional generation: FID \(↓\\downarrow\) at various NFE\. All fine\-tuned variants start from the same FM baseline \(200 epochs\) and train for 30 additional epochs at lr==5×\\times10\-5\(unless noted\)\. Bold: improvement over baseline\.Shaded: best configuration\.Q3: Transfer to high dimensions\.Table[2](https://arxiv.org/html/2605.06680#S7.T2)shows FID across configurations\. The critical finding is thefine\-tune control: 30 epochs of additional training*without*regularization produces FID indistinguishable from the original baseline \(e\.g\., 26\.0 vs\. 25\.8 at NFE==10\)\. This strongly suggests that the improvements observed in VFM fine\-tuning are due to Jacobian regularization rather than to additional training alone\.

The best overall configuration \(α=10−6\\alpha=10^\{\-6\},β=0\\beta=0, 30 epochs\) improves FID at NFE==10 from 26\.0 to 22\.2 \(−\-15% vs\. control\) and at NFE==50 from 15\.4 to 13\.6 \(−\-12%\), while preserving high\-NFE quality\.

Ablation insights\.

*\(i\) Regularization is necessary\.*The fine\-tune control rules out the hypothesis that improvement comes from additional training\. All VFM variants outperform the control at low\-to\-mid NFE, supporting the interpretation that Jacobian regularization is the main active ingredient in these gains\.

*\(ii\) Bias\-complexity tradeoff\.*Increasingα\\alphafrom10−610^\{\-6\}to2×10−62\\times 10^\{\-6\}improves low\-NFE FID \(46\.4 vs\. 47\.1 at NFE==5\) but degrades high\-NFE FID \(20\.2 vs\. 13\.9 at NFE==100\), confirming the tradeoff predicted by Proposition[8](https://arxiv.org/html/2605.06680#Thmtheorem8)\.

*\(iii\) Learned flows are naturally near\-irrotational\.*Theβ\\beta\-only configuration \(α=0\\alpha=0,β=10−6\\beta=10^\{\-6\}\) achieves Reg≈3×10−4\\approx 3\\times 10^\{\-4\}during training — three orders of magnitude smaller than theα\\alpha\-only Reg of4×10−24\\times 10^\{\-2\}— confirming that‖Ω‖F2≪‖S‖F2\\\|\\Omega\\\|\_\{F\}^\{2\}\\ll\\\|S\\\|\_\{F\}^\{2\}in learned velocity fields\. Yetβ\\beta\-only still improves over the control\. We interpret this cautiously: it may indicate an additional implicit regularization effect from the VJP\-based computational pathway \(which requiresrequires\_gradon intermediate states\), or another mechanism not captured by the present theory\. This discrepancy warrants further investigation\.

*\(iv\) Fine\-tuning duration matters\.*Extending from 30 to 50 epochs improves NFE==5 \(45\.9 vs\. 47\.1\) but degrades high\-NFE quality \(21\.0 vs\. 13\.9\), showing that the optimal fine\-tuning duration also follows a bias\-complexity tradeoff\.

*\(v\) Computational overhead\.*Each VFM fine\-tuning step requires one additional VJP \(viatorch\.autograd\.grad\) compared to standard FM\. Whenα=β\\alpha=\\beta, only‖∇v‖F2=𝔼z​\[‖J⊤​z‖2\]\\\|\\nabla v\\\|\_\{F\}^\{2\}=\\mathbb\{E\}\_\{z\}\[\\\|J^\{\\top\}z\\\|^\{2\}\]is needed \(1 VJP per probe vector, typically 1–2 probes\)\. Whenα≠β\\alpha\\neq\\beta, estimatingtr​\(J2\)\\mathrm\{tr\}\(J^\{2\}\)requires an additional finite\-difference JVP per probe \(1 extra forward pass\)\. Total overhead:∼2\.5×\\sim\\\!2\.5\\timestraining time for 1 probe \(∼\\sim200s/epoch vs\.∼\\sim80s/epoch on A5000\)\. Peak GPU memory increases by∼\\sim20% due to the computation graph retained for the VJP\. Inference cost is*identical*to standard FM — the regularizer is training\-only\.

Scope of results\.These results should be interpreted as a proof\-of\-concept\. Our model \(27M parameters, 200 epochs\) is significantly smaller than competitive baselines \(e\.g\., EDM2 at 280M–1\.1B parameters\)\. The purpose is to probe whether the theoretical picture remains informative in high dimensions, not to achieve the lowest possible FID\.

Dimensional scaling\.A notable practical finding is that the raw regularization weight must scale asα=O​\(1/d\)\\alpha=O\(1/d\)in high dimensions, since‖∇v‖F2=O​\(d\)\\\|\\nabla v\\\|\_\{F\}^\{2\}=O\(d\)\. On CIFAR\-10 \(d=3072d=3072\), the effective weightα=10−6\\alpha=10^\{\-6\}corresponds to normalizedα~=α​d≈0\.003\\tilde\{\\alpha\}=\\alpha d\\approx 0\.003, compared toα~≈0\.6\\tilde\{\\alpha\}\\approx 0\.6\(α=0\.3\\alpha=0\.3,d=2d=2\) in 2D\.

Jacobian statistics along trajectories\.We measured‖S‖F\\\|S\\\|\_\{F\}and‖Ω‖F\\\|\\Omega\\\|\_\{F\}along Euler trajectories for both the FM baseline and VFM fine\-tuned model using Hutchinson estimators\. Key findings: \(i\)∥S∥F/∥Ω∥F≈31×\\\|S\\\|\_\{F\}/\\\|\\Omega\\\|\_\{F\}\\approx 31\\timesfor both models, confirming that learned velocity fields are dominated by strain, consistent with 2D observations; \(ii\)‖S‖F\\\|S\\\|\_\{F\}increases from∼55\\sim\\\!55att=0t=0to∼208\\sim\\\!208att=0\.75t=0\.75, consistent with the Gaussian analysis predicting strain peaks near endpoints; \(iii\) VFM fine\-tuning reduces average‖S‖F\\\|S\\\|\_\{F\}by∼14%\\sim\\\!14\\%\(181→\\to156\), a modest but consistent reduction\. Direct measurement ofλmax​\(S\)\\lambda\_\{\\max\}\(S\)via power iteration at CIFAR\-10 scale \(d=3072d=3072\) is computationally expensive and is left to future work\.

### 7\.3Qualitative Sample Comparison

Figure[7](https://arxiv.org/html/2605.06680#S7.F7)compares unconditional CIFAR\-10 samples from the FM baseline, the matched fine\-tuning control, and the best VFM configuration across several NFEs using the same visualization protocol\. The clearest differences appear at low NFE, where the VFM model produces cleaner global structure and fewer obvious artifacts than either the FM baseline or the no\-regularization fine\-tuning control\. At higher NFE, sample quality is largely preserved, consistent with the quantitative trends in Table[2](https://arxiv.org/html/2605.06680#S7.T2)\.

![Refer to caption](https://arxiv.org/html/2605.06680v1/figures/cifar_qualitative_main.png)Figure 7:Qualitative CIFAR\-10 comparison across sampling budgets\. Rows show the FM baseline, the matched fine\-tuning control without regularization, and the best VFM model \(α=10−6,β=0\\alpha=10^\{\-6\},\\beta=0\)\. Columns correspond to NFE==5, 10, and 50\. The VFM model shows the most noticeable gains at low NFE, while maintaining competitive visual quality at higher NFE\.

## 8Related Work

Few\-step generation\.Rectified Flow\(Liu et al\.,[2023](https://arxiv.org/html/2605.06680#bib.bib2)\)straightens trajectories via iterative reflow; Sequential Reflow further refines this with coupling\-adapted objectives\. Consistency Models\(Song et al\.,[2023](https://arxiv.org/html/2605.06680#bib.bib3)\)and their extensions\(Geng et al\.,[2025](https://arxiv.org/html/2605.06680#bib.bib5)\)enforce temporal self\-consistency\. MeanFlow\(Geng et al\.,[2025](https://arxiv.org/html/2605.06680#bib.bib4)\)learns an average velocity field enabling one\-step generation via a JVP\-based identity\. Distillation methods \(Flow Generator Matching, Progressive Distillation\) train a student to mimic the teacher’s multi\-step output in a single step\. These methods are complementary to our analysis: our theoretical framework explains*why*straight trajectories help \(low strain→\\tono exponential error growth\) and predicts when they are insufficient \(high vorticity→\\topersistent truncation error\)\. We do not claim to outperform these methods, but rather to provide the analytical foundation that could inform their design\.

Training\-free sampling improvements\.Orthogonal to training\-time regularization, recent work improves ODE sampling through solver design: adaptive step\-size selection, backward\-error\-informed scheduling, and curvature\-aware solvers \(e\.g\., DPM\-Solver\)\. These methods reduce discretization error without modifying the velocity field\. Our analysis complements this line of work by identifying*which velocity field properties*determine solver performance\.

Jacobian regularization\.Jacobian norm penalties have been used for adversarial robustness\(Hoffman et al\.,[2019](https://arxiv.org/html/2605.06680#bib.bib9)\)and in neural ODEs for stability\(De Marinis et al\.,[2025](https://arxiv.org/html/2605.06680#bib.bib10)\)\. The latter work uses the logarithmic norm to control error growth in classification settings, establishing the relevance ofμ2\\mu\_\{2\}for neural ODE stability\. Our contribution builds on this by*decomposing*the Jacobian into strain and vorticity and proving their asymmetric roles specifically for generative ODE integration, yielding the predictionα\>β\\alpha\>\\beta\.

Regularization in flow\-based models\.Density\-weighted regularization\(Eguchi,[2025](https://arxiv.org/html/2605.06680#bib.bib13)\)suppresses velocity field oscillations in low\-density regions via a modified loss geometry\. Our approach is complementary: we regularize the Jacobian’s internal structure \(strain vs\. vorticity\) rather than the spatial distribution of‖v‖\\\|v\\\|\.

Numerical analysis of generative ODEs\.Recent work derives finite\-time convergence bounds for discretized ODE generation\(Benton et al\.,[2023](https://arxiv.org/html/2605.06680#bib.bib11)\), typically in terms of the Lipschitz constantLLor assumptions on score estimation error\. Our logarithmic norm analysis provides strictly tighter bounds \(μ\+≤L\\mu\_\{\+\}\\leq L, with equality only whenΩ=0\\Omega=0\) and reveals the structural decomposition underlying tightness\. The connection between OT displacement interpolation and improved Euler accuracy \(Theorem[5](https://arxiv.org/html/2605.06680#Thmtheorem5)\) appears to be new in this literature\.

## 9Discussion

Implications for model design\.Our analysis suggests that future work on few\-step generation could benefit from: \(i\) explicitly monitoring‖S‖F2\\\|S\\\|\_\{F\}^\{2\}as a diagnostic for integration stiffness \(sinceμ\+=λmax​\(S\)\\mu\_\{\+\}=\\lambda\_\{\\max\}\(S\)governs exponential error growth\); \(ii\) exploring gradient\-field parameterizations \(v=∇ϕv=\\nabla\\phi\) that enforce irrotationality by construction; \(iii\) designing time\-dependent regularization schedules informed by the Gaussian analysis of Section[6](https://arxiv.org/html/2605.06680#S6), which shows strain peaks neart=0t=0andt=1t=1\.

Practical guidance for choosingα\\alphaandβ\\beta\.We recommend targeting a regularization\-to\-loss ratio of Reg/ℒFM≈10​–​20%\\mathcal\{L\}\_\{\\text\{FM\}\}\\approx 10\\text\{\-\-\}20\\%\. Given the dimensional scaling‖∇v‖F2=O​\(d\)\\\|\\nabla v\\\|\_\{F\}^\{2\}=O\(d\), a useful starting point isα≈0\.15⋅ℒFM/‖∇v‖F2\\alpha\\approx 0\.15\\cdot\\mathcal\{L\}\_\{\\text\{FM\}\}/\\\|\\nabla v\\\|\_\{F\}^\{2\}, which can be estimated from a few training steps\. In our experiments, this heuristic yieldedα≈0\.3\\alpha\\approx 0\.3ford=2d=2andα≈10−6\\alpha\\approx 10^\{\-6\}ford=3072d=3072, both within the effective range\. Forβ\\beta, our default recommendation is to treat it as a secondary weight: start fromβ=0\\beta=0to test the strain\-only prediction, then compare against eitherβ=α\\beta=\\alpha\(equal\-budget Jacobian penalty\) or a smaller strain\-dominant choice such asβ∈\[0\.1​α,0\.5​α\]\\beta\\in\[0\.1\\alpha,\\,0\.5\\alpha\]\. This is exactly the logic behind our reported settings:β=0\\beta=0isolates the theoretically favored component,β=α\\beta=\\alphatests whether isotropic Jacobian penalization is competitive, and intermediate choices such as\(α,β\)=\(0\.3,0\.05\)\(\\alpha,\\beta\)=\(0\.3,0\.05\)probe whether modest vorticity suppression adds value without diluting strain control\.

Time\-dependent weighting\.The Gaussian analysis \(Section[6](https://arxiv.org/html/2605.06680#S6)\) shows that‖S​\(t\)‖F\\\|S\(t\)\\\|\_\{F\}peaks neart=0t=0andt=1t=1\. This suggests that time\-dependent weightsα​\(t\)\\alpha\(t\)concentrated at the endpoints could reduce strain where it matters most, while minimizing bias at intermediate times\. We did not experiment withα​\(t\)\\alpha\(t\)schedules in this work; this is a promising direction that could improve the bias\-complexity tradeoff, particularly for high\-NFE preservation\.

Higher\-order solvers\.Our analysis is specific to the Euler method\. For Heun’s method \(second\-order\), the local truncation error involves third derivatives rather than second, and the strain/vorticity decomposition of these higher\-order terms may yield a different asymmetry\. We conjecture that the qualitative conclusion \(strain matters more than vorticity\) persists, since the logarithmic norm governs error propagation regardless of truncation order, but the quantitative gap between regimes may narrow\. Empirical investigation with Heun/RK4 solvers is an important direction for future work\.

Relationship to training\-free sampling improvements\.Orthogonal to model\-side regularization, recent work improves sampling via solver\-side innovations: adaptive step\-size schedules, backward\-error\-informed scheduling, and curvature\-aware solvers\. These approaches reduce discretization error*without*modifying the velocity field\. Our analysis is complementary: it identifies which properties ofvvmake it amenable \(or resistant\) to efficient integration, regardless of the solver\. In principle, combining model\-side strain reduction with solver\-side adaptivity could yield compounding benefits\.

Statistical considerations\.FID scores are computed on 50K generated samples using a single random seed per configuration\. FID has inherent variance \(±0\.5​–​1\.0\\pm 0\.5\\text\{\-\-\}1\.0at our quality levels\), so differences smaller than∼2\\sim\\\!2points should be interpreted cautiously\. The key comparisons in Table[2](https://arxiv.org/html/2605.06680#S7.T2)\(e\.g\., control vs\.α\\alpha\-only: 26\.0 vs\. 22\.2 at NFE==10\) exceed this noise floor\. Future work should report confidence intervals via multiple seeds\.

Limitations\.Our theoretical analysis assumes Euler integration; extension to higher\-order solvers remains open \(see above\)\. The CIFAR\-10 experiments use a modest model \(27M parameters\); scaling to competitive baselines \(EDM2, DiT\) would strengthen the empirical evidence and enable direct comparison with few\-step methods such as Rectified Flow, Consistency Models, and MeanFlow under standardized protocols\. The regularizer operates on‖S‖F\\\|S\\\|\_\{F\}, while the amplification bound depends onλmax​\(S\)\\lambda\_\{\\max\}\(S\); as discussed in Section[5\.2](https://arxiv.org/html/2605.06680#S5.SS2), the Frobenius penalty is sufficient but not tight in high dimensions\. Theβ\\beta\-only result on CIFAR\-10 \(Table[2](https://arxiv.org/html/2605.06680#S7.T2)\) suggests an implicit regularization effect from the VJP computational pathway that is not captured by our theory and warrants further investigation\.

Broader impact\.This work provides analytical tools for understanding and improving generative model efficiency\. The theoretical insights \(strain/vorticity asymmetry, OT irrotationality\) are general and may find applications in other ODE\-based generative frameworks beyond flow matching\.

## 10Conclusion

We have shown that the symmetric and antisymmetric parts of the velocity Jacobian play fundamentally different roles in numerical integration error for flow matching: strain drives exponential error amplification while vorticity contributes only linearly\. This asymmetry, formalized through the logarithmic norm, helps explain why some velocity fields require many integration steps while others do not, and provides principled guidance for regularization\. We further showed that optimal transport flows are irrotational and have zero material derivative, yielding second\-order Euler accuracy in the Eulerian analysis; for exact displacement interpolation, the associated Lagrangian particle dynamics are integrated exactly by Euler, which explains the machine\-precision behavior observed on both Gaussian and nonlinear OT flows\. Experiments support the main predictions of the theory: on 2D benchmarks, strain regularization yields2\.7×2\.7\\timeserror reduction, while preliminary CIFAR\-10 experiments show consistent low\-NFE improvements under Jacobian regularization together with a matched fine\-tuning control\. We hope these analytical tools will complement engineering advances in few\-step generation and motivate further connections between numerical analysis, optimal transport, and generative modeling\.

## References

- Lipman et al\. \[2023\]Y\. Lipman, R\. T\. Q\. Chen, H\. Ben\-Hamu, M\. Nickel, and M\. Le\.Flow matching for generative modeling\.In*ICLR*, 2023\.
- Liu et al\. \[2023\]X\. Liu, C\. Gong, and Q\. Liu\.Flow straight and fast: Learning to generate and transfer data with rectified flow\.In*ICLR*, 2023\.
- Song et al\. \[2023\]Y\. Song, P\. Dhariwal, M\. Chen, and I\. Sutskever\.Consistency models\.In*ICML*, 2023\.
- Geng et al\. \[2025\]Z\. Geng, M\. Deng, X\. Bai, J\. Z\. Kolter, and K\. He\.Mean flows for one\-step generative modeling\.*arXiv:2505\.13447*, 2025\.
- Geng et al\. \[2025\]Z\. Geng, A\. Pokle, W\. Luo, J\. Lin, and J\. Z\. Kolter\.Consistency models made easy\.In*ICLR*, 2025\.
- Brenier \[1991\]Y\. Brenier\.Polar factorization and monotone rearrangement of vector\-valued functions\.*Comm\. Pure Appl\. Math\.*, 44\(4\):375–417, 1991\.
- Benamou & Brenier \[2000\]J\.\-D\. Benamou and Y\. Brenier\.A computational fluid mechanics solution to the Monge\-Kantorovich mass transfer problem\.*Numer\. Math\.*, 84:375–393, 2000\.
- Söderlind \[2006\]G\. Söderlind\.The logarithmic norm\. History and modern theory\.*BIT Numer\. Math\.*, 46\(3\):631–652, 2006\.
- Hoffman et al\. \[2019\]J\. Hoffman, D\. A\. Roberts, and S\. Yaida\.Robust learning with Jacobian regularization\.*arXiv:1908\.02729*, 2019\.
- De Marinis et al\. \[2025\]A\. De Marinis, N\. Guglielmi, S\. Sicilia, and F\. Tudisco\.Improving the robustness of neural ODEs with minimal weight perturbation\.*arXiv:2501\.10740*, 2025\.
- Benton et al\. \[2023\]J\. Benton, G\. Deligiannidis, and A\. Doucet\.Error bounds for flow matching methods\.*arXiv:2305\.16860*, 2023\.
- Villani \[2003\]C\. Villani\.*Topics in Optimal Transportation*\.AMS, 2003\.
- Eguchi \[2025\]S\. Eguchi\.Implicit geometric regularization in flow matching via density weighted Stein operators\.*arXiv:2512\.23956*, 2025\.

## Appendix AComplete Proofs: Separated Error Bound

### A\.1Frobenius Orthogonality

###### Proposition 9\.

ForSSsymmetric andΩ\\Omegaantisymmetric:tr​\(S⊤​Ω\)=0\\mathrm\{tr\}\(S^\{\\top\}\\Omega\)=0, hence‖S\+Ω‖F2=‖S‖F2\+‖Ω‖F2\\\|S\+\\Omega\\\|\_\{F\}^\{2\}=\\\|S\\\|\_\{F\}^\{2\}\+\\\|\\Omega\\\|\_\{F\}^\{2\}\.

###### Proof\.

tr​\(S​Ω\)=tr​\(\(S​Ω\)⊤\)=tr​\(Ω⊤​S\)=tr​\(−Ω​S\)=−tr​\(S​Ω\)\\mathrm\{tr\}\(S\\Omega\)=\\mathrm\{tr\}\(\(S\\Omega\)^\{\\top\}\)=\\mathrm\{tr\}\(\\Omega^\{\\top\}S\)=\\mathrm\{tr\}\(\-\\Omega S\)=\-\\mathrm\{tr\}\(S\\Omega\), where the last step uses the cyclic property of trace\. Hencetr​\(S​Ω\)=0\\mathrm\{tr\}\(S\\Omega\)=0\. ∎

### A\.2Logarithmic Norm

###### Proposition 10\.

μ2​\(A\)=λmax​\(\(A\+A⊤\)/2\)=λmax​\(SA\)\\mu\_\{2\}\(A\)=\\lambda\_\{\\max\}\(\(A\+A^\{\\top\}\)/2\)=\\lambda\_\{\\max\}\(S\_\{A\}\)\.

###### Proof\.

For the ODEy˙=A​y\\dot\{y\}=Ay:dd​t​‖y‖2=2​⟨y,A​y⟩=2​⟨y,S​y⟩\+2​⟨y,Ω​y⟩\\frac\{d\}\{dt\}\\\|y\\\|^\{2\}=2\\langle y,Ay\\rangle=2\\langle y,Sy\\rangle\+2\\langle y,\\Omega y\\rangle\. SinceΩ\\Omegais antisymmetric,⟨y,Ω​y⟩=0\\langle y,\\Omega y\\rangle=0\. Thereforedd​t​‖y‖2=2​⟨y,S​y⟩≤2​λmax​\(S\)​‖y‖2\\frac\{d\}\{dt\}\\\|y\\\|^\{2\}=2\\langle y,Sy\\rangle\\leq 2\\lambda\_\{\\max\}\(S\)\\\|y\\\|^\{2\}\. By Grönwall:‖y​\(t\)‖≤‖y​\(0\)‖​eλmax​\(S\)​t\\\|y\(t\)\\\|\\leq\\\|y\(0\)\\\|e^\{\\lambda\_\{\\max\}\(S\)t\}\. Since this holds for ally​\(0\)y\(0\):‖et​A‖≤et​λmax​\(S\)\\\|e^\{tA\}\\\|\\leq e^\{t\\lambda\_\{\\max\}\(S\)\}\. ∎

### A\.3Local Truncation Error Decomposition

The local truncation error at stepnnis:τn=h22​\[∂tv\+\(∇v\)​v\]\+O​\(h3\)=h22​\[∂tv\+S​v\+Ω​v\]\+O​\(h3\)\\tau\_\{n\}=\\frac\{h^\{2\}\}\{2\}\\left\[\\partial\_\{t\}v\+\(\\nabla v\)v\\right\]\+O\(h^\{3\}\)=\\frac\{h^\{2\}\}\{2\}\\left\[\\partial\_\{t\}v\+Sv\+\\Omega v\\right\]\+O\(h^\{3\}\)\. By triangle inequality:‖τn‖≤h22​\(Mt\+MS\+MΩ\)\+O​\(h3\)\\\|\\tau\_\{n\}\\\|\\leq\\frac\{h^\{2\}\}\{2\}\(M\_\{t\}\+M\_\{S\}\+M\_\{\\Omega\}\)\+O\(h^\{3\}\)\.

### A\.4Discrete Grönwall Lemma

###### Lemma 11\.

Ifan\+1≤\(1\+δ\)​an\+Ba\_\{n\+1\}\\leq\(1\+\\delta\)a\_\{n\}\+Bwitha0=0a\_\{0\}=0, thenan≤Bδ​\(\(1\+δ\)n−1\)a\_\{n\}\\leq\\frac\{B\}\{\\delta\}\(\(1\+\\delta\)^\{n\}\-1\)\.

###### Proof\.

By induction\. Base:a0=0a\_\{0\}=0\. Step:ak\+1≤\(1\+δ\)​\[Bδ​\(\(1\+δ\)k−1\)\]\+B=Bδ​\(\(1\+δ\)k\+1−1\)a\_\{k\+1\}\\leq\(1\+\\delta\)\[\\frac\{B\}\{\\delta\}\(\(1\+\\delta\)^\{k\}\-1\)\]\+B=\\frac\{B\}\{\\delta\}\(\(1\+\\delta\)^\{k\+1\}\-1\)\. ∎

### A\.5Proof of Theorem[1](https://arxiv.org/html/2605.06680#Thmtheorem1)

Taking norms:‖en\+1‖≤‖I\+h​∇v‖​‖en‖\+‖τn‖≤\(1\+h​μ\+\)​‖en‖\+h22​\(Mt\+MS\+MΩ\)\\\|e\_\{n\+1\}\\\|\\leq\\\|I\+h\\nabla v\\\|\\\|e\_\{n\}\\\|\+\\\|\\tau\_\{n\}\\\|\\leq\(1\+h\\mu\_\{\+\}\)\\\|e\_\{n\}\\\|\+\\frac\{h^\{2\}\}\{2\}\(M\_\{t\}\+M\_\{S\}\+M\_\{\\Omega\}\)\.

Apply Grönwall withδ=h​μ\+\\delta=h\\mu\_\{\+\},B=h22​\(Mt\+MS\+MΩ\)B=\\frac\{h^\{2\}\}\{2\}\(M\_\{t\}\+M\_\{S\}\+M\_\{\\Omega\}\),N=T/hN=T/h:‖eN‖≤h​\(Mt\+MS\+MΩ\)2​μ\+​\(eμ\+​T−1\)\+O​\(h2\)\\\|e\_\{N\}\\\|\\leq\\frac\{h\(M\_\{t\}\+M\_\{S\}\+M\_\{\\Omega\}\)\}\{2\\mu\_\{\+\}\}\(e^\{\\mu\_\{\+\}T\}\-1\)\+O\(h^\{2\}\)\.

## Appendix BComplete Proofs: Optimal Transport

### B\.1Proof of Theorem[4](https://arxiv.org/html/2605.06680#Thmtheorem4)

The displacement interpolation givesφt​\(x\)=\(1−t\)​x\+t​∇Ψ​\(x\)\\varphi\_\{t\}\(x\)=\(1\-t\)x\+t\\nabla\\Psi\(x\), with Jacobian∇xφt=\(1−t\)​I\+t​∇2Ψ\\nabla\_\{x\}\\varphi\_\{t\}=\(1\-t\)I\+t\\nabla^\{2\}\\Psi\. SettingH=∇2ΨH=\\nabla^\{2\}\\Psi\(symmetric p\.d\.\), the Eulerian velocity Jacobian is:

∇yvO​T=\(H−I\)​\[\(1−t\)​I\+t​H\]−1\\nabla\_\{y\}v^\{OT\}=\(H\-I\)\[\(1\-t\)I\+tH\]^\{\-1\}\.

SinceHHis symmetric, it has spectral decompositionH=Q​Λ​Q⊤H=Q\\Lambda Q^\{\\top\}\. Then:

\(H−I\)​\[\(1−t\)​I\+t​H\]−1=Q​\(Λ−I\)​\[\(1−t\)​I\+t​Λ\]−1​Q⊤=Q​diag​\(λi−1\(1−t\)\+t​λi\)​Q⊤\(H\-I\)\[\(1\-t\)I\+tH\]^\{\-1\}=Q\(\\Lambda\-I\)\[\(1\-t\)I\+t\\Lambda\]^\{\-1\}Q^\{\\top\}=Q\\,\\mathrm\{diag\}\\\!\\left\(\\frac\{\\lambda\_\{i\}\-1\}\{\(1\-t\)\+t\\lambda\_\{i\}\}\\right\)Q^\{\\top\},

which is symmetric\. HenceΩO​T=0\\Omega^\{OT\}=0\.∎

### B\.2Proof of Theorem[5](https://arxiv.org/html/2605.06680#Thmtheorem5)

In Lagrangian coordinates,φ˙t​\(x\)=∇Ψ​\(x\)−x\\dot\{\\varphi\}\_\{t\}\(x\)=\\nabla\\Psi\(x\)\-xis time\-independent\. Thereforeφ¨t​\(x\)=0\\ddot\{\\varphi\}\_\{t\}\(x\)=0\. The material derivative equals the Lagrangian acceleration:D​vD​t\|\(t,φt​\(x\)\)=φ¨t​\(x\)=0\\frac\{Dv\}\{Dt\}\\big\|\_\{\(t,\\varphi\_\{t\}\(x\)\)\}=\\ddot\{\\varphi\}\_\{t\}\(x\)=0\. Since theO​\(h2\)O\(h^\{2\}\)truncation error term ish22​D​vD​t=0\\frac\{h^\{2\}\}\{2\}\\frac\{Dv\}\{Dt\}=0, the leading error isO​\(h3\)O\(h^\{3\}\), yielding second\-order global convergence\.∎