Integrals: From Riemann Sums to Neural ODEs

A residual layer adds a small change to its current state. Let those changes become a continuous velocity field and the network becomes a differential equation. Recovering its final state requires accumulating change over time.

An area is a limit of sums

Partition $[a,b]$ into $n$ equal intervals of width $\Delta x=(b-a)/n$. Pick $x_i^*$ within each strip. For a continuous function,

$\int_a^b f(x)\,dx=\lim_{n\to\infty}\sum_{i=1}^n f(x_i^*)\Delta x.$

The integral is signed: contributions below the axis are negative. To measure unsigned area, integrate $|f|$. For $x^2$ on $[0,1]$, the right sum is $n^{-3}\sum_{i=1}^n i^2=(n+1)(2n+1)/(6n^2)$, tending to $1/3$.

The fundamental theorem

Define $F(x)=\int_a^x f(t)\,dt$ for continuous f. Then

$\frac{F(x+h)-F(x)}h=\frac1h\int_x^{x+h}f(t)\,dt.$

The average value over the shrinking interval tends to $f(x)$ by continuity. Thus $F'=f$. Conversely, if $A'=f$, the definite integral is $A(b)-A(a)$.

This converts accumulation into differentiation in reverse. It also differentiates integrals with moving boundaries: $d[\int_0^{x^2}f(t)dt]/dx=2xf(x^2)$ by the chain rule.

Substitution and integration by parts

For a differentiable change $x=g(u)$, substitution gives $\int f(g(u))g'(u)du=\int f(x)dx$. Definite integrals need transformed bounds. The formula comes directly from the chain rule.

Integrating the product rule gives $\int u\,dv=uv-\int v\,du$. For example, $\int_0^1 xe^x dx=[xe^x-e^x]_0^1=1$.

In multiple dimensions, volume changes by the absolute Jacobian determinant. Polar coordinates have area element $r\,dr\,d\theta$. This factor is essential in the Gaussian integral. Let $I=\int_{-\infty}^{\infty}e^{-x^2}dx$. Then, using positivity to combine integrals,

$I^2=\int_{\mathbb R^2}e^{-(x^2+y^2)}dxdy=\int_0^{2\pi}\int_0^\infty e^{-r^2}r\,dr\,d\theta=\pi.$

Since $I>0$, $I=\sqrt\pi$. Rescaling gives the standard normal normalising constant $\sqrt{2\pi}$.

Numerical accumulation

For sufficiently smooth functions and a fixed interval, left and right rules generally have error $O(n^{-1})$; midpoint and trapezoid have $O(n^{-2})$; composite Simpson has $O(n^{-4})$ and requires an even number of strips. Endpoint singularities can spoil these rates.

A Monte Carlo estimator uses uniform $U_i$ on $[a,b]$: $(b-a)n^{-1}\sum_i f(U_i)$. It is unbiased when the expectation exists, and has standard error proportional to $n^{-1/2}$ when the variance is finite. Its strength is handling high dimensions, rather than fast convergence for smooth one-dimensional curves.

From residual steps to an ODE

Euler's method for $z'=f(z,t)$ is $z_{k+1}=z_k+h f(z_k,t_k)$. A residual update has this form. For $z'=z$ with $z(0)=1$, n equal steps to time 1 give $(1+1/n)^n$, approaching $e$. Finer discretisation improves this example, but some ODEs require much smaller steps for stability.

Read beyond

Read the equation in context

Neural ODEs define the hidden-state evolution through a learned derivative and use an ODE solver to obtain later states. The integral contains the unknown trajectory z(t); it is not generally an ordinary integral of a known fixed curve. Solver tolerance and computational cost become part of evaluating the model.

Your turn

Compare left, midpoint, trapezoid and Simpson estimates. Reach an absolute error below 0.001 for three integrals with at most sixteen strips each.

Coding problems