Bi-Lipschitz Normalizing Flows Explained

• Bi-Lipschitz normalizing flows are defined as invertible maps with both forward and inverse Lipschitz bounds, ensuring controlled local volume distortion.
• They stabilize likelihood optimization by bounding Jacobian determinants, though such constraints can limit expressivity when modeling sharp or separated distributions.
• Architectural realizations—including residual, monotone-operator, and spline-based flows—balance stability with model capacity, addressing practical trade-offs in generative modeling.

Searching arXiv for recent and foundational papers on bi-Lipschitz normalizing flows. arXiv query: bi-Lipschitz normalizing flows expressivity monotone operators B-spline flows Bi-Lipschitz normalizing flows are normalizing flows whose forward map and inverse are both globally Lipschitz with finite bounds. In the standard formulation, an invertible map $$F : X \subset \mathbb{R}^d \to Z \subset \mathbb{R}^d$$ is $(L_1,L_2)$-bi-Lipschitz if

$$\forall x,y \in X:\ \|F(x)-F(y)\| \le L_1\|x-y\|, \qquad \forall u,v \in Z:\ \|F^{-1}(u)-F^{-1}(v)\| \le L_2\|u-v\|,$$

equivalently,

$$\frac{1}{L_2}\|x-y\| \le \|F(x)-F(y)\| \le L_1\|x-y\|.$$

With a base density $q$ and differentiable bijection $F$, the induced density is $\hat p(x)=q(F(x))\,|\det D F(x)|$ (Verine et al., 2021). This class is central because many modern flow architectures are either bi-Lipschitz by design or trained under Lipschitz control to ensure numerical stability, tractable inverses, and robust optimization, yet the same regularity imposes nontrivial expressivity constraints (Verine et al., 2021). Subsequent work has refined both sides of this picture: negative results under uniform Lipschitz constraints, positive $L^1$-density results when Lipschitz bounds are allowed to depend on the approximation target, and architectural constructions that preserve explicit bi-Lipschitz control while enlarging the class of realizable transports (Iske et al., 7 May 2026).

1. Definitions, geometric meaning, and likelihood structure

For differentiable bi-Lipschitz maps, the Jacobian is uniformly controlled. If $F$ is differentiable and $L_1$-Lipschitz, then $(L_1,L_2)$0 for all $(L_1,L_2)$1, while $(L_1,L_2)$2 being $(L_1,L_2)$3-Lipschitz implies $(L_1,L_2)$4. Hence

$(L_1,L_2)$5

for all $(L_1,L_2)$6 (Verine et al., 2021). Equivalent statements appear in broader flow-theoretic treatments: if a map $(L_1,L_2)$7 is bi-Lipschitz with constants $(L_1,L_2)$8, then singular values of $(L_1,L_2)$9 lie in $$\forall x,y \in X:\ \|F(x)-F(y)\| \le L_1\|x-y\|, \qquad \forall u,v \in Z:\ \|F^{-1}(u)-F^{-1}(v)\| \le L_2\|u-v\|,$$0 almost everywhere and

$$\forall x,y \in X:\ \|F(x)-F(y)\| \le L_1\|x-y\|, \qquad \forall u,v \in Z:\ \|F^{-1}(u)-F^{-1}(v)\| \le L_2\|u-v\|,$$1

(Koehler et al., 2020, Liao et al., 2021, Hong et al., 2023).

These determinant bounds encode bounded local volume distortion. Regions cannot be expanded arbitrarily and cannot be collapsed arbitrarily. In generative modeling terms, the pushforward density cannot exhibit unlimited dynamic range relative to the base distribution because both the latent density term and the Jacobian term are uniformly constrained (Verine et al., 2021). This is the geometric core of why bi-Lipschitzness stabilizes inversion and likelihood optimization while also restricting approximation of sharply concentrated or highly separated targets.

The likelihood objective follows the usual change-of-variables formula,

$$\forall x,y \in X:\ \|F(x)-F(y)\| \le L_1\|x-y\|, \qquad \forall u,v \in Z:\ \|F^{-1}(u)-F^{-1}(v)\| \le L_2\|u-v\|,$$2

with $$\forall x,y \in X:\ \|F(x)-F(y)\| \le L_1\|x-y\|, \qquad \forall u,v \in Z:\ \|F^{-1}(u)-F^{-1}(v)\| \le L_2\|u-v\|,$$3 the transport to latent space (Koehler et al., 2020, Liao et al., 2021). In a bi-Lipschitz flow, the log-determinant term is automatically bounded: $$\forall x,y \in X:\ \|F(x)-F(y)\| \le L_1\|x-y\|, \qquad \forall u,v \in Z:\ \|F^{-1}(u)-F^{-1}(v)\| \le L_2\|u-v\|,$$4 (Koehler et al., 2020, Liao et al., 2021). This mitigates numerical instability during maximum-likelihood training. A related theoretical result states that for given base and target distributions, the Jacobian determinant mapping is unique almost everywhere, and for the Quasi-Linear Flow class the globally optimal value of the log-likelihood can be written explicitly in terms of eigenvalues of the pointwise auto-correlation matrix (Liao et al., 2021). That result concerns likelihood structure rather than expressivity per se, but it sharpens the role of determinant control in flow design.

Total variation distance is a preferred metric in the expressivity analysis of uniformly bi-Lipschitz flows. For distributions $$\forall x,y \in X:\ \|F(x)-F(y)\| \le L_1\|x-y\|, \qquad \forall u,v \in Z:\ \|F^{-1}(u)-F^{-1}(v)\| \le L_2\|u-v\|,$$5 with densities $$\forall x,y \in X:\ \|F(x)-F(y)\| \le L_1\|x-y\|, \qquad \forall u,v \in Z:\ \|F^{-1}(u)-F^{-1}(v)\| \le L_2\|u-v\|,$$6,

$$\forall x,y \in X:\ \|F(x)-F(y)\| \le L_1\|x-y\|, \qquad \forall u,v \in Z:\ \|F^{-1}(u)-F^{-1}(v)\| \le L_2\|u-v\|,$$7

Pinsker’s inequality gives $$\forall x,y \in X:\ \|F(x)-F(y)\| \le L_1\|x-y\|, \qquad \forall u,v \in Z:\ \|F^{-1}(u)-F^{-1}(v)\| \le L_2\|u-v\|,$$8, so lower bounds in total variation immediately imply nontrivial Kullback–Leibler lower bounds (Verine et al., 2021). The same paper notes that precision/recall style support-based analyses become trivial when both supports are full-dimensional, which motivates TV-based lower bounds.

2. Why bi-Lipschitzness is pervasive in flow architectures

Many widely used architectures enforce or approximately enforce bi-Lipschitz behavior. In invertible residual networks and Residual Flows, blocks are of the form $$\forall x,y \in X:\ \|F(x)-F(y)\| \le L_1\|x-y\|, \qquad \forall u,v \in Z:\ \|F^{-1}(u)-F^{-1}(v)\| \le L_2\|u-v\|,$$9 with $$\frac{1}{L_2}\|x-y\| \le \|F(x)-F(y)\| \le L_1\|x-y\|.$$0, which implies for $$\frac{1}{L_2}\|x-y\| \le \|F(x)-F(y)\| \le L_1\|x-y\|.$$1 blocks

$$\frac{1}{L_2}\|x-y\| \le \|F(x)-F(y)\| \le L_1\|x-y\|.$$2

(Verine et al., 2021). The same scaling appears in practical parameterizations of iResNets, where residual contraction is enforced by spectral normalization (Iske et al., 7 May 2026).

Glow-like models bound forward and inverse Lipschitz constants through products of operator norms of invertible linear transforms and their inverses (Verine et al., 2021). In affine coupling layers, the Jacobian is triangular, and

$$\frac{1}{L_2}\|x-y\| \le \|F(x)-F(y)\| \le L_1\|x-y\|.$$3

so bounded scale outputs keep the layer bi-Lipschitz and the log-determinant bounded (Koehler et al., 2020). One practical prescription is to bound $$\frac{1}{L_2}\|x-y\| \le \|F(x)-F(y)\| \le L_1\|x-y\|.$$4 so that $$\frac{1}{L_2}\|x-y\| \le \|F(x)-F(y)\| \le L_1\|x-y\|.$$5, thereby keeping singular values of the triangular Jacobian in the same interval (Koehler et al., 2020). Spectral normalization, weight clipping, and related operator-norm constraints are standard tools for controlling these constants during training (Verine et al., 2021, Iske et al., 7 May 2026).

This architectural ubiquity arises from a basic trade-off. Small per-layer Lipschitz constants stabilize inverse evaluation, control Jacobian singular values, and prevent exploding determinants, but they also restrict how much a model can locally expand or compress probability mass (Verine et al., 2021). The literature repeatedly frames this as a conditioning–expressivity tension. The depth and conditioning analysis of affine couplings shows that shallow universality in Wasserstein distance is possible only when ill-conditioning is allowed; as approximation error is driven toward zero in the constructive result, the condition number blows up at least as fast as $$\frac{1}{L_2}\|x-y\| \le \|F(x)-F(y)\| \le L_1\|x-y\|.$$6, so strict bi-Lipschitz control is violated (Koehler et al., 2020). Conversely, under bounded Lipschitzness and limited layer width, depth lower bounds appear for approximating complex distributions (Koehler et al., 2020).

A notable architectural response is the monotone-operator formulation. “Monotone Flows” parameterize a 1-Lipschitz Cayley operator $$\frac{1}{L_2}\|x-y\| \le \|F(x)-F(y)\| \le L_1\|x-y\|.$$7 with $$\frac{1}{L_2}\|x-y\| \le \|F(x)-F(y)\| \le L_1\|x-y\|.$$8 and define

$$\frac{1}{L_2}\|x-y\| \le \|F(x)-F(y)\| \le L_1\|x-y\|.$$9

The induced block is globally invertible, $q$0-strongly monotone, and $q$1-Lipschitz with

$q$2

so $q$3 is bi-Lipschitz with explicit forward and inverse bounds (Ahn et al., 2022). This construction is designed precisely to enlarge expressivity relative to contraction-based residual flows while keeping explicit global conditioning control.

3. Expressivity limits under uniform bi-Lipschitz constraints

The most direct negative theory studies approximation from a standard Gaussian base under a uniformly bi-Lipschitz transport. Several lower bounds in total variation identify target distributions that are intrinsically hard to approximate when $q$4 and $q$5 are fixed (Verine et al., 2021).

A general dense-subset bound states that if $q$6 is $q$7-Lipschitz and $q$8 is the average target density on a measurable set $q$9, then

$F$0

Thus if some region has average density exceeding $F$1, the total variation distance cannot vanish (Verine et al., 2021). This formalizes the intuition that a uniformly $F$2-Lipschitz map cannot create arbitrarily dense regions from a standard Gaussian base.

A more geometric dense-ball version uses the fact that an $F$3-Lipschitz map sends a ball $F$4 into a latent ball of radius at most $F$5. The resulting tight closed-form lower bound is

$F$6

where $F$7 is the lower incomplete gamma function (Verine et al., 2021). A looser but dimension-uniform version replaces the Gaussian mass term with $F$8.

The complementary sparse-ball bound uses the inverse Lipschitz constant. If $F$9 is $\hat p(x)=q(F(x))\,|\det D F(x)|$0-Lipschitz, then around the special point $\hat p(x)=q(F(x))\,|\det D F(x)|$1,

$\hat p(x)=q(F(x))\,|\det D F(x)|$2

This states that the model must allocate at least a certain amount of mass near the inverse image of the latent mode, because the inverse cannot compress too aggressively there (Verine et al., 2021).

These two mechanisms generate a catalogue of unfavorable target distributions. High-mass spikes, thin ridges, tightly clustered modes, or mass concentrated near low-dimensional manifolds activate the forward-Lipschitz obstruction. Well-separated multimodal mixtures with low-density valleys activate the inverse-Lipschitz obstruction, because the Gaussian base places maximal density near $\hat p(x)=q(F(x))\,|\det D F(x)|$3 while the target may place little mass near $\hat p(x)=q(F(x))\,|\det D F(x)|$4 (Verine et al., 2021). A corollary for two separated target regions $\hat p(x)=q(F(x))\,|\det D F(x)|$5 at distance at least $\hat p(x)=q(F(x))\,|\det D F(x)|$6 gives

$\hat p(x)=q(F(x))\,|\det D F(x)|$7

together with an explicit upper bound on maximum precision,

$\hat p(x)=q(F(x))\,|\det D F(x)|$8

(Verine et al., 2021).

Dimension modifies the force of these bounds. For the dense-ball lower bound, the Gaussian mass term $\hat p(x)=q(F(x))\,|\det D F(x)|$9 decays rapidly in high dimension unless $L^1$0 is comparable to the typical Gaussian radius, which the paper describes as approximately $L^1$1 (Verine et al., 2021). This means that matching a fixed dense ball generally requires larger $L^1$2 as dimension grows. By contrast, for the sparse-ball bound, the relevant Gaussian mass decreases with dimension when the effective radius $L^1$3 is small, so the inverse-Lipschitz limitation is more salient in low dimension for fixed radius ratio (Verine et al., 2021).

Empirical examples in the same work illustrate the practical content of these lower bounds. On CIFAR-10, for small per-block $L^1$4 such as $L^1$5, at least about 10 residual blocks are needed to drive the dense-ball lower bound near zero; for a one-dimensional Gaussian with small variance, at least 7 layers with per-block $L^1$6 are suggested; on 2D Circles, the inner circle can force the lower bound toward $L^1$7; and for separated multimodal targets such as the two-cluster setup and 8-Gaussians, small per-block $L^1$8 may require at least 15 layers to overcome the sparse-ball lower bound in practice (Verine et al., 2021).

4. Universal approximation from a diffusion perspective

A later line of work argues that the preceding negative results do not show that bi-Lipschitzness itself is an intrinsic barrier. Rather, the obstruction is the imposition of a uniform Lipschitz budget independent of the target accuracy (Iske et al., 7 May 2026). The key construction uses the probability flow ODE associated with variance-preserving diffusion.

For the constant-rate variance-preserving process,

$L^1$9

the associated probability flow ODE is

$F$0

(Iske et al., 7 May 2026). If the score satisfies

$F$1

for all $F$2 and $F$3 for every $F$4, then the solution map $F$5 is a bi-Lipschitz $F$6-diffeomorphism for all $F$7 (Iske et al., 7 May 2026). Moreover,

$F$8

and

$F$9

(Iske et al., 7 May 2026).

This gives a deterministic ODE bridge from Gaussian to target whose transport maps are bi-Lipschitz whenever the score is regular enough. The induced Lipschitz constants are explicit through Grönwall: $L_1$0 with the same bound for the inverse map (Iske et al., 7 May 2026).

The same paper verifies score regularity for several classes of targets. For compactly supported densities, explicit bounds on $L_1$1 follow from a covariance representation of the Hessian of $L_1$2 (Iske et al., 7 May 2026). For log-concave densities, Brascamp–Lieb yields

$L_1$3

on any interval $L_1$4 (Iske et al., 7 May 2026). For Gaussian convolutions of compactly supported measures, and therefore for finite Gaussian mixtures, the score gradient is uniformly bounded on every $L_1$5, which in turn yields KL convergence and $L_1$6 convergence without early stopping (Iske et al., 7 May 2026).

The culminating statement is a universal approximation theorem: Gaussian pullbacks induced by bi-Lipschitz variance-preserving transport maps are $L_1$7-dense among all probability densities. Formally, for every density $L_1$8 and every $L_1$9, there exists a bi-Lipschitz VP flow $(L_1,L_2)$00 such that

$(L_1,L_2)$01

(Iske et al., 7 May 2026). The paper’s interpretation is explicit: bi-Lipschitz constraints per se do not destroy universal approximation power; the negative results arise from demanding uniform bounds that do not scale with the target or the desired error. A plausible implication is that the apparent contradiction between lower bounds and universality is resolved by distinguishing fixed-constant bi-Lipschitz classes from non-uniform unions of such classes.

5. Architectural realizations and design strategies

Several concrete architectures instantiate bi-Lipschitz normalizing flows with different trade-offs in smoothness, tractability, and expressivity.

A first family is based on affine couplings and invertible residual networks. Affine coupling layers preserve triangular Jacobians and analytic log-determinants; the analysis of depth and conditioning shows that constant-depth affine couplings can exactly represent permutations and invertible $(L_1,L_2)$02 convolutions, including the linear maps used in Glow, with $(L_1,L_2)$03 alternating triangular factors in the theorem statement (Koehler et al., 2020). At the same time, at least two full coupling layers are necessary for certain linear maps, and depth lower bounds remain under bounded Lipschitz constraints (Koehler et al., 2020). These results imply that partition choice is not the fundamental representational bottleneck, whereas conditioning control still exacts a depth cost.

A second family is based on monotone operators. In Monotone Flows, each block is built from the monotone formulation

$(L_1,L_2)$04

with inverse

$(L_1,L_2)$05

and the implicit relation

$(L_1,L_2)$06

(Ahn et al., 2022). The Jacobian has the closed form

$(L_1,L_2)$07

and the log-determinant admits the identity

$(L_1,L_2)$08

(Ahn et al., 2022). Because $(L_1,L_2)$09, both forward and inverse global Lipschitz constants are known explicitly. The paper combines spectral normalization with the 1-Lipschitz CPila activation,

$(L_1,L_2)$10

with $(L_1,L_2)$11, $(L_1,L_2)$12, and $(L_1,L_2)$13 in experiments (Ahn et al., 2022).

A third family emphasizes higher smoothness. Neural diffeomorphic non-uniform B-spline flows construct at least twice continuously differentiable coupling transforms while remaining bi-Lipschitz and retaining analytic inverses for the cubic case (Hong et al., 2023). For a non-uniform B-spline transform of order $(L_1,L_2)$14, the paper gives the sufficient condition

$(L_1,L_2)$15

which implies

$(L_1,L_2)$16

on the interval and hence yields a one-dimensional bi-Lipschitz map with explicit forward constant $(L_1,L_2)$17 and inverse constant $(L_1,L_2)$18 (Hong et al., 2023). In the coupling setting, the Jacobian determinant is the product of scalar derivatives,

$(L_1,L_2)$19

and second derivatives are available in closed form, which is important in force-matching and Boltzmann-generator applications (Hong et al., 2023).

The architectural landscape therefore spans triangular flows with bounded scales, residual-style flows with contraction or monotone-operator parameterizations, and spline-based diffeomorphisms with explicit derivative control. Across these variants, the recurring design principle is the same: keep Jacobian singular values bounded away from both $(L_1,L_2)$20 and $(L_1,L_2)$21, but decide how much expressive power is spent on exact likelihoods, analytic inverses, smoothness class, or numerical stability.

6. Trade-offs, empirical behavior, and unresolved questions

The strongest practical trade-offs concern base distribution choice, depth, and the strength of Lipschitz control. Changing the Gaussian base variance from $(L_1,L_2)$22 to $(L_1,L_2)$23 rescales the two principal lower bounds in opposite directions: $(L_1,L_2)$24

$(L_1,L_2)$25

Increasing $(L_1,L_2)$26 helps the dense-ball obstruction but worsens the sparse-ball obstruction, and decreasing $(L_1,L_2)$27 does the opposite (Verine et al., 2021). The same work therefore argues that changing variance alone does not generally remove expressivity limits.

Using more complex base distributions can mitigate these constraints. With a $(L_1,L_2)$28-component Gaussian mixture base, a favorable version of the dense-ball bound acquires a $(L_1,L_2)$29 factor, and mixture components can be placed near low-density valleys to reduce forced mass near $(L_1,L_2)$30 (Verine et al., 2021). The stated caveat is that choosing $(L_1,L_2)$31 and training such mixtures is challenging; mixture-of-flows and partition-based variants can alleviate limitations but often sacrifice simplicity, exact likelihoods, or training stability (Verine et al., 2021).

Empirically, monotone-operator parameterizations show that one can relax the restrictive $(L_1,L_2)$32 regime of contraction-based residual flows without surrendering global bi-Lipschitz guarantees. On density-estimation benchmarks, Monotone Flows report 0.928 bpd on MNIST, 3.215 bpd on CIFAR-10, 3.961 bpd on ImageNet32, and 3.734 bpd on ImageNet64 under uniform dequantization, with additional ablations isolating gains from the monotone formulation and the CPila activation (Ahn et al., 2022). These are concrete performance claims about one architecture, not generic properties of all bi-Lipschitz flows.

For applications requiring $(L_1,L_2)$33 structure, non-uniform cubic B-spline flows report runtime per sample for inversion of 1.12 ms versus 19.8 ms for smooth normalizing flows and 0.59 ms for RQ-spline flows, while in a 60-D alanine dipeptide force-matching task the non-uniform B-spline with force matching achieved approximately $(L_1,L_2)$34 FME versus approximately $(L_1,L_2)$35 for the smooth baseline (Hong et al., 2023). The stated interpretation is that analytic inversion and closed-form second derivatives matter when forces depend on Hessians.

Several open questions remain explicit in the literature. One is computational: verifying density or sparsity conditions behind the TV lower bounds in high dimension requires approximating suprema over all balls, which is costly (Verine et al., 2021). Another is statistical and architectural: translating functional existence results from the VP probability flow ODE into explicit, architecture-aware training prescriptions with provable end-to-end performance remains open (Iske et al., 7 May 2026). Additional open directions include extending lower-bound analyses to non-invertible or locally bi-Lipschitz architectures and to other metrics such as Wasserstein distance (Verine et al., 2021), understanding when a learned score can be faithfully distilled into an explicit invertible neural network (Iske et al., 7 May 2026), and deriving scalable regularizers that directly control inverse conditioning or minimum singular values during training (Liao et al., 2021).

Taken together, the modern theory of bi-Lipschitz normalizing flows is not reducible to a single verdict about “expressive” or “inexpressive.” Under fixed global bounds, there are explicit approximation gaps for dense spikes, deep valleys, separated modes, and related geometries (Verine et al., 2021). Under non-uniform bounds induced by score-regular probability flow ODEs, Gaussian pullbacks by bi-Lipschitz transports are $(L_1,L_2)$36-dense among all densities (Iske et al., 7 May 2026). Architectural work correspondingly focuses not on escaping bi-Lipschitzness altogether, but on deciding how those bounds are realized, aggregated, and traded against depth, smoothness, analytic tractability, and training stability (Koehler et al., 2020, Ahn et al., 2022, Hong et al., 2023).