Normalizing Flows: Invertible Transformations for Generative Modeling

Suppose you are training a model on 50,000 images of handwritten Devanagari characters, and you want it to do two things at once: (1) generate a brand-new, plausible-looking character it has never seen, and (2) given any image — real or fake — tell you exactly how likely that image was under the data it learned from, as an actual number, not a rough score. A GAN can do the first job well but has no honest way to do the second; it never learns a density, only a generator. A VAE gives you an approximate number for the second job, bounded below by something called the ELBO, but the approximation gap is real and usually unmeasured. Neither model can answer "what is the exact probability density of this exact image?"

Normalizing flows are the one major class of generative model that does both jobs exactly, with the same set of learned parameters, using nothing more exotic than a change of variables you can derive yourself from the definition of a cumulative distribution function. That derivation — and the clever architectural trick that makes it computationally survivable in high dimensions — is the entire content of this chapter.

The Core Idea in One Dimension

Start with the simplest possible version of the problem. Let X be a random variable with a known density p_X, and let Y = f(X) for some function f that is strictly increasing (so it is invertible) and differentiable. What is the density of Y?

You are not allowed to guess "just substitute x = f⁻¹(y) into p_X" — that is wrong, and seeing exactly why it is wrong is the whole point. Work from the definition of a CDF instead:

F_Y(y) = P(Y <= y)
       = P(f(X) <= y)                since Y = f(X)
       = P(X <= f^{-1}(y))           f is increasing, so f(X)<=y  <=>  X<=f^{-1}(y)
       = F_X(f^{-1}(y))

Differentiate both sides with respect to y, using the chain rule on the right:

p_Y(y) = d/dy F_Y(y) = p_X(f^{-1}(y)) * (f^{-1})'(y)

If f were decreasing instead of increasing, the same steps give a minus sign, which is fixed by taking an absolute value. The general 1-D change-of-variables formula is:

p_Y(y) = p_X(x) * |dx/dy|,   where x = f^{-1}(y)

The factor |dx/dy| is not decoration — it is the entire reason this formula is not "just plug in." Stretching or compressing the input axis changes how much probability mass gets packed into a unit interval of the output axis, and the density has to compensate so that total probability still integrates to exactly 1.

Concrete check: let X be Uniform on [0, 1], so p_X(x) = 1 on that interval, and let Y = f(X) = X², which is increasing on [0,1] and has inverse x = f⁻¹(y) = √y. Then dx/dy = 1/(2√y), so:

p_Y(y) = 1 * 1/(2*sqrt(y)) = 1/(2*sqrt(y)),   0 < y < 1

This blows up near y = 0 — which makes sense, because squaring crushes the interval near 0 into an even smaller interval, so probability mass piles up there. Verify it is a genuine density by integrating: ∫₀¹ 1/(2√y) dy = [√y]₀¹ = 1. It checks out. Notice that nowhere did we need to know anything about neural networks — this is pure single-variable calculus, and it is the seed from which the entire theory of normalizing flows grows.

From One Variable to Many: the Jacobian Determinant

Now let X and Y be vectors in ℝⁿ, related by an invertible, differentiable map Y = f(X). The single derivative dx/dy generalizes to the Jacobian matrix of partial derivatives, and the "how much does a small volume get stretched" role is played by the absolute value of its determinant:

p_Y(y) = p_X(x) * |det J_{f^{-1}}(y)|,   where x = f^{-1}(y)

Equivalently, working in the forward direction (this is the form you will use constantly, so fix it in memory):

p_X(x) = p_Y(f(x)) * |det J_f(x)|

where J_f(x) is the n×n matrix with entry (i,j) = ∂f_i/∂x_j. This is the multivariable change-of-variables formula from calculus — the same one that lets you convert a double integral from Cartesian to polar coordinates using the factor r, which is exactly |det J| for the polar-to-Cartesian map. Nothing new is being invented here; it is being repurposed.

A normalizing flow builds a generative model out of this formula directly. Pick a simple base distribution p_Z — almost always a standard multivariate Gaussian, because its density is trivial to write down and to sample from. Then learn an invertible, differentiable map f (built from a neural network) such that Z = f(X) pushes your complicated real data distribution toward that simple Gaussian. The name is literal: f "normalizes" the data distribution into a standard one, and a deep f is built by chaining — letting probability mass "flow" through — several simpler invertible layers, f = f_L ∘ f_{L-1} ∘ ... ∘ f_1. Once trained:

• Density evaluation (exact, not approximate): run data x forward through f to get z, evaluate the easy Gaussian density at z, and multiply by the accumulated |det J| of every layer.
• Sampling / generation: draw z from the Gaussian, run it backward through f⁻¹ to produce a new x.

The same weights do both jobs. This is strictly more than a GAN or a VAE offers, and it is why flows matter enough to have their own architecture family (NICE, RealNVP, Glow) rather than being a historical footnote.

The Determinant Is the Bottleneck

There is an obstacle standing between this formula and a working model: computing det J for an arbitrary n×n matrix costs O(n³) operations (Gaussian elimination / LU decomposition), and computing the full Jacobian of a generic neural network with n input and n output units in the first place is already expensive. For an image with just 28×28 = 784 pixels, n³ is on the order of 5×10⁸ — per layer, per training step. That is not survivable at scale. A generic invertible neural network is therefore useless as a flow layer even though it satisfies the math; you also need its Jacobian determinant to be cheap.

The fix is architectural, not mathematical: design f so that its Jacobian is triangular by construction. The determinant of a triangular matrix is just the product of its diagonal entries — an O(n) computation, however dense and nonlinear the rest of the matrix is. This single design constraint is the entire reason coupling layers, the workhorse of modern flows, look the way they do.

The Coupling-Layer Trick

Split the input vector into two pieces, x = (x₁, x₂). Define the output as:

y1 = x1
y2 = x2 * exp(s(x1)) + t(x1)

where s and t are arbitrary functions of x1 — in practice, ordinary feed-forward neural networks with no constraints on their weights at all. This is a deliberately asymmetric design: the first half passes through untouched, and it is used to compute a per-element scale (exp(s(x1))) and shift (t(x1)) applied to the second half.

Invertibility. Given (y1, y2), recovery is immediate: x1 = y1, and then since x1 is now known, s(x1) and t(x1) are just numbers you can compute, so x2 = (y2 - t(x1)) * exp(-s(x1)). No matrix inversion, no iterative solving — just plugging numbers back in. This works regardless of how complicated s and t are, because x2 only ever appears through a scale-and-shift (an affine map in x2), and an affine map with a guaranteed-nonzero, positive scale — exp(·) is never zero or negative — is always invertible. s and t themselves never need to be inverted at all.

The Jacobian. Treat x1, x2 as scalars first, to see the mechanics cleanly, before generalizing to blocks. Differentiate each output with respect to each input:

∂y1/∂x1 = 1
∂y1/∂x2 = 0

∂y2/∂x1 = x2 * exp(s(x1)) * s'(x1) + t'(x1)     (product rule on the first term,
                                                   plus the derivative of t(x1))
∂y2/∂x2 = exp(s(x1))

The ∂y2/∂x1 entry needs both pieces: differentiating x2·exp(s(x1)) with respect to x1 by the product rule gives x2·exp(s(x1))·s'(x1) (chain rule inside, since the exponent itself depends on x1), and differentiating the additive t(x1) term contributes its own derivative t'(x1) on top, by the sum rule. Assembled as a matrix:

       [ ∂y1/∂x1   ∂y1/∂x2 ]   [        1                    0        ]
J_f =  [                    ] = [                                      ]
       [ ∂y2/∂x1   ∂y2/∂x2 ]   [ x2*exp(s(x1))*s'(x1)+t'(x1)  exp(s(x1)) ]

Because the top-right entry is exactly 0, this matrix is lower triangular, and the determinant of a triangular matrix is the product of its diagonal entries only — every off-diagonal entry, however messy, is irrelevant to the determinant:

det(J_f) = (1) * exp(s(x1)) = exp(s(x1))

The bottom-left entry — the one with the product rule and the t'(x1) term — never appears in this product. This is the payoff of the triangular design: s and t can be ten-layer neural networks with millions of parameters, and the Jacobian determinant is still exp(s(x1)), computable in constant time, because the only quantity that matters is the diagonal.

Generalizing from scalars to blocks — x1 ∈ ℝ^{d1}, x2 ∈ ℝ^{d2}, with s, t : ℝ^{d1} → ℝ^{d2} — the Jacobian becomes block lower-triangular:

       [    I_{d1}              0        ]
J_f =  [                                  ]
       [  ∂y2/∂x1     diag(exp(s(x1)))    ]

and the determinant of a block-triangular matrix is the product of the determinants of its diagonal blocks:

det(J_f) = det(I_{d1}) * det(diag(exp(s(x1)))) = ∏_{i=1}^{d2} exp(s(x1)_i) = exp(Σ s(x1))

A single number: the exponential of the sum of the outputs of the scale network. That is the whole trick that makes normalizing flows tractable on real, high-dimensional data.

A Fully Worked Numerical Example

To make this concrete rather than symbolic, trace one coupling layer through actual numbers, using simple linear stand-ins s(x1) = 0.5·x1 and t(x1) = 2·x1 (in a real flow these would be small neural networks; linear functions are used here purely so the arithmetic can be checked by hand).

Take x1 = 1, x2 = 2.

s(x1) = 0.5 * 1 = 0.5        exp(0.5) ≈ 1.6487
t(x1) = 2 * 1 = 2

y1 = x1 = 1
y2 = x2 * exp(s(x1)) + t(x1) = 2 * 1.6487 + 2 ≈ 5.2974

det J = exp(s(x1)) = exp(0.5) ≈ 1.6487
log det J = s(x1) = 0.5

Check the inverse recovers the original input exactly:

x1 = y1 = 1
x2 = (y2 - t(x1)) * exp(-s(x1)) = (5.2974 - 2) * exp(-0.5) ≈ 3.2974 * 0.6065 ≈ 2.0000  ✓

Now verify this in code, and use it to compute an actual exact log-density under a standard Gaussian base distribution — the quantity a flow is trained to maximize:

import numpy as np

def s(x1): return 0.5 * x1
def t(x1): return 2.0 * x1

def coupling_forward(x1, x2):
    y1 = x1
    y2 = x2 * np.exp(s(x1)) + t(x1)
    log_det = s(x1)              # log(exp(s(x1))) = s(x1)
    return y1, y2, log_det

def coupling_inverse(y1, y2):
    x1 = y1
    x2 = (y2 - t(x1)) * np.exp(-s(x1))
    return x1, x2

x1, x2 = 1.0, 2.0
y1, y2, log_det = coupling_forward(x1, x2)
print(y1, y2, log_det)          # 1.0 5.297442541400256 0.5

# Round trip check
rx1, rx2 = coupling_inverse(y1, y2)
print(rx1, rx2)                 # 1.0 2.0000000000000004  (matches input)

# Exact log-density of x under base N(0, I), using change of variables:
# log p_X(x) = log p_Z(y) + log_det
def log_std_normal(v):
    return -0.5 * np.log(2 * np.pi) - 0.5 * v**2

log_pZ = log_std_normal(y1) + log_std_normal(y2)   # independent dims -> sum
log_pX = log_pZ + log_det
print(log_pX)                   # about -15.87

Every line above is checkable by hand from the formulas already derived: log_det = s(x1) = 0.5 exactly reproduces the triangular-determinant result, and log_pX is the change-of-variables formula from the "From One Variable to Many" section applied with an actual number coming out the other end — not a symbolic promise.

Stacking Layers: Why You Must Alternate

A single coupling layer has a glaring weakness: y1 = x1 is left completely untouched. If you stacked ten coupling layers all splitting the vector the same way, the first half of every output would still equal the first half of the original input, and the model could never reshape that half of the distribution at all. The fix is to alternate which half passes through unchanged — swap the two halves (a "flip" or fixed permutation) between successive coupling layers, so that a dimension left untouched by layer k gets transformed by layer k+1:

def flip(x1, x2):
    return x2, x1   # swap which half is "held fixed" for the next layer

def flow_forward(x1, x2, layers):
    total_log_det = 0.0
    for i in range(layers):
        x1, x2, ld = coupling_forward(x1, x2)
        total_log_det += ld
        x1, x2 = flip(x1, x2)
    return x1, x2, total_log_det

Because determinants of composed (chained) transformations multiply, and log turns products into sums, the total log-determinant of an L-layer flow is just the sum of each layer's own log det — which is why the accumulator above is a running sum rather than anything more complicated. This additivity is what makes deep flows, with dozens of coupling layers, exactly as tractable per-layer-cost as a single one.

Training. Unlike a GAN, a flow is trained with an entirely conventional objective: maximum likelihood. Given a batch of real data x⁽¹⁾, ..., x⁽ᵐ⁾, adjust the weights inside every s and t network to maximize Σᵢ log p_X(x⁽ⁱ⁾) = Σᵢ [log p_Z(f(x⁽ⁱ⁾)) + log det J_f(x⁽ⁱ⁾)] — precisely the quantity computed in the code above, summed over a batch, then differentiated with backpropagation like any other neural network loss. There is no adversarial game, no discriminator, and no approximate lower bound: the objective being optimized is the exact log-likelihood.

The Diagram

Correcting a Common Misconception

Because CBSE's Relations and Functions chapter introduces invertibility mainly through simple, mostly linear or monotone-textbook examples (f(x) = 2x+3, f(x) = eˣ, f(x) = log x), it is natural — and wrong — to assume that "invertible" secretly means "simple" or "close to linear." Students carry this assumption into coupling layers and get confused: how can y2 = x2·exp(s(x1)) + t(x1) be invertible if s and t are deep, wildly nonlinear neural networks with no invertibility of their own?

The resolution is that invertibility of the whole coupling layer never depended on s or t being invertible at all. Look again at the inverse formula: x2 = (y2 - t(x1))·exp(-s(x1)). Here, s(x1) and t(x1) are being evaluated, not inverted — once x1 = y1 is known, they are just two numbers you plug in. What actually needs to be invertible is the map from x2 to y2 for a fixed x1, which is x2 ↦ x2·(a positive constant) + (another constant) — an affine map with nonzero slope, guaranteed invertible by nothing more sophisticated than solving a linear equation for x2. The wild nonlinearity of s and t as functions of x1 is completely free — it never has to be undone, only recomputed. This is the real content of the coupling-layer trick: push all the expressive, hard-to-invert nonlinearity into a role (computing coefficients) where invertibility is simply never required, and keep the one part that must be inverted (the map in x2) deliberately as simple as an affine transform.

Where This Fits in Your Exams

The individual mathematical skills used in this chapter are squarely inside your existing syllabus, even though "normalizing flow" as a named model is not. Computing the determinant of the 2×2 Jacobian above is the same operation as any CBSE Class 12 Determinants-chapter question — you are filling a matrix with numbers (here, derivatives instead of plain constants) and applying the same rules, including the property that a triangular matrix's determinant is the product of its diagonal entries. The definitions of one-one, onto, and invertible functions come directly from the Class 12 Relations and Functions chapter, and recognizing why an affine map with positive slope is always invertible is exactly that chapter's content applied to a new setting. The differentiation used to build the Jacobian entries — product rule, chain rule, sum rule, applied together in one expression — is standard JEE Main/BITSAT-level differentiation practice, just assembled into a matrix instead of left as a single answer.

The 1-D change-of-variables derivation via the CDF is a genuine step up in rigor from board-exam probability, which stays mostly discrete and combinatorial (conditional probability, Bayes' theorem, binomial distributions). Continuous transformation-of-random-variables problems like the one worked out here are more characteristic of entrance exams that test calculus-based probability directly, such as the ISI B.Stat/B.Math and CMI entrance papers — if you are preparing for those, this derivation is worth redoing from scratch without looking. Beyond school and entrance exams, this exact machinery — change of variables, Jacobian determinants, log-likelihood as a training objective — reappears as-is in undergraduate coursework on deep generative modeling, so working through the derivation now rather than memorizing the final formula pays off later regardless of which specific course or programme you eventually take it in.

One historical note: KVPY (Kishore Vaigyanik Protsahan Yojana), which some older material still lists as a target exam for this kind of content, was discontinued starting the 2022 cycle and replaced by the INSPIRE scheme, which does not use a comparable written-test format — so it is not listed here as a live pathway.

Active Recall

1. Let X be Uniform on [0, 2] and Y = f(X) = X³. Derive p_Y(y) from the CDF method used in this chapter (state the domain of y too), then check your answer integrates to 1 over that domain.
2. For a coupling layer with s(x1) = -0.3·x1 and t(x1) = 0, and inputs x1 = 2, x2 = 5: compute y2, det J, and log det J by hand.
3. Explain in two or three sentences why stacking two coupling layers that both split the vector the same way, with no flip in between, fails to build a useful flow — which dimension never gets transformed, and why?
4. A classmate says the coupling layer "can't be that expressive because it's basically just multiplying by a number and adding a number." Using the misconception correction above, explain what part of the layer is actually doing the expressive work, and why the scale-and-shift part is deliberately kept simple.

Answer key — (1) x = y^{1/3}, dx/dy = (1/3)y^{-2/3}, so p_Y(y) = (1/2)·(1/3)y^{-2/3} = 1/(6·y^{2/3}) for 0 < y < 8; integrating, ∫₀⁸ (1/6)y^{-2/3} dy = (1/6)·[3y^{1/3}]₀⁸ = (1/6)(3·2) = 1. (2) exp(s(x1)) = exp(-0.6) ≈ 0.5488; y2 = 5·0.5488 + 0 ≈ 2.744; det J = exp(-0.6) ≈ 0.5488; log det J = -0.6. (3) The first half, x1, is copied straight to y1 in every layer since neither layer ever swaps which half is held fixed, so no transformation ever touches it — it stays exactly equal to the original input forever. (4) The expressive, hard-to-invert nonlinearity lives entirely inside s(x1) and t(x1), which can be arbitrarily deep networks; the scale-and-shift applied to x2 is kept deliberately affine only so that that one step stays trivially invertible, not because the whole layer is limited to affine expressiveness.

Think About It

Think about this: How would you explain normalizing flows: invertible transformations for generative modeling to a friend who has never seen a computer? What real-world analogy would you use? Imagine you had to build a system using these concepts — what would be your first step? Try this: before moving on, write down three things you learned and one question you still have.

Key Takeaways — Summary and Recap

Let us recap what we covered: the core ideas behind normalizing flows: invertible transformations for generative modeling, how they connect to real-world applications, and why they matter for your journey in computer science. Remember these key points as you move forward. For competitive exam preparation (CBSE, JEE, BITSAT), focus on understanding the WHY behind each concept, not just the WHAT.