Interactive Calculus

Definition. A function \(f:X\to Y\) is a rule assigning to each element of the domain \(X\) exactly one element of the codomain \(Y\). The condition "exactly one" is the essential part — if one input corresponds to two values, it is not a function.

A model is a function, in two senses. A neural network \(f_\theta(x)\) is a function of the input \(x\) and at the same time a function of the parameters \(\theta\). At inference we fix \(\theta\) and treat \(x\) as the variable; during training we fix \(x\) and treat \(\theta\) as the variable. This structure — the same expression doing entirely different work depending on which argument is regarded as variable — is exactly the same kind of perspective shift as "probability vs likelihood" in the probability article.

Composition is everything. The "deep" in deep learning means the depth of composition. Simple functions are stacked layer upon layer as \(f = f_L\circ\cdots\circ f_2\circ f_1\) to build a complicated one. Each layer is a combination of a linear transformation and a nonlinear activation, and if even one nonlinearity is missing the whole thing collapses into a single linear transformation however deep it is stacked.

Why smoothness is needed. Gradient descent works only if the function is differentiable. So the functions used in deep learning are almost all continuous and (almost everywhere) differentiable. Functions with a kink at one point, such as ReLU, are used too, but they are smooth away from that single point, so it causes no trouble in practice.

\[\begin{aligned} f &: X \to Y \\ (g\circ f)(x) &= g(f(x)) \\ f_\theta &: \mathbb{R}^{d} \to \mathbb{R}^{k} \quad \text{(model)} \end{aligned}\]

Pick a function and move its parameters. Even for the same expression, the graph responds differently depending on which variable you move. Turn on composition mode and f(x) and g(f(x)) are drawn together, showing how the shape is deformed when functions are stacked.

Definition. \(\lim_{x\to a}f(x)=L\) means that by bringing \(x\) sufficiently close to \(a\) (but with \(x\ne a\)) we can make \(f(x)\) as close to \(L\) as we like. What value the function takes at \(x=a\), or even whether it is defined there, is irrelevant to the limit.

Why the detour is necessary. A slope is \(\frac{\Delta y}{\Delta x}\), but to obtain the slope "at a point" we would need \(\Delta x=0\), which gives \(0/0\). A limit avoids this dead end — it sends \(\Delta x\) to 0 without ever setting it to 0. Shrinking \(h\) in the demo shows the slope of the secant approaching a value, and that value is the slope of the tangent.

Continuity. A function is continuous at \(a\) when \(\lim_{x\to a}f(x)=f(a)\); that is, when the limit and the function value agree. Without continuity there can be no derivative, but continuity does not guarantee differentiability — \(|x|\) is continuous at the origin, yet the slopes approaching from the left and the right are \(-1\) and \(+1\), so it is not differentiable.

In machine learning. Limits appear when we talk about the convergence of training. The statement "sending the learning rate to 0 makes SGD converge to a differential equation called gradient flow" is one example, and it is also the starting point of the next article. It matters numerically too — when approximating a gradient by finite differences rather than autodiff, too large an \(h\) gives a large approximation error while too small an \(h\) lets floating-point rounding error dominate.

\[\begin{aligned} \lim_{x\to a} f(x) &= L \\ f'(a) &= \lim_{h\to 0}\frac{f(a+h)-f(a)}{h} \\ \text{continuous} &\iff \lim_{x\to a}f(x) = f(a) \end{aligned}\]

Set the h slider large, then shrink it. The secant (orange) through two points gradually merges with the tangent (green). You can watch the secant slope in the readout converge to the true derivative, and also see the numerical error grow again when h becomes extremely small.

Definition. \(f'(a)=\lim_{h\to 0}\frac{f(a+h)-f(a)}{h}\). If this limit exists, \(f\) is said to be differentiable at \(a\). The notations \(\frac{df}{dx}\), \(\dot{f}\) and \(Df\) are several, but the meaning is one.

Geometrically: the tangent. \(f'(a)\) is the slope of the line touching the graph at \((a,f(a))\). That tangent \(y=f(a)+f'(a)(x-a)\) is the first-degree function best approximating \(f\) near \(a\), and it is the first-order term of the Taylor expansion two sections from here.

What the sign and magnitude say. \(f'>0\) means increasing, \(f'<0\) decreasing, and \(f'=0\) stationary (a critical point). The larger the absolute value, the steeper. In a minimization problem we move opposite to \(f'\) — which is \(\theta \leftarrow \theta - \eta f'(\theta)\), gradient descent. In the derivative graph below you can see the peaks and valleys of the original function line up exactly with the points where \(f'=0\).

A critical point is not a minimum. \(f'=0\) holds at minima, maxima and saddle points alike. Telling them apart requires second-order information, which is the subject of the Hessian section. In high-dimensional loss landscapes saddle points are known to be far more common than local minima, and this fact directly influences the design of optimization algorithms.

Automatic differentiation. What PyTorch does is not finite differences. It knows the derivative rule for each primitive operation and applies the chain rule mechanically along the computation graph. So the error is at the level of machine precision rather than numerical approximation, and every partial derivative is obtained in a single backward traversal even for hundreds of millions of parameters.

\[\begin{aligned} f'(a) &= \lim_{h\to 0}\frac{f(a+h)-f(a)}{h} \\ y &= f(a) + f'(a)(x-a) \quad \text{(tangent)} \\ \theta &\leftarrow \theta - \eta\, f'(\theta) \quad \text{(gradient descent)} \end{aligned}\]

Drag the point a. A tangent is drawn on the upper graph and the corresponding value is marked on the derivative graph below. Check the correspondence: where the original function flattens, the derivative crosses 0. The "One step of gradient descent" button shows the actual descent toward a minimum.

Definition. If \(y=g(f(x))\) then \(\frac{dy}{dx}=g'(f(x))\cdot f'(x)\). In Leibniz notation, with \(u=f(x)\), this reads \(\frac{dy}{dx}=\frac{dy}{du}\cdot\frac{du}{dx}\), a form that looks as though fractions were cancelling.

Intuition: amplification of rates. Move \(x\) by 1 and \(u\) moves by \(f'(x)\); move \(u\) by that much and \(y\) moves by a further factor of \(g'(u)\). It is like meshing gear ratios in series, so the overall ratio is the product of the individual ratios.

Backpropagation is exactly this. In \(L = \ell(f_L(\cdots f_1(x)))\), the gradient with respect to \(\theta_k\) is the product of all the derivatives from layer \(k\) to the output. Storing each layer's local derivative on the forward pass and multiplying them out in a single backward sweep yields the gradients for every layer in one traversal. What would be \(O(L^2)\) work if computed layer by layer becomes \(O(L)\).

Vanishing and exploding gradients. Being a product has a price. If each stage's derivative averages 0.5, then after 20 layers it has vanished to \(0.5^{20}\approx 10^{-6}\); at 1.5 it explodes to \(1.5^{20}\approx 3300\). The fact that the maximum slope of a sigmoid is 0.25 blocked early deep-network training, and ReLU (slope 1 on the positive side), residual connections (preserving a gradient of 1 along an identity path) and normalization are all devices for keeping this product near 1.

Extension to several variables. When inputs and outputs are vectors, each stage's derivative is a Jacobian matrix rather than a number, and the chain rule becomes matrix multiplication: \(J_{g\circ f}=J_g\,J_f\). What backpropagation actually does is multiply these matrices, but without ever forming them — computing only vector–Jacobian products (VJPs) to save memory.

\[\begin{aligned} \frac{dy}{dx} &= \frac{dy}{du}\cdot\frac{du}{dx} = g'(f(x))\,f'(x) \\ \frac{\partial L}{\partial \theta_k} &= \frac{\partial L}{\partial h_L}\prod_{i=k+1}^{L}\frac{\partial h_i}{\partial h_{i-1}}\cdot\frac{\partial h_k}{\partial \theta_k} \\ J_{g\circ f} &= J_g\, J_f \end{aligned}\]

Adjust each layer's gradient factor with the slider and increase the number of layers. Below 1, the accumulated gradient vanishes rapidly toward 0; above 1, it explodes. Only near 1 does a deep network stay in a trainable range — which is why ResNets and normalization techniques exist.

Definition. Near \(a\), \(f(x)\approx \sum_{n=0}^{N}\frac{f^{(n)}(a)}{n!}(x-a)^n\). For \(N=1\) this is the tangent line, for \(N=2\) a tangent parabola, and raising \(N\) brings it closer to the original function (within the radius of convergence).

Why divide by \(n!\)? Differentiating \((x-a)^n\) \(n\) times produces an \(n!\). Dividing it out in advance is what makes the \(n\)-th coefficient match \(f^{(n)}(a)\) exactly at \(x=a\). In other words, a Taylor polynomial is a polynomial engineered to agree with the original function at \(a\) through the \(N\)-th derivative.

Gradient descent = trusting the first-order approximation. In \(f(\theta+\Delta)\approx f(\theta)+f'(\theta)\Delta\), reducing the loss means giving \(\Delta\) the opposite sign to \(f'\). But the first-order approximation is only valid near \(\theta\), so we cannot move far in one go — which is exactly why a learning rate exists. The learning rate sets "how far do we trust this linear approximation."

Newton's method = trusting the second-order approximation. Minimizing \(f(\theta+\Delta)\approx f+f'\Delta+\tfrac{1}{2}f''\Delta^2\) over \(\Delta\) gives \(\Delta=-f'/f''\). No learning rate need be chosen; the curvature sets the step — large where the ground is flat (small \(f''\)), small where it is steep. In several variables \(f''\) becomes the Hessian matrix, and the cost of inverting it is why deep learning uses approximations (L-BFGS, Adam's diagonal approximation, and so on).

It does not work everywhere. A Taylor expansion is a local approximation. It diverges rapidly away from \(a\), and for some functions the radius of convergence is itself finite. Moving the expansion point in the demo shows the interval of good agreement travelling with it.

\[\begin{aligned} f(x) &\approx \sum_{n=0}^{N}\frac{f^{(n)}(a)}{n!}(x-a)^n \\ \Delta_{\text{GD}} &= -\eta\, f'(\theta) &&\text{(first order)} \\ \Delta_{\text{Newton}} &= -\frac{f'(\theta)}{f''(\theta)} &&\text{(second order)} \end{aligned}\]

Raise the degree N from 0. The approximation wraps around the original function in the order constant → tangent → parabola → … Drag the expansion point a and the interval of good agreement follows. The error in the readout shows how quickly it breaks down as you move away from a.

Definition (definite integral). \(\int_a^b f(x)\,dx\) is the limit of a sum of areas of rectangles obtained by chopping the interval finely — the limit of a Riemann sum. Geometrically it is the signed area between the \(x\)-axis and the curve: where the curve lies below the axis it counts as negative.

The fundamental theorem of calculus. Differentiation and integration are inverse operations. Defining \(F(x)=\int_a^x f(t)\,dt\) gives \(F'(x)=f(x)\), and therefore \(\int_a^b f = F(b)-F(a)\). This discovery — that the problem of finding an area and the problem of finding a slope are the same problem — is the heart of calculus. The cumulative curve below rising and falling according to the sign of the curve above shows the theorem directly.

Why area and slope are inverse operations. The theorem can be accepted from a single picture before following any proof. Read \(F(x)=\int_a^x f\) as "the area accumulated so far." If we extend \(x\) by a tiny \(dx\), how much does the area grow? The newly added sliver is a thin rectangle of width \(dx\) and height \(f(x)\), so \(dF\approx f(x)\,dx\), and dividing both sides by \(dx\) gives \(F'(x)=f(x)\). That is, the rate at which the accumulated quantity grows is precisely the height being laid down at this instant. It is like a reservoir whose rate of filling is exactly the flow currently coming out of the tap — which is why "accumulating" and "measuring a rate" undo each other.

So \(F(b)-F(a)\) is natural too. Chop the interval finely and add up the area \(dF\) gained in each piece, and the intermediate terms cancel telescopically, leaving only the start and the end. In the reservoir analogy: "the total that flowed in over the interval = final level − initial level." The familiar procedure of finding an antiderivative and subtracting at the two ends is a shortcut that checks only the beginning and end of the accumulated quantity instead of counting the area directly.

In machine learning: mostly incomputable. The expectation \(\mathbb{E}[f(X)]=\int f(x)p(x)\,dx\), marginalization \(p(x)=\int p(x,z)\,dz\), and the evidence term of Bayes' rule are all integrals. Yet in high dimensions these integrals almost never have a closed form. So there are two standard responses in practice — Monte Carlo (approximate by a sample mean) and variational inference (replace by a computable lower bound, that is the ELBO).

Numerical integration. In low dimensions we can compute directly. The trapezoid rule approximates each interval by a straight line and Simpson's rule by a parabola, the latter converging much faster (error \(O(h^4)\) against \(O(h^2)\)). But as the dimension rises the number of grid points grows exponentially and the approach becomes unusable — which is why Monte Carlo wins in high dimensions.

\[\begin{aligned} \int_a^b f(x)\,dx &= \lim_{n\to\infty}\sum_{i=1}^{n} f(x_i)\,\Delta x \\ \frac{d}{dx}\int_a^x f(t)\,dt &= f(x) \quad \text{(fundamental theorem)} \\ \mathbb{E}_p[f(X)] &= \int f(x)p(x)\,dx \approx \frac{1}{N}\sum_i f(x^{(i)}) \end{aligned}\]

Move the endpoints of the integration interval with the sliders and adjust the number of subdivisions n. You can watch the sum of rectangles converge to the true value, and see intervals where the curve dips below the axis subtract from the area. That the slope of the cumulative curve below matches the value of the curve above is the fundamental theorem.

Definition. The double integral \(\iint_A f(x,y)\,dA\) of a two-variable function is the limit obtained by chopping the region \(A\) into small rectangles, placing on each a column of volume \(f(x,y)\,\Delta x\,\Delta y\), and adding them all. What was "collect rectangles into an area" in one variable becomes "collect boxes into a volume" here.

Computation is ultimately repeated single-variable integration. By Fubini's theorem a double integral is computed as an iterated integral, working from the inside out. Integrating \(y\) first gives, for each \(x\), the area of a vertical cross-section \(A(x)=\int f(x,y)\,dy\), and integrating that cross-sectional area over \(x\) gives the volume — like slicing a loaf thinly, finding the area of each slice, and multiplying by the thickness before adding. Swapping the order of integration in the demo confirms that the final value is the same.

When the region is not a rectangle. The limits of the inner integral become functions of the outer variable. For instance, if \(x\in[0,1]\) and for each \(x\) the variable \(y\) runs from 0 to \(2-2x\), we write \(\int_0^1\!\!\int_0^{2-2x} f\,dy\,dx\). Swapping the order then requires re-deriving the limits, and it is common for an integral that resists one order to fall out easily in the other.

In machine learning. The condition that a joint distribution integrates to 1, marginalization \(p(x)=\int p(x,y)\,dy\), and expectations are all multiple integrals. The probability article's description of marginalization as "an operation that flattens along one axis" reveals its true identity here as one concrete inner integration. And as before, as the dimension rises this grid computation becomes exponentially expensive and one must switch to Monte Carlo.

\[\begin{aligned} \iint_A f(x,y)\,dA &= \int_a^b\!\!\left[\int_{y_1(x)}^{y_2(x)} f(x,y)\,dy\right]dx \\ &= \int_c^d\!\!\left[\int_{x_1(y)}^{x_2(y)} f(x,y)\,dx\right]dy \\ p(x) &= \int p(x,y)\,dy \quad \text{(marginalization = the inner integral)} \end{aligned}\]

Move the slice position and the cross-section at that location is drawn on the lower graph, with its area displayed. Collecting and adding those areas is the total volume. Check that swapping the order of integration gives the same value.

Definition. \(\frac{\partial f}{\partial x}(a,b)=\lim_{h\to 0}\frac{f(a+h,b)-f(a,b)}{h}\). Treat the other variables as constants and differentiate with respect to one only. The symbol is \(\partial\) rather than \(d\) to signal "there are other variables, but they are held fixed for now."

Geometrically: the slope of a cross-section. The graph of a two-variable function is a surface. Cutting at \(y=b\) slices the surface into a single curve, and the slope of that curve is \(\partial f/\partial x\). Cutting at \(x=a\) gives a different curve whose slope is \(\partial f/\partial y\). That the slope at the same point differs depending on the direction of the cut is the starting point of multivariable differentiation.

Generalizing to directional derivatives. There is no reason for the coordinate-axis directions to be special. The rate of change along an arbitrary unit vector \(u\) is \(D_u f = \nabla f\cdot u\), and partial derivatives are the special cases where \(u\) is \(e_1\) or \(e_2\). This formula determines the meaning of the gradient in the next section.

In machine learning. The gradient each parameter receives from the loss \(L(\theta_1,\ldots,\theta_d)\) is \(\partial L/\partial \theta_i\). With hundreds of millions of parameters we need hundreds of millions of partial derivatives, and obtaining them one at a time by finite differences would require hundreds of millions of forward passes. That backpropagation solves this in a single backward traversal is the decisive fact that made deep learning computable.

\[\begin{aligned} \frac{\partial f}{\partial x}(a,b) &= \lim_{h\to 0}\frac{f(a+h,b)-f(a,b)}{h} \\ D_u f &= \nabla f \cdot u, \quad \|u\|=1 \\ \frac{\partial L}{\partial \theta_i} &\;\; i=1,\ldots,d \end{aligned}\]

Drag the point on the contour plot. The cross-sectional slopes along the two axis directions are marked with arrows on the upper figure, and the two graphs below draw the cross-section curve in each direction with its tangent. Note that at the same point the slopes in x and y can be entirely different.

Definition. \(\nabla f = \left(\frac{\partial f}{\partial x_1},\ldots,\frac{\partial f}{\partial x_d}\right)\). A vector whose components are all the partial derivatives of a scalar function; it lives in the same dimension as the input.

Why is it the steepest-ascent direction? The problem is to maximize the directional derivative \(D_u f=\nabla f\cdot u\) over unit vectors \(u\). Since the inner product is \(\nabla f\cdot u=\|\nabla f\|\cos\theta\), it is largest at \(\theta=0\), that is when \(u\) points the same way as \(\nabla f\). So to decrease fastest one must go in the direction \(-\nabla f\), and that is gradient descent. A single line of inner product from linear algebra fixes the basic direction of optimization.

It is orthogonal to the contours. Moving along a contour leaves the function value unchanged, so the directional derivative in that direction is 0, that is \(\nabla f\cdot u=0\). The gradient is therefore always perpendicular to the contours. Moving the point in the demo shows the arrow always pointing across the contour lines.

The gradient is a "direction," not a "destination." \(-\nabla f\) is only the steepest direction at the current position, not the direction toward the minimum. In a long narrow valley the gradient bounces zigzag against the valley walls and fails to travel along the floor. Momentum, Adam, and second-order methods are all attempts to correct this defect — making the landscape elongated in the demo lets you watch the zigzag directly.

A note on notation. \(\nabla f\) is defined only for scalar functions. The derivative of a vector-valued function is not a gradient but the Jacobian matrix of the next section. A neural network's loss is a scalar, which is why \(\nabla_\theta L\) is well defined and why training becomes a problem of following a single vector.

\[\begin{aligned} \nabla f &= \left(\tfrac{\partial f}{\partial x_1},\ldots,\tfrac{\partial f}{\partial x_d}\right) \\ D_u f &= \nabla f\cdot u = \|\nabla f\|\cos\theta \\ \theta &\leftarrow \theta - \eta\,\nabla_\theta L \end{aligned}\]

Drag the starting point and press "Run descent." The trajectory crosses the contours on its way down to the minimum. Use the anisotropy slider to elongate the valley and the trajectory starts bouncing zigzag — the fundamental weakness of pure gradient descent.

Definition. The Jacobian of \(F:\mathbb{R}^n\to\mathbb{R}^m\) is the \(m\times n\) matrix whose \((i,j)\) entry is \(\partial F_i/\partial x_j\). For \(m=1\) it is the transpose of the gradient, and for \(n=m=1\) it is just the derivative — the Jacobian is the most general form of the notion of a derivative.

Locally, everything is a linear transformation. Near a point \(a\), \(F(a+\Delta)\approx F(a)+J_F(a)\,\Delta\). That is, even a nonlinear function acts as a single matrix when viewed closely enough. This is the opening claim of this article, "look closely enough and everything is a straight line," carried into several dimensions, and the deformation of the grid seen in the linear transformation article is reproduced here exactly.

Chain rule = matrix multiplication. \(J_{g\circ f}(x)=J_g(f(x))\,J_f(x)\). Backpropagation is a chain of these matrix products, but real implementations never build the matrices. Since the loss is a scalar, taking successive vector–Jacobian products \(v^\top J\) from the left suffices, and this approach (reverse mode) is overwhelmingly favourable for deep learning, where parameters are many and the output is one.

The determinant: a volume factor. \(|\det J|\) says by what factor an infinitesimal volume changes near that point. It means exactly what the determinant meant as an area factor in the linear algebra article, differing only in that it varies from point to point. At points where \(\det J=0\) the dimension collapses locally and the function has no inverse there.

In machine learning: normalizing flows. Transforming a random variable by \(z\mapsto x=F(z)\) changes the density to \(p_x(x)=p_z(z)\,|\det J_F(z)|^{-1}\) — the density thins in proportion to how much the volume expanded. Normalizing flows are the family of models that design \(F\) so that \(\det J\) is cheap to compute (triangular Jacobians, coupling layers, and so on) and thereby obtain an exact log-likelihood. The probability-density flow of the next article is ultimately the continuous-time version of this relationship.

\[\begin{aligned} (J_F)_{ij} &= \frac{\partial F_i}{\partial x_j}, \qquad J_F \in \mathbb{R}^{m\times n} \\ F(a+\Delta) &\approx F(a) + J_F(a)\,\Delta \\ J_{g\circ f} &= J_g\,J_f \\ p_x(x) &= p_z(z)\,\bigl|\det J_F(z)\bigr|^{-1} \end{aligned}\]

Drag the point on the left-hand grid and the figure shows where it is carried on the right, and into what parallelogram the small square around it is deformed. The area ratio of that parallelogram is |det J|. Change the transformation and hunt for places that stretch and places that fold (det J = 0).

Definition. \((H_f)_{ij}=\frac{\partial^2 f}{\partial x_i\partial x_j}\). The \(d\times d\) matrix collecting the second partial derivatives; if \(f\) is smooth enough the order of mixed partials does not matter (Schwarz's theorem), so it is symmetric. Symmetry is what guarantees that the eigenvalues are all real and the eigenvectors can be chosen orthogonal.

The coefficients of the second-order approximation. \(f(a+\Delta)\approx f(a)+\nabla f\cdot\Delta+\tfrac{1}{2}\Delta^\top H\Delta\). This is the multivariable version of the Taylor expansion, and the final quadratic form determines the bowl shape near that point. It is also why contours look like ellipses, and the axes of those ellipses are the eigenvectors of \(H\) — the same geometry as the covariance ellipse of the probability article, reappearing here.

Classifying critical points. At a point where \(\nabla f=0\), the signs of the Hessian's eigenvalues give the answer. All positive means a local minimum (curving upward in every direction), all negative a local maximum, and mixed signs a saddle point. In high dimensions the probability that all eigenvalues share a sign drops sharply, leading to the conclusion that most critical points of a deep learning loss landscape are saddle points.

Why the signs of the eigenvalues give the answer. At a critical point \(\nabla f=0\), so the first-order term of the second-order approximation vanishes and only \(f(a+\Delta)-f(a)\approx\tfrac12\Delta^\top H\Delta\) remains. The question therefore reduces to "is this quadratic form positive in every direction?" But \(H\) is symmetric, so we may change coordinates to its orthogonal eigenvector axes, and in those axes all cross terms disappear, leaving the simple sum \(\tfrac12(\lambda_1 c_1^2+\lambda_2 c_2^2+\cdots)\). Since \(c_i^2\) is always non-negative, the sign of each term is decided solely by \(\lambda_i\). The Hessian, in short, records with a single number \(\lambda_i\) per eigenvector direction whether the landscape curves up or down there, and the classification rule is just a way of reading it.

The name "saddle." Mixed signs means a valley along one axis and a peak along another. The shape resembles a horse's saddle, rising front to back and falling side to side, hence the name. Such a point has \(\nabla f=0\) and so looks "flat" to gradient descent, yet it is not a minimum; escaping requires finding a direction with \(\lambda<0\). This is why pure gradient descent lingers near saddle points, and also why the noise in SGD helps by pushing it along such directions.

The condition number = the difficulty of optimization. The ratio of largest to smallest eigenvalue, \(\kappa=\lambda_{\max}/\lambda_{\min}\), is the condition number. A large \(\kappa\) makes the contours a long narrow valley and gradient descent falls into the zigzag seen in the previous section, because the stable learning rate is limited by \(\lambda_{\max}\) while the speed of convergence is determined by \(\lambda_{\min}\). One reason normalization techniques speed up training is that they lower this condition number.

In machine learning. With \(d\) in the hundreds of millions, \(H\) has \(d^2\) entries and cannot even be formed. So practice goes around it — computing only Hessian–vector products (possible without the full matrix), approximating only the diagonal as Adam does, or building a low-rank approximation from recent gradients as L-BFGS does. Meanwhile the eigenvalue distribution of the Hessian is actively used as a measure of the flatness (sharpness) of a loss landscape in generalization research.

\[\begin{aligned} (H_f)_{ij} &= \frac{\partial^2 f}{\partial x_i \partial x_j} = (H_f)_{ji} \\ f(a+\Delta) &\approx f(a) + \nabla f^\top \Delta + \tfrac{1}{2}\Delta^\top H \Delta \\ \kappa &= \lambda_{\max}/\lambda_{\min} \end{aligned}\]

Move between bowl, saddle and valley with the presets and drag the point. The two eigenvector axes (red and violet) and their eigenvalues are displayed, and the type of critical point is determined by the combination of signs. Check how large the condition number becomes on the valley preset, and how elongated the contours are then.

Definition. A vector field is a function \(F:\mathbb{R}^n\to\mathbb{R}^n\), and the essential point is that the input and output dimensions are equal. Read it as assigning to each point \(x\) an arrow \(F(x)\) attached at that point. The velocity field of a fluid, a gravitational field, and the \(-\nabla L\) of the previous section are all vector fields.

A gradient field is a special vector field. A vector field expressible as the gradient of some scalar function is called conservative. Such a field has a "potential," so the arrows always point uphill, and going around a closed loop produces no net change, so it does not circulate. A rotating vector field, by contrast, cannot be written as the gradient of any scalar. Gradient descent on a loss does not spin in circles but eventually descends somewhere because \(-\nabla L\) is conservative — whereas the vector field created by the two players of a GAN is not conservative, so circulation arises and training oscillates.

Streamlines. The curve drawn by continuing to follow the arrow direction at each point is a streamline, and it is the solution of the differential equation \(\dot{x}=F(x)\). That is, drawing a vector field and writing down a differential equation are the same act. Clicking anywhere in the demo draws the streamline starting from that point.

In machine learning. Training dynamics is itself a flow on a vector field. In the limit of the learning rate going to 0, gradient descent becomes the differential equation \(\dot\theta=-\nabla L(\theta)\) (gradient flow). Neural ODEs make the network itself the vector field and solve \(\dot{h}=f_\theta(h,t)\), and the reverse process of a diffusion model and the Probability Flow ODE are likewise flows carrying noise to data along a learned vector field.

\[\begin{aligned} F &: \mathbb{R}^n \to \mathbb{R}^n \\ \dot{x} &= F(x) \quad \text{(the streamline is the solution)} \\ \dot\theta &= -\nabla L(\theta) \quad \text{(gradient flow)} \end{aligned}\]

Compare several vector fields with the presets. Clicking anywhere on the figure draws the streamline starting there. In a gradient field the streamlines cross the contours and must converge somewhere, while in a rotational field they circle forever — that difference is the subject of the next section on divergence and curl.

Definition. \(\nabla\cdot F=\frac{\partial F_1}{\partial x_1}+\cdots+\frac{\partial F_n}{\partial x_n}\). It equals the trace of the Jacobian, and it takes a vector field in and returns a scalar field.

Geometrically: sources and sinks. Place a very small box around a point and take the limit of the difference between what enters and what leaves, divided by the box's volume — that is the divergence. If \(\nabla\cdot F>0\), flow wells up at that point (a source); if \(<0\), it is drawn in (a sink); if \(=0\), as much leaves as enters. In the demo the red regions are sources and the blue regions sinks.

It is a rate of volume change. Think of a small droplet drifting along the vector field, and the relative rate of change of its volume is exactly the divergence: \(\frac{d}{dt}\log V = \nabla\cdot F\). Where the previous section's \(|\det J|\) was a discrete volume factor, divergence is its continuous-time version — indeed \(\frac{d}{dt}\log\det J = \nabla\cdot F\) holds.

The continuity equation. If matter is neither created nor destroyed, the density \(\rho\) satisfies \(\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho F)=0\). It is the conservation law "the rate at which density decreases = the amount flowing out," and in the next article this equation is used verbatim to describe how probability density flows. In the end the Fokker–Planck equation and the Probability Flow ODE are variations on this one line.

In machine learning. In a Continuous Normalizing Flow the change in log-density is given by \(\frac{d\log p}{dt}=-\nabla\cdot f_\theta\) — no determinant need be computed, only a trace, and even that trace can be cheaply approximated by a Hutchinson estimator, which is what made this family of models practical.

\[\begin{aligned} \nabla\cdot F &= \sum_i \frac{\partial F_i}{\partial x_i} = \operatorname{tr}(J_F) \\ \frac{d}{dt}\log V &= \nabla\cdot F \\ \frac{\partial \rho}{\partial t} + \nabla\cdot(\rho F) &= 0 \quad \text{(continuity equation)} \end{aligned}\]

The background colour is the divergence — red where flow wells up, blue where it is drawn in. Drag the point and the figure also shows how a small circle around it grows or shrinks as it follows the flow. In a rotational field the arrows spin fiercely and yet the divergence is 0 everywhere — rotation and divergence measure different things.

Definition. \(\Delta f=\nabla\cdot(\nabla f)=\sum_i \frac{\partial^2 f}{\partial x_i^2}\). Form the gradient (a vector field) and then take its divergence (a scalar); it equals the trace of the Hessian: \(\Delta f=\operatorname{tr}(H_f)\). It is also written \(\nabla^2 f\).

The key intuition: the difference from the surrounding average. The Laplacian is proportional to "how much smaller this point's value is than the average around it." Precisely, the mean value over a small sphere of radius \(r\) is \(\bar{f}\approx f(x)+\frac{r^2}{2n}\Delta f(x)\). So \(\Delta f>0\) means the point is a dip below its surroundings, and \(\Delta f<0\) means a bump above them. This one sentence explains every equation in which the Laplacian appears.

Which is where diffusion comes from. The heat equation \(\frac{\partial u}{\partial t}=\alpha\Delta u\) is a direct transcription of "places colder than their surroundings warm up and hotter places cool down." As time passes the bumps are shaved away and the dips filled in until everything is flat. The heat and diffusion equations of the next article, and the forward process of a diffusion model, are all this smoothing.

Harmonic functions. A function with \(\Delta f=0\) is called harmonic, and its value at every point equals exactly the average around it. Such a function can have no extremum in the interior (the maximum principle), and the boundary values completely determine the interior.

In machine learning. The Laplacian filter in image processing detects edges, because a point whose value departs greatly from its surrounding average is precisely a boundary. The graph Laplacian \(L=D-A\) carries this operator onto a graph and forms the theoretical basis of spectral clustering and GNNs — its eigenvectors correspond to "frequencies" on the graph. This term also appears in the Fokker–Planck equation corresponding to the forward SDE of a diffusion model, where it plays the role of blurring the data distribution into a Gaussian.

\[\begin{aligned} \Delta f &= \nabla\cdot(\nabla f) = \sum_i \frac{\partial^2 f}{\partial x_i^2} = \operatorname{tr}(H_f) \\ \bar{f}_{\text{sphere}} &\approx f(x) + \tfrac{r^2}{2n}\Delta f(x) \\ \frac{\partial u}{\partial t} &= \alpha\,\Delta u \quad \text{(heat equation)} \end{aligned}\]

The background colour is the Laplacian — blue for dips below the surroundings (to be filled in), red for bumps above them (to be shaved away). Press "Play diffusion" to watch the contours actually flatten. The order in which the bumps collapse first and the dips fill in matches the sign of the Laplacian exactly.