Mathematics for science

Functions

Science is largely the study of how one quantity depends on another. The pressure of a gas depends on its temperature, the energy of two atoms depends on the distance between them, the number of undecayed nuclei in a sample depends on how long you have waited. The mathematical object that captures “depends on” is the function.

A function is a rule that takes an input and produces exactly one output. If the rule is called $f$ and the input is $x$, the output is written $f(x)$, read “f of x”. The insistence on exactly one output is the whole point: if you know the input, the output is determined. The rule “square the number” is a function, because every number has one square. The rule “give me a number whose square is $x$” is not, because 4 has two square roots, $2$ and $-2$, and the rule does not say which to choose.

It helps to picture a function as a machine. A number goes in on one side, the machine does something to it, and a number comes out of the other side. The machine does not remember previous inputs and does not behave differently on different days. The figure below lets you feed numbers into such a machine. Each time a number comes out, the pair (input, output) is recorded as a point on a graph, and if you feed many inputs the points join up into a curve. That curve is the graph of the function, and it is nothing more than the record of every input–output pair.

Some vocabulary makes it possible to talk about functions precisely. The set of inputs a function accepts is its domain. The square root of a negative number is not a real number, so the domain of $\sqrt{x}$ is $x \ge 0$; the function $1/x$ accepts every number except zero. The set of outputs a function can produce is its range: $x^2$ never produces a negative number, so its range is $y \ge 0$. The input is often called the independent variable and the output the dependent variable, because the output depends on the input and not the other way round.

Functions can be described in four equivalent ways, and moving between them is a skill worth practising. A formula such as $f(x) = 1.8x + 32$ gives the rule symbolically. A table lists selected input–output pairs, which is how experimental data usually arrive. A graph shows every pair at once as a curve. And a description in words says what the rule does: “convert a temperature from degrees Celsius to degrees Fahrenheit”. Each form is good for something different. The formula is good for calculation, the table for recording measurements, the graph for seeing behaviour at a glance, and the words for remembering what the function means.

Functions in science carry units. In $T_F = 1.8\,T_C + 32$ the input is in °C and the output in °F, so the number 1.8 has units of °F per °C and the 32 is in °F. Keeping track of units is the single most effective way to catch mistakes in a formula: if two terms that are added together have different units, the formula is wrong.

Reading a graph

A graph is drawn on two perpendicular axes. The horizontal axis carries the input, conventionally $x$, and the vertical axis the output, $y = f(x)$. The point where the axes cross is the origin, $(0, 0)$. Every point on the curve has coordinates $(x, f(x))$.

Several features of a function can be read directly from its graph. Where the curve crosses the vertical axis, the input is zero, so that height is $f(0)$, the intercept. Where the curve crosses the horizontal axis the output is zero; these inputs are the roots or zeros of the function and solving the equation $f(x) = 0$ means finding them. Where the curve rises from left to right the function is increasing, and where it falls it is decreasing. The highest and lowest points of a hump or a dip are maxima and minima, and they matter enormously in science: a stable arrangement of atoms is a minimum of the energy, and the barrier a chemical reaction must cross is a maximum along the reaction path. A curve may also approach a straight line without ever reaching it, as $1/x$ approaches the horizontal axis when $x$ becomes large; such a line is an asymptote.

Graphs can be transformed in predictable ways, and recognising these transformations makes new functions look familiar. Replacing $f(x)$ by $f(x) + c$ moves the whole graph up by $c$. Replacing $f(x)$ by $f(x - h)$ moves it to the right by $h$, which surprises many people at first: the output that used to occur at $x = 0$ now occurs at $x = h$. Multiplying the output, $a\,f(x)$, stretches the graph vertically by the factor $a$, and a negative $a$ also flips it upside down. Multiplying the input, $f(kx)$, squeezes the graph horizontally by the factor $k$. A Gaussian peak centred at $\mu$ with width $\sigma$, for example, is just the standard bell curve shifted by $\mu$ and stretched by $\sigma$.

Linear functions

The simplest non-trivial function is the straight line,

The slope is the defining property of a straight line: between any two points $(x_1, y_1)$ and $(x_2, y_2)$ on it, the ratio $(y_2 - y_1)/(x_2 - x_1)$ is the same number $m$. “Rise over run” is the usual mnemonic. A positive slope means $y$ grows with $x$, a negative slope means it shrinks, and a zero slope gives a horizontal line. The slope has units, output units per input unit, and it often has a physical meaning of its own: in a graph of position against time the slope is a velocity, in a graph of charge against voltage it is a capacitance.

When $c = 0$ the line passes through the origin and $y$ is proportional to $x$, written $y \propto x$: doubling $x$ doubles $y$. Many laws of physics take this form near equilibrium. Hooke’s law says the restoring force of a spring is proportional to its extension, $F = -kx$. Ohm’s law says the current through a resistor is proportional to the voltage across it. Linear behaviour is so common because almost any smooth function looks like a straight line if you only look at a small enough piece of it, an idea that becomes the derivative in the next chapter.

Linear functions are also the workhorse of data analysis. Plotting data so that a law appears as a straight line, and then measuring the slope and intercept, is one of the oldest and most reliable ways of extracting physical constants from experiment. The next sections show how to do this even when the underlying law is not linear at all.

Powers and polynomials

A power function has the form $y = a x^n$. For whole-number powers this is familiar: $x^2$ is a parabola, symmetric about the vertical axis, and $x^3$ rises steeply on the right and falls steeply on the left. Negative powers describe quantities that fall off with distance. The electrostatic energy of two charges varies as $1/r$, the gravitational force as $1/r^2$, and the repulsion between two atoms when their electron clouds overlap rises roughly as $1/r^{12}$ in the Lennard-Jones model. Fractional powers also appear: $x^{1/2} = \sqrt x$ describes, for instance, how far a diffusing molecule typically travels in time $x$.

Power functions have a scaling property that makes them special. If $y = a x^n$ and you multiply $x$ by a factor $\lambda$, then $y$ is multiplied by $\lambda^n$, whatever the starting value of $x$. Doubling the distance between two charges halves their energy, everywhere. This is why power laws describe phenomena that have no built-in length or time scale.

A sum of whole-number powers, $p(x) = a_0 + a_1 x + a_2 x^2 + \dots + a_N x^N$, is a polynomial of degree $N$. Polynomials are easy to evaluate, differentiate and integrate, and a polynomial of degree $N$ has at most $N$ real roots and at most $N - 1$ turning points. They are the raw material of approximation: Taylor series build complicated functions from polynomials, and the harmonic approximation of a chemical bond replaces its true energy curve by a parabola near the minimum.

Exponential functions

Suppose a quantity changes at a rate proportional to how much of it there is. A population of bacteria with plenty of food grows at a rate proportional to the number of bacteria; a sample of radioactive nuclei decays at a rate proportional to the number of nuclei left; the light passing through an absorbing liquid loses intensity at a rate proportional to the intensity. All of these are described by the exponential function:

The defining feature of an exponential is that adding a fixed amount to the input multiplies the output by a fixed factor. For a straight line, each step of one unit in $x$ adds $m$ to $y$; for an exponential, each step multiplies $y$ by $e^k$. As a result, an exponential always has a fixed doubling time (for growth) or half-life (for decay), equal to $\ln 2/|k|$. Carbon-14 has a half-life of about 5730 years: after 5730 years half of the original nuclei remain, after 11 460 years a quarter, after 17 190 years an eighth, and so on, regardless of how much you started with. The golden dots in the figure above mark these successive halvings when you select the exponential family.

Why the strange base $e$? Any base could be used, since $2^x$, $10^x$ and $e^x$ are all exponentials and each can be converted into the others. The number $e$ is chosen because it makes the mathematics of change simplest: the slope of $e^x$ at any point equals its height. This is proved in the chapter on derivatives, and it is why nature’s exponentials, which arise from rates proportional to amounts, appear with base $e$.

Exponentials are the most important functions in statistical mechanics. The probability that a molecule has energy $E$ at temperature $T$ is proportional to the Boltzmann factor $e^{-E/k_BT}$, and the rate of a chemical reaction with activation energy $E_a$ follows the Arrhenius law $k = A\,e^{-E_a/k_BT}$. Both are decaying exponentials of energy: each additional $k_BT$ of energy makes a state $e \approx 2.7$ times less likely. These ideas are developed in the Boltzmann distribution.

Logarithms

The logarithm undoes the exponential. If $y = e^x$, then $x = \ln y$, the natural logarithm of $y$. More generally, $\log_b y$ is the power to which the base $b$ must be raised to give $y$, so $\log_{10} 1000 = 3$ because $10^3 = 1000$. Since exponentials only produce positive numbers, logarithms are only defined for positive inputs.

Three rules follow directly from the rules for powers, and they are used constantly:

Before electronic calculators, these rules made logarithms the standard tool for multiplying large numbers: look up two logarithms, add them, and look up the number whose logarithm is the sum. Today their importance is different. Logarithms compress enormous ranges into manageable ones. The acidity of a solution spans more than fourteen orders of magnitude in hydrogen-ion concentration, and the pH scale, $\mathrm{pH} = -\log_{10}[\mathrm{H}^+]$, turns this into numbers between about 0 and 14. Sound intensity in decibels and earthquake magnitudes are logarithmic for the same reason. In statistical mechanics, the entropy $S = k_B \ln \Omega$ is the logarithm of a count of microstates that is itself astronomically large, and the free energy $F = -k_BT \ln Z$ is the logarithm of a sum of Boltzmann factors.

A logarithm grows without limit, but extraordinarily slowly: $\ln x$ reaches 10 only when $x$ exceeds 22 000, and 100 only when $x$ exceeds $10^{43}$. Each multiplication of $x$ by $e$ adds exactly one to $\ln x$, which is the logarithmic version of the exponential’s constant doubling time.

Logarithmic axes and laws from data

The rules of logarithms give a powerful way of discovering laws in data. Take the logarithm of both sides of a power law, $y = A x^n$:

This is a straight line, not in $y$ against $x$ but in $\log y$ against $\log x$, with slope $n$. Plotting data on log–log axes, where both axes are spaced logarithmically, turns every power law into a straight line, and the slope of the line reveals the power. Similarly, taking the logarithm of an exponential, $y = A e^{kx}$, gives $\ln y = \ln A + kx$, which is a straight line on semi-logarithmic axes with only the vertical axis logarithmic.

The figure below contains real data on the orbits of the eight planets: the size of each orbit (its semi-major axis, in units of the Earth’s) and the time it takes to go round the Sun (in years). On ordinary axes the points form a curve. Switch to log–log axes and they fall on a straight line with slope 1.5. This is Kepler’s third law, $T \propto a^{3/2}$, found by Kepler in 1619 by exactly this kind of search for a pattern in numbers, and explained seventy years later by Newton’s law of gravitation. The second data set is a simulated radioactive sample, whose counts fall on a straight line only on semi-log axes, as an exponential should. The third is the repulsive part of the Lennard-Jones potential, a pure power law with slope $-12$.

This technique is routine in research. The decay of a correlation function is plotted on semi-log axes to extract a relaxation time; the mean-squared displacement of a diffusing particle is plotted on log–log axes to check that it grows linearly with time; the error of a numerical method is plotted against its step size on log–log axes to measure its order of accuracy, as you will see in the chapter on derivatives.

Periodic functions

Many phenomena repeat: vibrations of atoms, alternating currents, the seasons, waves of light. The basic repeating functions are the sine and cosine. Imagine a point moving anticlockwise at a steady rate around a circle of radius 1 centred on the origin. If the angle it has turned through, measured from the positive horizontal axis, is $\theta$, then its height is $\sin\theta$ and its horizontal position is $\cos\theta$. As the point goes round and round, its height rises and falls forever, repeating every full turn of $2\pi$ radians (360°).

A general oscillation has an amplitude, a frequency and a phase:

Angles in science are measured in radians, defined so that an angle of $\theta$ radians cuts off an arc of length $\theta$ on a circle of radius 1. A full turn is $2\pi$ radians. Radians are not an arbitrary convention: as with the base $e$ for exponentials, they are the unit for which the derivative of $\sin\theta$ is exactly $\cos\theta$, with no awkward conversion factor.

Sinusoids describe any system that is pulled back towards equilibrium with a force proportional to its displacement, such as a mass on a spring or an atom vibrating in a molecule, as the chapter on motion shows. A remarkable theorem of Fourier says that essentially any periodic signal, however jagged, can be written as a sum of sines and cosines of different frequencies, which is why vibrational spectra, X-ray diffraction and much of signal processing are analysed in terms of sinusoids.

The Gaussian

The last family in this chapter is the bell-shaped Gaussian,

Choose the bell curve in the families figure: the shaded band covers $\mu \pm \sigma$ and the teal bar marks the full width at half maximum. The Gaussian is an exponential of a negative square, so it falls to zero much faster than an ordinary exponential. It is the shape of the distribution of velocities of molecules in a gas, of the spread of random measurement errors, of a spectral line broadened by thermal motion, and of the smearing functions used to turn atomic positions into smooth densities in machine-learning descriptors. The chapter on probability explains why it is so common: the sum of many small independent random effects always tends to a Gaussian.

Combining functions

New functions are built from old ones in a few standard ways. They can be added, subtracted, multiplied and divided. They can also be composed: feeding the output of one function into another, written $(f \circ g)(x) = f(g(x))$. The Gaussian is a composition: first square the input, then multiply by $-1/2\sigma^2$, then exponentiate. The Lennard-Jones energy of two atoms, $V(r) = 4\varepsilon[(\sigma/r)^{12} - (\sigma/r)^6]$, is a difference of two power laws that together produce a potential well with a steep repulsive wall and a gentle attractive tail. Recognising how a complicated function is assembled from simple pieces is usually the key to understanding its graph.

A function that undoes another is its inverse, written $f^{-1}$. The logarithm is the inverse of the exponential, and the square root (on non-negative numbers) is the inverse of squaring. The graph of an inverse is the mirror image of the original graph in the diagonal line $y = x$, since inverting a function simply swaps the roles of input and output. A function only has an inverse if it never gives the same output for two different inputs.

Finally, many quantities depend on more than one input. The potential energy of a molecule depends on the positions of all its atoms, and the loss of a neural network on millions of weights. A function of two variables, $f(x, y)$, can be drawn as a surface or as a contour map, and most of the mathematics of this chapter carries over. The chapter on gradients develops these functions of many variables.

References and further reading

Average rates of change

A car leaves at noon and has travelled 120 km by 2 pm. Its average speed was 60 km/h. That number is a ratio of two changes, the change in position divided by the change in time, and it says nothing about what happened in between: the car may have waited at traffic lights and then driven at 90 km/h on the motorway. The speedometer, by contrast, shows the speed at this instant. Making the idea of an instantaneous rate precise is the purpose of the derivative, and it is one of the most useful ideas in all of mathematics.

The average rate of change of a function $f$ between two inputs $x$ and $x + h$ is

Geometrically this is the slope of the straight line joining the two points $(x, f(x))$ and $(x + h, f(x + h))$ on the graph, called a secant line. For a straight-line function the secant slope is the same for every pair of points, which is what makes a line straight. For a curve it depends on which two points you choose.

The derivative as a limit

The way to get an instantaneous rate is to make the interval shorter and shorter. Take the secant line through $x$ and $x + h$ and slide the second point towards the first. As $h$ shrinks, the secant line swings round and settles towards a definite line that just grazes the curve at $x$: the tangent line. Its slope is the derivative.

There is a second way to see the same thing, and it is the more illuminating one. Zoom in on a smooth curve around a point. As the magnification increases, the curve looks less and less curved, until at high magnification it is indistinguishable from a straight line. The slope of that line is the derivative. A function that has this property at a point, of looking straight when magnified enough, is called differentiable there. Most functions in physics are differentiable almost everywhere; the exceptions are corners (like $|x|$ at zero), jumps and infinite spikes.

Try it: with a magnification of a few hundred the orange curve in the right panel lies exactly on the dashed tangent, and as $h$ is reduced the secant slope agrees with the true derivative to more and more decimal places. Both views describe one fact. Over a small enough region, any smooth function behaves linearly, $f(x + h) \approx f(x) + f'(x)\,h$, and the derivative is the coefficient of that linear behaviour.

Let us compute one derivative from the definition. For $f(x) = x^2$,

So the slope of the parabola $x^2$ at any point $x$ is $2x$: zero at the bottom, positive on the right, negative on the left, and steeper the further you go from the centre.

The derivative as a function

Since the derivative can be computed at every point, it is itself a function, $f'(x)$, and it has its own graph. Reading the two graphs together is one of the most important skills in applied mathematics. Where $f$ is rising, $f'$ is positive. Where $f$ is falling, $f'$ is negative. Where $f$ has a peak or a valley, its tangent is horizontal and $f'$ is zero. Such points are called stationary points, and finding them by solving $f'(x) = 0$ is how minima and maxima are located.

The derivative of the derivative is the second derivative, $f''(x) = \mathrm d^2 f/\mathrm dx^2$. It measures how the slope itself changes, which is the curvature of the graph. Where $f'' > 0$ the slope is increasing and the curve bends upwards like a cup; where $f'' < 0$ it bends downwards like a cap. This gives a test for stationary points: if $f'(x) = 0$ and $f''(x) > 0$, the point is a minimum; if $f''(x) < 0$, it is a maximum. In the figure, tick the box for the second derivative and check that it is positive at the minimum of $x^3 - 3x$ and negative at the maximum.

When the variable is time, the derivatives have familiar names. The derivative of position $x(t)$ is the velocity $v = \mathrm dx/\mathrm dt$, and the derivative of velocity is the acceleration $a = \mathrm dv/\mathrm dt = \mathrm d^2x/\mathrm dt^2$. Newton’s second law, $F = ma$, is therefore a statement about a second derivative, and the chapter on motion shows how it determines a trajectory.

The rules of differentiation

Computing every derivative from the limit definition would be tedious. Instead, a small number of rules, each proved once from the definition, handle almost every function you will meet. The derivatives of the basic functions are:

The power rule generalises the calculation done above for $x^2$. The rule for $e^x$ is the property that singles out the number $e$: of all exponentials $b^x$, only $e^x$ has a slope equal to its height everywhere. For a general base, the derivative of $b^x$ is $(\ln b)\,b^x$, and $\ln b = 1$ precisely when $b = e$.

Functions built from these pieces are handled by four combination rules:

The chain rule deserves special emphasis because it is everywhere. If $y$ depends on $u$ and $u$ depends on $x$, then a small change in $x$ changes $u$ by $u'(x)$ times as much, and that change in $u$ changes $y$ by $y'(u)$ times as much, so $\mathrm dy/\mathrm dx = (\mathrm dy/\mathrm du)(\mathrm du/\mathrm dx)$. For example, the derivative of the Gaussian $e^{-x^2/2}$ is $e^{-x^2/2}$ times the derivative of $-x^2/2$, that is $-x\,e^{-x^2/2}$. The training of every neural network is one enormous application of the chain rule, carried out link by link through the network, which is what backpropagation means.

Forces and equilibrium

In physics the most important derivative is the one connecting energy and force. For a particle moving along a line in a potential energy $V(x)$, the force is the negative slope of the potential:

Select the atom-pair energy in the derivative figure above. It is the Lennard-Jones potential, in reduced units $V(r) = 4(r^{-12} - r^{-6})$, which models the interaction of two noble-gas atoms. Using the power rule,

and setting $V'(r) = 0$ gives $r^6 = 2$, so the equilibrium separation is $r_0 = 2^{1/6} \approx 1.122$ (in units of the atomic diameter $\sigma$), where the energy is $V(r_0) = -1$, the depth of the well. To the left of $r_0$ the slope is negative and the force $-V'$ is positive, pushing the atoms apart; to the right it pulls them together. The second derivative at the minimum, $V''(r_0) = 72/2^{1/3} \approx 57.1$, is the stiffness of the bond, its effective spring constant. This number sets the frequency at which the atoms vibrate about their equilibrium spacing, which is how derivatives of energy surfaces connect to measured vibrational spectra. The same reasoning in many dimensions gives forces as gradients, developed in the chapter on gradients and used throughout molecular dynamics.

Optimisation

Finding the input that makes a function as small or as large as possible is called optimisation, and it is the central computational task of both chemistry (finding stable structures) and machine learning (finding the best parameters). The calculus answer is to solve $f'(x) = 0$ and check the sign of $f''$. This works beautifully when the derivative can be solved by hand, as for the Lennard-Jones minimum. When it cannot, the derivative still tells you which way to go: if $f'(x) > 0$, decreasing $x$ lowers $f$, and if $f'(x) < 0$, increasing $x$ does. Repeatedly stepping a little way in the downhill direction,

with a small step size $\eta$, walks down to a minimum. This is gradient descent, the optimisation method behind nearly all of deep learning, studied in detail in its own chapter. A faster method for smooth functions is Newton’s method, which uses the second derivative as well: it approximates the function near $x$ by a parabola and jumps straight to the parabola’s minimum, $x_{\text{new}} = x - f'(x)/f''(x)$. Near a minimum it roughly doubles the number of correct digits at every step.

Taylor series

The idea that a smooth function looks straight up close can be pushed further. A straight line matches a function’s value and slope at a point. A parabola can also match its curvature, a cubic the rate of change of curvature, and so on. Matching all the derivatives at a point $a$ gives the Taylor series:

To see why the coefficients must be $f^{(n)}(a)/n!$, differentiate the polynomial $n$ times and set $x = a$: every term except the $n$-th vanishes, and the $n$-th gives $n!$ times its coefficient, which must equal $f^{(n)}(a)$.

Some Taylor series are worth knowing by heart, all expanded about $a = 0$:

The first two converge for every $x$: with enough terms the polynomial matches the function on any finite interval, as the figure shows. The third converges only for $-1 < x \le 1$. The distance from the expansion point within which a series converges is its radius of convergence, and it is limited by the nearest point where the function misbehaves, here the logarithm’s singularity at $x = -1$.

Truncated Taylor series are behind a large fraction of physical approximations. The small-angle approximation $\sin\theta \approx \theta$ is the first-order series, and it turns the pendulum’s nonlinear equation into the simple harmonic oscillator. The harmonic approximation of a chemical bond expands the potential energy about its minimum, where $V'(r_0) = 0$:

a parabola with spring constant $k = V''(r_0)$. Choose the bond energy in the figure and set the order to 2: the parabola hugs the true curve near the minimum but misses the soft attractive tail and the steep wall. Adding the cubic term captures the asymmetry, which is what makes real crystals expand when heated. Normal-mode analysis of molecules, phonons in crystals and vibrational spectroscopy are all built on the second-order Taylor expansion of the energy.

Derivatives on a computer

Simulations and data analysis often need derivatives of functions known only as numbers, or as code too complicated to differentiate by hand. The definition suggests an approximation: choose a small but finite $h$ and compute the forward difference $[f(x + h) - f(x)]/h$. Taylor’s theorem tells us how accurate this is. Since $f(x + h) = f(x) + f'(x) h + \tfrac12 f''(x) h^2 + \dots$, the forward difference equals $f'(x) + \tfrac12 f''(x) h + \dots$, so its error is proportional to $h$. A cleverer formula, the central difference $[f(x + h) - f(x - h)]/2h$, cancels the $h^2$ terms of the two expansions and has an error proportional to $h^2$, which is far smaller for small $h$.

It would seem, then, that the best choice is the smallest possible $h$. The figure below shows why this is wrong. It plots the actual error of both formulas for the derivative of $\sin x$ at $x = 1$, computed in the same floating-point arithmetic a real program uses, against $h$ on log–log axes.

For large $h$ the errors fall along straight lines of slope 1 (forward) and 2 (central): this is the truncation error predicted by the Taylor series. But below a certain $h$ the curves turn round and the error grows again, becoming erratic. The cause is rounding. A computer stores about 16 significant decimal digits in double precision (the chapter on numbers in a computer explains why). When $h$ is tiny, $f(x + h)$ and $f(x)$ agree in almost all their digits, and subtracting them leaves only the few digits that differ, most of which are rounding noise. Dividing by the tiny $h$ then magnifies that noise. The round-off error grows like $\epsilon/h$, where $\epsilon \approx 2.2 \times 10^{-16}$ is the machine precision.

The best step balances the two errors. For the forward difference the optimum is near $h \approx \sqrt{\epsilon} \approx 10^{-8}$, giving an error around $10^{-8}$; for the central difference it is near $h \approx \epsilon^{1/3} \approx 10^{-5}$, giving an error of $10^{-10}$ or better. Tick the single-precision box, as used on many GPUs, and the optimum moves to much larger $h$ with an error of only around $10^{-4}$ to $10^{-7}$. Finite differences therefore never reach full machine precision, and they need two or more function evaluations per input variable, which is prohibitive for a function of a million variables.

Modern machine-learning software sidesteps both problems with automatic differentiation. Rather than approximating, it applies the chain rule exactly to every elementary operation the program performs, carrying derivatives along with values. The result is exact up to rounding, and in its reverse mode the full gradient of a function of any number of variables costs only a small multiple of one function evaluation. It is how neural networks are trained and how differentiable molecular dynamics computes the derivatives of a whole simulation.

References and further reading

Accumulation

The derivative answers the question “how fast is this changing?”. The integral answers the reverse question: “if I know how fast something changes at every moment, how much has accumulated in total?”. If a car travels at a constant 60 km/h for 2 hours, it covers $60 \times 2 = 120$ km. On a graph of speed against time, that product is the area of a rectangle 2 hours wide and 60 km/h tall. The units confirm it: (km/h) × h = km. Area under a rate curve is accumulated quantity.

When the speed varies, the region under the graph is no longer a rectangle, but the principle is the same. Over a short enough time interval the speed is nearly constant, so the distance covered in that interval is nearly speed × duration, the area of a thin vertical strip. Adding up the strips gives the total distance. The integral is the exact value this sum approaches as the strips become infinitely thin.

The same idea appears whenever a total is built from a varying rate or density. The mass of a rod whose density varies along its length is the integral of density over length. The charge that has flowed through a wire is the integral of the current over time. The amount of product formed in a chemical reaction is the integral of the reaction rate. The work done by a varying force is the integral of force over distance. And, as later chapters show, every average and every probability in statistical mechanics is an integral.

Riemann sums and the definite integral

To make this precise, divide the interval from $a$ to $b$ into $n$ strips of equal width $\Delta x = (b - a)/n$. In each strip pick a sample point $x_i$, for instance its left edge, and form the sum of rectangle areas:

The figure below draws the sum. With a handful of strips the rectangles overshoot in some places and undershoot in others, and the total is visibly wrong. Press “Double n” a few times and the staircase hugs the curve ever more closely while the error shrinks. For any continuous function the sum converges to the same limit whichever sample points are used, which is why the integral is well defined. The choice still matters for a finite number of strips, as comparing the left-edge, right-edge, midpoint and trapezoid rules shows.

Notice that where the function is negative, the strips hang below the axis and contribute negative amounts. The definite integral is a signed area: area above the axis counts positively and area below counts negatively. For $\sin x$ over a full period from $0$ to $2\pi$ the positive and negative lobes cancel exactly and the integral is zero. Physically this is the right behaviour. If a velocity is negative, the object is moving backwards and the distance it covers should be subtracted from its net displacement.

The fundamental theorem of calculus

Computing integrals as limits of sums is laborious. The great discovery of Newton and Leibniz in the 1660s and 1670s was that integration is the reverse of differentiation, which turns many integrals into a lookup in a table of derivatives.

Fix the lower limit $a$ and let the upper limit vary. The accumulated area up to $x$ is then a function of $x$,

How fast does $A$ grow as $x$ increases? Moving the upper limit from $x$ to $x + h$ adds a thin strip of width $h$ and height approximately $f(x)$, so $A(x + h) - A(x) \approx f(x)\,h$. Dividing by $h$ and letting $h \to 0$ gives the theorem:

The figure below makes the theorem visible. The upper panel shades the area under a rate curve up to a moving point $x$. The lower panel plots that accumulated area $A(x)$ as a function of $x$. Drag the point and compare the teal tangent to the lower curve with the height of the upper curve: they are always equal. Where the rate is large, $A$ climbs steeply; where the rate is zero, $A$ is momentarily flat; where the rate is negative, $A$ goes down.

With the theorem, integrals become calculations. To find $\int_0^3 x^2\,\mathrm dx$, look for a function whose derivative is $x^2$: the power rule run backwards gives $F(x) = x^3/3$. The integral is $F(3) - F(0) = 9$. Every row of the table of derivatives, read from right to left, is an integral:

Techniques and the Gaussian integral

Two techniques extend the table a long way. Substitution is the chain rule in reverse: if the integrand has the form $f(g(x))\,g'(x)$, then it is the derivative of $F(g(x))$, so $\int f(g(x)) g'(x)\,\mathrm dx = F(g(x))$. For example, $\int x\,e^{-x^2/2}\,\mathrm dx = -e^{-x^2/2}$, because the factor $x$ is (up to sign) the derivative of the exponent. Integration by parts is the product rule in reverse: $\int u\,v'\,\mathrm dx = uv - \int u'\,v\,\mathrm dx$. It moves a derivative from one factor to the other and is used, among many other places, to derive the virial expression for pressure in molecular simulation.

Not every integral has an antiderivative that can be written with elementary functions. The most important example in all of science is the Gaussian integral. The function $e^{-x^2}$ has no elementary antiderivative, yet its integral over the whole line has an exact value, found by a beautiful trick: square the integral, write the square as a double integral over the plane, and switch to polar coordinates, where the extra factor of $r$ makes the integrand the derivative of $-\tfrac12 e^{-r^2}$.

The normalisation of the normal distribution, the partition function of an ideal gas, the momentum integrals that give the Maxwell–Boltzmann distribution and the zero-point properties of a harmonic oscillator all reduce to this one result. In the strips figure, choose $e^{-x^2}$: the area between $-3$ and $3$ is already $1.7724$, within $0.003\,\%$ of $\sqrt\pi = 1.77245$, because the tails beyond $|x| = 3$ are tiny.

Integrals in physics

Several integrals recur so often in the physical sciences that they are worth recognising on sight.

The work done by a force $F(x)$ that moves an object from $a$ to $b$ is $W = \int_a^b F(x)\,\mathrm dx$. When the force comes from a potential, $F = -\mathrm dV/\mathrm dx$, the fundamental theorem gives $W = -[V(b) - V(a)]$: the work done by the force equals the drop in potential energy. This is the integral form of the statement that force is the slope of energy, and it is why energy is conserved for such forces. Choose the spring force in the accumulation figure: the accumulated work is a downward parabola, minus the familiar $\tfrac12 kx^2$.

The average of a function over an interval is its integral divided by the length of the interval, $\bar f = \frac{1}{b - a}\int_a^b f(x)\,\mathrm dx$: the height of the rectangle with the same area. A time average from a simulation is exactly this, with the integral replaced by a sum over saved frames.

A probability density $p(x)$ describes a continuous random quantity: the probability that it lies between $a$ and $b$ is $\int_a^b p(x)\,\mathrm dx$, and the total integral is 1. The average of any quantity $g(x)$ is then $\langle g \rangle = \int g(x)\,p(x)\,\mathrm dx$. The chapter on probability develops this, and statistical mechanics is built on it: the average energy of a system is an integral of energy weighted by the Boltzmann probability over all positions and momenta.

Integrals extend naturally to several dimensions. The volume under a surface $f(x, y)$ is the double integral $\iint f\,\mathrm dx\,\mathrm dy$, the limit of a sum of thin columns, and the mass of a three-dimensional object is the triple integral of its density over its volume. The partition function of $N$ classical atoms is an integral over all $3N$ position coordinates and $3N$ momentum coordinates, a $6N$-dimensional integral. Computing such integrals is the central difficulty of statistical mechanics, and it is the reason for the last section of this topic.

Numerical integration

When an antiderivative is unavailable, or the integrand is known only as data, integrals are computed numerically by a quadrature rule: a weighted sum of function values. The rules in the strips figure are the simplest. Their accuracy can be read off from Taylor’s theorem, just as for finite differences. With $n$ strips of width $h$, the left- and right-edge rules have errors proportional to $h$; the midpoint and trapezoid rules, which are symmetric within each strip, have errors proportional to $h^2$; and Simpson’s rule, which fits a parabola through each pair of strips, has error proportional to $h^4$. Doubling $n$ therefore halves the error of the left-edge rule, quarters that of the trapezoid rule and divides that of Simpson’s rule by sixteen. You can check the first two in the figure: compare the errors before and after pressing “Double n”.

The same ideas, applied to the equations of motion, give the integrators of molecular dynamics. Advancing a trajectory by one time step means integrating the acceleration to update the velocity and integrating the velocity to update the position. The velocity Verlet algorithm used in nearly every MD code is a carefully symmetrised version of the midpoint idea, as the chapter on integrators explains.

Monte Carlo integration

Grid-based quadrature has a fatal weakness in many dimensions. A grid with 10 points along each axis needs $10^d$ function evaluations in $d$ dimensions: 100 in two dimensions, but $10^{30}$ in thirty, far beyond any computer. Phase-space integrals over hundreds of atoms are hopeless on a grid. This exponential cost is called the curse of dimensionality.

The way out is to sample at random. To estimate the area of a region, enclose it in a box of known area, scatter $N$ random points uniformly in the box, and count the fraction that fall inside. That fraction times the box area estimates the region’s area. The figure below estimates π this way: the quarter-circle of radius 1 has area $\pi/4$, so four times the fraction of random points inside it approaches π.

The error falls as $1/\sqrt N$: a hundred times more darts buys one more correct digit. That is slow compared with Simpson’s rule in one dimension. The crucial point is that the $1/\sqrt N$ law does not depend on the dimension at all. The same number of random samples gives the same statistical error whether the integral is over two variables or six thousand. This is why the integrals of statistical mechanics are evaluated by sampling, either with Monte Carlo methods that draw configurations at random with the right probabilities, or with molecular dynamics, which generates them by following the motion. The error bar on every simulated average is ultimately this $1/\sqrt N$ law, with $N$ the number of effectively independent samples, as the chapter on probability explains.

References and further reading

Vectors

Some quantities are fully described by a single number: a temperature, a mass, an energy. These are scalars. Others need a direction as well as a size. A displacement of 3 m means nothing until you say which way; a force of 10 N pushes in some particular direction; the velocity of an atom has a speed and a heading. Such quantities are vectors.

A vector can be pictured as an arrow. Its length is the vector’s size, or magnitude, and its direction is where it points. Two arrows with the same length and direction represent the same vector wherever they are drawn. The same vector can also be written as a list of numbers. Draw the arrow starting at the origin of a coordinate system; the coordinates of its tip are the vector’s components. In two dimensions,

and in three dimensions a third component $v_z$ is added. Vectors are written in bold, $\mathbf v$, or with an arrow, $\vec v$.

The two pictures, arrow and list, are equally important. The arrow is how we think about vectors geometrically; the list is how we compute with them, and how a computer stores them. The power of linear algebra comes from moving freely between the two.

The list picture also frees vectors from the three dimensions of physical space. The positions of all $N$ atoms in a molecule form a single vector with $3N$ components, a point in a $3N$-dimensional configuration space. A material can be described by a vector of features: its composition, density, band gap and so on. An image of $28 \times 28$ pixels is a vector with 784 components. We cannot draw arrows in 784 dimensions, but every rule in this chapter works in any number of dimensions, and the two-dimensional pictures are faithful guides to what happens there.

Adding and scaling vectors

Two operations define what a vector is. Addition is done component by component, $\mathbf a + \mathbf b = (a_x + b_x,\, a_y + b_y)$, and geometrically it means placing the arrows tip to tail: walk along $\mathbf a$, then along $\mathbf b$, and the sum is the arrow from start to finish. Because it does not matter which arrow is walked first, the two routes form a parallelogram whose diagonal is the sum. Forces combine this way: the net force on an atom is the vector sum of the forces from all its neighbours.

Scalar multiplication multiplies every component by the same number, $c\,\mathbf a = (c\,a_x,\, c\,a_y)$. Geometrically it stretches the arrow by the factor $c$, and a negative $c$ also reverses its direction. Subtraction combines the two: $\mathbf a - \mathbf b = \mathbf a + (-1)\mathbf b$, the arrow from the tip of $\mathbf b$ to the tip of $\mathbf a$. The vector from atom $i$ to atom $j$ is $\mathbf r_j - \mathbf r_i$, and its length is the interatomic distance.

A sum of scaled vectors, $c_1\mathbf a + c_2\mathbf b$, is called a linear combination. Every vector in the plane can be written as a linear combination of the two unit vectors $\hat{\mathbf i} = (1, 0)$ and $\hat{\mathbf j} = (0, 1)$: the components of a vector are simply the coefficients, $\mathbf v = v_x\hat{\mathbf i} + v_y\hat{\mathbf j}$. A set of vectors from which every vector can be built in exactly one way is a basis, and the number of vectors in a basis is the dimension of the space. The pair $\hat{\mathbf i}, \hat{\mathbf j}$ is the standard basis, but any two vectors that do not point along the same line also form a basis of the plane. Crystallographers use the three lattice vectors of a unit cell as a basis, so that atomic positions are given as fractional coordinates along the cell edges.

Length and the dot product

The length of a vector follows from Pythagoras’ theorem: $|\mathbf v| = \sqrt{v_x^2 + v_y^2}$, with a third term $v_z^2$ in three dimensions and $n$ terms in $n$ dimensions. A vector of length 1 is a unit vector, and any non-zero vector can be turned into one by dividing it by its length.

The most useful way of combining two vectors into a number is the dot product:

The equality of the two formulas, one algebraic and one geometric, is the law of cosines in disguise. It means the angle between two vectors can be computed from their components: $\cos\theta = \mathbf a \cdot \mathbf b / (|\mathbf a||\mathbf b|)$. Bond angles in a molecule are computed exactly this way from atomic coordinates.

The geometric meaning is a projection. $\mathbf a \cdot \mathbf b / |\mathbf a|$ is the length of the shadow that $\mathbf b$ casts onto the line through $\mathbf a$. In the figure, select the dot product and rotate $\mathbf b$ around: the shadow is longest when the vectors are aligned, vanishes when they are perpendicular, and becomes negative when $\mathbf b$ points backwards. Physics uses this constantly. The work done by a force $\mathbf F$ over a displacement $\mathbf d$ is $\mathbf F \cdot \mathbf d$, since only the component of the force along the motion does work. Machine learning uses it as a measure of similarity: two feature vectors with a large positive dot product point in similar directions, and the cosine similarity $\cos\theta$ compares them regardless of their lengths. The attention mechanism of transformers decides how much one word should attend to another by taking dot products of vectors that represent them.

Matrices as transformations

A matrix is a rectangular array of numbers. Matrices can be thought of as tables of data, but their real meaning is as linear transformations: functions that take a vector in and give a vector out, while keeping straight lines straight, keeping the origin fixed and keeping evenly spaced grid lines evenly spaced. Rotations, stretches, reflections and shears are all linear transformations.

The key observation is that a linear transformation is completely determined by what it does to the basis vectors. Since every vector is $\mathbf v = x\hat{\mathbf i} + y\hat{\mathbf j}$, and a linear transformation respects addition and scaling, the image of $\mathbf v$ must be $x$ times the image of $\hat{\mathbf i}$ plus $y$ times the image of $\hat{\mathbf j}$. Write the images of $\hat{\mathbf i}$ and $\hat{\mathbf j}$ as the two columns of a matrix, and the rule for multiplying a matrix by a vector follows:

The figure below shows a matrix acting on the entire plane at once. The faint grey grid is the original plane; the blue grid is where every grid line goes. Drag the green and red arrows, which are the two columns of the matrix, and the whole grid follows. Try the presets and press play to watch the plane morph from its original state into the transformed one.

Several transformations are worth recognising by their matrices. The identity $\begin{pmatrix}1&0\\0&1\end{pmatrix}$ does nothing. A rotation by an angle $\theta$ is $\begin{pmatrix}\cos\theta&-\sin\theta\\\sin\theta&\cos\theta\end{pmatrix}$, since it sends $\hat{\mathbf i}$ to $(\cos\theta, \sin\theta)$. A diagonal matrix $\begin{pmatrix}s_x&0\\0&s_y\end{pmatrix}$ stretches each axis separately. The shear $\begin{pmatrix}1&1\\0&1\end{pmatrix}$ slides each horizontal line sideways by an amount proportional to its height, as when a deck of cards is pushed; it is the transformation of a simulation box under shear flow. And a matrix whose two columns lie on the same line squashes the whole plane onto that line.

Composition and matrix multiplication

Applying one transformation and then another is again a linear transformation, so it has a matrix of its own. That matrix is the product of the two. If $\mathbf A$ is applied first and $\mathbf B$ second, the combined transformation is $\mathbf B\mathbf A$ (read from right to left, like function composition $f(g(x))$). To find its columns, apply $\mathbf B$ to each column of $\mathbf A$; the familiar rule “row times column” for each entry follows directly.

Two consequences are important. Matrix multiplication is not commutative: in general $\mathbf A\mathbf B \ne \mathbf B\mathbf A$, because rotating and then shearing is not the same as shearing and then rotating. Try it with a book on a table. And a transformation that can be undone has an inverse matrix $\mathbf A^{-1}$ with $\mathbf A^{-1}\mathbf A = \mathbf I$, the identity. For a $2 \times 2$ matrix,

which exists only when $ad - bc \ne 0$. That quantity has a geometric meaning of its own.

The determinant

Watch the gold square in the transformation figure. It starts as the unit square with area 1, and the matrix turns it into a parallelogram. Every other region of the plane is scaled by the same factor, since the transformed grid cells are all identical parallelograms. That factor is the determinant:

In the figure, the rotations have determinant 1 (areas are preserved), the stretch has determinant 1 because it doubles one axis and halves the other, the reflection has $-1$, and “squash flat” has 0: the parallelogram collapses to a line. In three dimensions the determinant of the matrix whose columns are the lattice vectors of a crystal is the volume of the unit cell, which is how simulation codes compute the volume of a triclinic box. In statistical mechanics, determinants appear whenever variables are changed in an integral, as the factor that accounts for how the change stretches volume elements.

Linear systems

Many problems ask the reverse question: given the output $\mathbf b$, which input $\mathbf x$ produced it? That is, solve $\mathbf M\mathbf x = \mathbf b$. Written out, this is a set of simultaneous linear equations. If $\det\mathbf M \ne 0$ the transformation can be undone and there is exactly one solution, $\mathbf x = \mathbf M^{-1}\mathbf b$. If $\det\mathbf M = 0$, the plane has been squashed, so most $\mathbf b$ cannot be reached at all and those that can are reached by infinitely many $\mathbf x$.

In practice, inverses are rarely computed. Computers solve linear systems by Gaussian elimination, subtracting multiples of equations from one another until the system becomes triangular and can be solved from the bottom up; in matrix language this is a factorisation $\mathbf M = \mathbf L\mathbf U$ into lower- and upper-triangular pieces. Fitting a straight line or a polynomial to data by least squares, the first algorithm of machine learning, comes down to solving a linear system, the normal equations $\mathbf X^{\mathsf T}\mathbf X\,\mathbf w = \mathbf X^{\mathsf T}\mathbf y$. Gaussian processes, Kriging and kernel methods require solving systems with thousands of unknowns, and the Cholesky factorisation used for them is a specialised version of elimination for symmetric matrices.

Eigenvectors and eigenvalues

Most vectors are knocked off their line by a transformation: in the figure, the blue vector $\mathbf v$ and its image $\mathbf M\mathbf v$ usually point in different directions. But for many matrices there are special directions that are only stretched or shrunk, never turned. A vector $\mathbf u$ along such a direction satisfies

Tick “show eigenvectors” in the transformation figure. For the symmetric matrix $\begin{pmatrix}2&1\\1&2\end{pmatrix}$, the two dashed purple lines along the diagonals are the eigen-directions: vectors on the line $y = x$ are stretched by 3, vectors on the line $y = -x$ are left alone (eigenvalue 1), and every other vector is a mixture and gets turned. The shear has only one eigen-direction, the horizontal axis. A rotation by 90° has none at all: every direction is turned, and the eigenvalues are complex numbers. Tick “what happens to a circle” as well, and you will see the unit circle become an ellipse whose axes, for a symmetric matrix, lie exactly along the eigenvectors, with semi-axes equal to the eigenvalues.

Symmetric matrices, those equal to their own transpose ($M_{ij} = M_{ji}$), are especially well behaved, and they are the ones that matter most in science. Their eigenvalues are always real, and their eigenvectors can always be chosen to be mutually perpendicular. Using the eigenvectors as a new basis, the matrix becomes diagonal: the transformation is just a stretch along each eigen-direction. This spectral decomposition, $\mathbf M = \mathbf Q\boldsymbol\Lambda\mathbf Q^{\mathsf T}$ with $\mathbf Q$ the matrix of eigenvectors and $\boldsymbol\Lambda$ the diagonal matrix of eigenvalues, appears across the sciences:

The normal modes of a molecule are the eigenvectors of its Hessian, the matrix of second derivatives of the energy with respect to atomic positions (mass-weighted), and the eigenvalues give the squared vibrational frequencies. A negative eigenvalue signals a saddle point rather than a minimum, which is how transition states are recognised, as the chapter on gradients shows.

Principal component analysis finds the eigenvectors of the covariance matrix of a data set. The first eigenvector is the direction in which the data vary most, and projecting onto the first few eigenvectors gives the best low-dimensional summary of high-dimensional data, whether those are atomic trajectories or materials properties.

In quantum mechanics the allowed energies of a system are the eigenvalues of its Hamiltonian matrix, and the stationary states are its eigenvectors. Every electronic-structure calculation, from a Hückel model of benzene to density-functional theory, is at its heart a large symmetric eigenvalue problem.

Matrices in machine learning and simulation

A single layer of a neural network takes an input vector $\mathbf x$, multiplies it by a weight matrix $\mathbf W$, adds a bias vector and applies a nonlinear function to each component: $\mathbf h = \sigma(\mathbf W\mathbf x + \mathbf b)$. The matrix rotates, stretches and shears the input space; the nonlinearity bends it. Stacking such layers lets a network reshape space until classes that were tangled become separable, as the chapter on how networks bend space shows with the same kind of grid picture used above. Training a network means adjusting the entries of these matrices, and running one means doing matrix multiplications, which is why graphics processors, built to multiply matrices quickly, became the engines of deep learning.

In molecular simulation, the simulation box is described by a $3 \times 3$ cell matrix whose columns are the box vectors; converting between Cartesian and fractional coordinates is multiplication by this matrix or its inverse, and the box volume is its determinant. Rotation matrices orient molecules, the stress and strain of a material are $3 \times 3$ symmetric matrices whose eigenvectors are the principal axes, and the moment of inertia tensor, diagonalised, gives a molecule’s principal axes of rotation. In all of these the same few ideas recur: a matrix is a transformation, its columns are the images of the basis, its determinant scales volume, and its eigenvectors are the directions it leaves unturned.

References and further reading

Functions of several variables

The energy of a diatomic molecule depends on one number, the bond length. The energy of a water molecule depends on two bond lengths and an angle, and the energy of a protein on the positions of thousands of atoms. The loss of a neural network depends on every one of its weights, often millions of them. Almost every function that matters in simulation and machine learning has many inputs, and this chapter extends the calculus of derivatives to them.

A function of two variables, $f(x, y)$, assigns a number to each point of the plane. It can be drawn as a surface hovering over the plane, with height $f(x, y)$ above each point, which is why such functions are often called landscapes. A surface is hard to read on a flat page, so the more useful picture is a contour map, like the maps hikers use. Each contour line joins points of equal height, $f(x, y) = \text{constant}$. Where contours crowd together the ground is steep; where they are far apart it is nearly flat. Closed rings of contours surround summits and valley floors, and a figure-of-eight pattern marks a mountain pass.

The landscapes in the figure are chosen because they have the features found on real potential energy surfaces. The bowl has a single minimum. The long curved valley is gently sloped along its floor and very steep across it, a shape that defeats simple optimisation methods. The landscape with two valleys has two stable states separated by a pass, the simplest model of a chemical reaction or a structural transition. The three-basin landscape is a sum of Gaussian wells in the spirit of the Müller–Brown surface, a standard two-dimensional test problem for reaction-path methods.

Partial derivatives

The ordinary derivative measures how a function changes when its single input changes. With two inputs, there are many directions in which to move. The simplest are along the coordinate axes. Hold $y$ fixed and let only $x$ vary: the function becomes an ordinary function of one variable, whose graph is a slice through the landscape along the $x$ direction. Its slope is the partial derivative with respect to $x$:

For example, if $f(x, y) = x^2 y + 3y^2$, then treating $y$ as a constant gives $\partial f/\partial x = 2xy$, and treating $x$ as a constant gives $\partial f/\partial y = x^2 + 6y$. Tick “slices” in the figure: the two small panels show the landscape cut along $x$ and along $y$ through the walker, each with its tangent line. The slopes of those tangents are the two partial derivatives.

Partial derivatives can themselves be differentiated, giving second partial derivatives such as $\partial^2 f/\partial x^2$, $\partial^2 f/\partial y^2$ and the mixed derivative $\partial^2 f/\partial x\,\partial y$. For any smooth function the mixed derivatives do not depend on the order of differentiation, $\partial^2 f/\partial x\,\partial y = \partial^2 f/\partial y\,\partial x$. This innocent-looking fact, Schwarz’s theorem, is the origin of the Maxwell relations of thermodynamics, which relate measurable quantities such as $(\partial S/\partial V)_T = (\partial P/\partial T)_V$.

The gradient

Collect the partial derivatives into a vector and something remarkable appears. The gradient of $f$ is

The directional derivative formula comes from the linear approximation. Near a point, a smooth landscape looks like a tilted plane: moving by a small displacement $(\Delta x, \Delta y)$ changes $f$ by approximately $\frac{\partial f}{\partial x}\Delta x + \frac{\partial f}{\partial y}\Delta y = \nabla f \cdot \Delta\mathbf r$. Walking a distance $s$ in the direction $\mathbf u$ means $\Delta\mathbf r = s\,\mathbf u$, so the slope is $\nabla f \cdot \mathbf u$.

Three geometric facts follow from $\nabla f \cdot \mathbf u = |\nabla f|\cos\theta$, and you can verify each in the figure by turning the violet compass arrow. The slope is largest when $\mathbf u$ points along $\nabla f$ ($\theta = 0$): the gradient points in the direction of steepest ascent, and its length $|\nabla f|$ is that steepest slope. The slope is most negative in the opposite direction, $-\nabla f$, the direction of steepest descent. And the slope is zero when $\mathbf u$ is perpendicular to $\nabla f$: walking at right angles to the gradient keeps you at constant height, so the gradient is perpendicular to the contour lines. Tick “arrow field” to see $-\nabla f$ drawn across the whole map: the arrows cut every contour at right angles and grow longer where the contours crowd together.

Forces as gradients

For a system of atoms with potential energy $V(\mathbf r_1, \dots, \mathbf r_N)$, the force on atom $i$ is the negative gradient of the energy with respect to that atom’s position:

This is the many-dimensional version of $F = -\mathrm dV/\mathrm dx$, and it is the single most important formula in molecular simulation: every MD step computes it. For a pair potential that depends only on the distance $r_{ij} = |\mathbf r_j - \mathbf r_i|$ between two atoms, the chain rule gives $\mathbf F_i = V'(r_{ij})\,(\mathbf r_j - \mathbf r_i)/r_{ij}$, a force along the line joining the atoms, and the force on atom $j$ is exactly opposite: Newton’s third law follows automatically from the energy depending only on distances. It is also why machine-learned interatomic potentials are built to predict an energy and then obtain forces by differentiating it, which guarantees that the forces are consistent with energy conservation, as the chapter on MLIPs explains.

Stationary points and the Hessian

Points where the gradient vanishes are stationary points. In one dimension they are minima or maxima. In two or more dimensions there is a third possibility, the saddle point: a minimum along one direction and a maximum along another, like the middle of a horse’s saddle or the highest point of a mountain pass on the lowest route between two valleys. Drag the walker to the middle of the two-valley landscape, at the origin: the gradient arrows disappear, yet the point is clearly neither a summit nor a valley floor.

To tell these apart one needs the second derivatives, arranged in a symmetric matrix called the Hessian:

The second formula is the Taylor series in several variables. At a stationary point the gradient term vanishes and the landscape is locally a quadratic bowl, dome or saddle described by $\mathbf H$. Because $\mathbf H$ is symmetric, it has real eigenvalues and perpendicular eigenvectors (the chapter on matrices explains why). Tick “curvature” in the figure: the walker now carries two perpendicular bars along the Hessian’s eigenvectors, gold for positive curvature and rose for negative. At the bottom of either valley of the two-valley landscape both bars are gold, with eigenvalues 8 and 2.4. At the pass the bar across the valley is rose with eigenvalue $-4$ and the bar along it is gold: one downhill direction, the signature of a saddle.

In chemistry these ideas have direct physical meaning. Minima of the potential energy surface are stable or metastable structures. The eigenvalues of the mass-weighted Hessian at a minimum are the squares of the vibrational frequencies, and the eigenvectors are the normal modes. A first-order saddle point, with exactly one negative eigenvalue, is a transition state: the highest point on the lowest-energy path between two minima. Its energy above the reactant minimum is the activation barrier that enters the Arrhenius law, and the eigenvector of its negative eigenvalue points along the reaction coordinate.

Following the gradient downhill

A ball released on a landscape, with enough friction that it never builds up speed, follows the steepest-descent direction $-\nabla f$ at every moment. Its path is a steepest-descent path, and it ends at a minimum (or, if released exactly on a ridge, at a saddle). Press “roll downhill” in the figure from several starting points. Every start on one side of the pass ends in one valley and every start on the other side ends in the other: the set of starting points that drain into a minimum is that minimum’s basin of attraction, and the boundary between basins runs through the saddle.

Turning this into an algorithm gives gradient descent: repeat $\mathbf r \leftarrow \mathbf r - \eta\,\nabla f(\mathbf r)$ with a small step $\eta$. It is the method used to minimise the energy of a structure before a simulation, and, with the energy replaced by a loss function, the method by which every neural network is trained (the chapter on gradient descent). Its weakness shows on the long curved valley. There the gradient points mostly across the valley rather than along it, because the walls are much steeper than the floor, so steepest descent zigzags from wall to wall and crawls along the floor. The cause is visible in the Hessian: its eigenvalues differ by a large factor, a situation called ill-conditioning. Better methods either rescale the directions (Newton’s method uses $\mathbf H^{-1}\nabla f$, and quasi-Newton methods such as L-BFGS build up an approximation to it) or add momentum, which averages out the zigzag.

The steepest-descent path from a transition state down into the two minima it connects is the minimum energy path, also called the intrinsic reaction coordinate. Methods such as the nudged elastic band find it by relaxing a chain of images of the system towards the path, and it provides the reaction coordinate for many rate calculations.

The multivariable chain rule

Suppose a quantity $f$ depends on $x$ and $y$, and both depend on time: $f(x(t), y(t))$. The rate of change of $f$ along the path collects a contribution from each input:

The rate of change is the gradient dotted with the velocity. Applied to the potential energy of a moving system, $\mathrm dV/\mathrm dt = \nabla V\cdot\mathbf v = -\mathbf F\cdot\mathbf v$, the negative of the power delivered by the forces, which is half of the proof that total energy is conserved in classical mechanics.

When many variables depend on many others, the partial derivatives of the outputs with respect to the inputs form a matrix, the Jacobian, and the chain rule becomes multiplication of Jacobians. A neural network is a long chain of such dependencies, layer by layer, and the gradient of its loss with respect to the weights of an early layer is a product of Jacobians of all the later layers. Computing that product efficiently from the output end backwards is backpropagation, developed in its own chapter. Every gradient used in modern machine learning, and every force in a machine-learned potential, is computed this way.

References and further reading

Chance and frequency

Toss a coin and you cannot say whether it will land heads. Toss it ten thousand times and you can say with great confidence that close to half the tosses will be heads. This combination, unpredictability in each case and regularity in the aggregate, is what probability theory describes, and it is everywhere in science. A single molecule in a gas moves erratically, but the pressure of $10^{23}$ molecules is steady. A single radioactive nucleus decays at an unpredictable moment, but a sample decays with a precise half-life. A single run of a stochastic training algorithm gives a slightly different network each time, but the distribution of results is reproducible.

The probability of an outcome is the fraction of times it occurs in a long sequence of identical, independent trials. This frequency interpretation is made precise by the law of large numbers: as the number of trials grows, the observed fraction converges to the probability. The convergence is not smooth. The figure below tracks the running fraction of heads, or of sixes, as trials accumulate. Early on it swings wildly; later it settles into a funnel that narrows steadily.

The width of the funnel is the first appearance of a law that governs all statistical estimates: the typical deviation of an observed fraction from the true probability is $\sqrt{p(1 - p)/n}$, shrinking in proportion to $1/\sqrt n$. Ten times more data narrows the uncertainty by only a factor of $\sqrt{10} \approx 3.2$. The loaded die in the figure shows why this matters. After 30 rolls, the extra sixes may not yet be distinguishable from chance; after 1000 they are unmistakable.

Probability has a second interpretation, as a degree of belief given the available information. “The probability that this candidate material is stable is 30 %” is not a statement about repeated trials, since the material either is or is not stable; it expresses uncertainty. This Bayesian interpretation obeys exactly the same mathematical rules, and it is the natural language for inference from data, for uncertainty in machine-learning predictions (Gaussian processes) and for deciding which experiment to do next (active learning).

The rules of probability

Probabilities are numbers between 0 (impossible) and 1 (certain). Two rules generate everything else.

For a fair die, the probability of rolling a 1 or a 2 is $1/6 + 1/6 = 1/3$, because the outcomes are exclusive. The probability of two sixes in two independent rolls is $1/6 \times 1/6 = 1/36$. The probability that a roll is even given that it is greater than 3 is $2/3$, because of the three outcomes 4, 5 and 6, two are even.

When all outcomes are equally likely, a probability is a count: the number of favourable outcomes divided by the total number. Counting is therefore central, and its basic tool is the number of ways of choosing $k$ items from $n$ without regard to order,

There are $\binom{10}{3} = 120$ ways to choose which 3 of 10 coin tosses are heads. Statistical mechanics is built on this kind of counting: the entropy of an Einstein solid is the logarithm of the number of ways of distributing energy quanta among oscillators.

Random variables and distributions

A random variable is a number whose value is determined by chance: the number of heads in ten tosses, the speed of a randomly chosen molecule, the error in a measurement. Its behaviour is described by its distribution, which says how probable each value is.

For a random variable that takes separate values, the distribution is a list of probabilities $P(k)$ that add up to 1. The most important example is the binomial distribution: the number of successes $k$ in $n$ independent trials, each succeeding with probability $p$, has

A Galton board makes the binomial distribution physical. Balls fall through rows of pegs; at each peg a ball bounces left or right at random. The bin a ball lands in is the number of rightward bounces, a binomial random variable with $n$ equal to the number of rows. Each individual path is unpredictable, yet the bins fill into a precise shape.

As the number of rows grows, the binomial histogram approaches a smooth bell curve, even when $p \neq 1/2$ and the histogram is lopsided at first. This is a first glimpse of the central limit theorem, discussed below.

Many random variables are continuous: a speed or a position can take any value in a range. The probability of any one exact value is then zero, and the distribution is described by a probability density $p(x)$, defined so that the probability of a value between $x$ and $x + \mathrm dx$ is $p(x)\,\mathrm dx$. Probabilities become areas under the density, $P(a \le X \le b) = \int_a^b p(x)\,\mathrm dx$, and the total area is 1, as described in the chapter on integrals. A density can exceed 1, since it is a probability per unit of $x$; only its integral is bounded.

Expectation and variance

Two numbers summarise most distributions. The expectation or mean is the probability-weighted average of the values, and the variance is the expected squared distance from the mean:

The mean is linear: $\langle aX + bY\rangle = a\langle X\rangle + b\langle Y\rangle$ for any random variables, independent or not. The variance is not linear, but for independent variables the variances add: $\mathrm{Var}(X + Y) = \mathrm{Var}(X) + \mathrm{Var}(Y)$. This one rule explains the $1/\sqrt n$ law. The sum of $n$ independent copies of a variable with variance $\sigma^2$ has variance $n\sigma^2$, so their average, which is the sum divided by $n$, has variance $n\sigma^2/n^2 = \sigma^2/n$ and standard deviation $\sigma/\sqrt n$. Check this against the Galton board: the observed spread of the bins converges to $\sqrt{np(1-p)}$, the square root of the summed variances of $n$ independent bounces.

In statistical mechanics, means are the observable quantities and variances carry physics too. The average energy of a system in contact with a heat bath is its internal energy, and the variance of the energy is proportional to its heat capacity, $\langle E^2\rangle - \langle E\rangle^2 = k_BT^2 C_V$. Fluctuations measured in a simulation are therefore a route to thermodynamic response functions, explored in the chapter on fluctuations.

The normal distribution and the central limit theorem

The bell curve, or normal distribution, has the density

The normal distribution owes its ubiquity to the central limit theorem: the sum (or average) of many independent random variables with finite variance is approximately normally distributed, whatever the distribution of the individual variables. The figure below lets you test this. Choose a distribution for single measurements that looks nothing like a bell: flat, sharply skewed, two separate humps, or just 0 and 1. Then average $n$ of them, many times over, and histogram the averages.

For $n = 1$ the histogram simply reproduces the original shape. By $n = 5$ the averages of uniform numbers already look Gaussian; the skewed exponential needs perhaps 30; the coin needs enough terms for its discrete steps to blur. In every case the width shrinks as $\sigma/\sqrt n$.

This theorem is why measurement errors, which are the sum of many small independent disturbances, are usually Gaussian; why the velocity of a heavy particle buffeted by many collisions is Gaussian; and why the bins of the Galton board fill into a bell. It is also the basis of every error bar. When you report the average of $n$ independent measurements as $\bar x \pm \sigma/\sqrt n$, the standard error of the mean, you are using the central limit theorem to say that the true mean lies within that range with about 68 % confidence, and within twice that range with about 95 %.

Other distributions in science

A few other distributions appear so often that they are worth recognising.

The Poisson distribution, $P(k) = \lambda^k e^{-\lambda}/k!$, describes the number of independent random events in a fixed interval when they occur at an average rate $\lambda$ per interval: radioactive decays per minute, photons arriving at a detector, defects in a length of wire. Its mean and variance are both $\lambda$, so the relative uncertainty of a count is $1/\sqrt\lambda$. A detector that registers 10 000 counts has a statistical uncertainty of 1 %.

The exponential distribution, $p(t) = \lambda e^{-\lambda t}$ for $t \ge 0$, describes the waiting time until the next such event. It is memoryless: having waited a long time does not make the event more imminent. Kinetic Monte Carlo simulations draw the waiting time between rare events from this distribution.

The uniform distribution, constant on an interval, is what random-number generators produce. Every other distribution used in a simulation is made from uniform random numbers by a transformation, such as $-\ln(1 - u)/\lambda$ for an exponential, or by an accept–reject scheme. Monte Carlo methods are built on this.

The Boltzmann distribution $p(s) \propto e^{-E_s/k_BT}$ is the central distribution of statistical mechanics. Its consequence for velocities is the Maxwell–Boltzmann distribution, in which each velocity component is normally distributed with variance $k_BT/m$.

Bayes’ theorem

The product rule can be written in two ways, $P(A \text{ and } B) = P(A \mid B)P(B) = P(B \mid A)P(A)$. Equating them gives Bayes’ theorem, which turns a conditional probability around:

The theorem is simple, but its consequences routinely surprise people. Suppose a fast computational screen is used to search a large set of candidate materials, of which only 5 % are truly good. The screen flags 90 % of the good ones and, wrongly, 10 % of the bad ones. If the screen flags a candidate, what is the probability that it is good? Intuition suggests something close to 90 %. The figure shows a thousand candidates.

Of 1000 candidates, 50 are good and the screen flags 45 of them. Of the 950 bad ones it flags 95. So 140 are flagged, and only 45 of them, 32 %, are good. When the thing being looked for is rare, even a fairly accurate test produces mostly false positives, because the few true positives are outnumbered by the small fraction of a huge number of negatives. The lesson applies to high-throughput screening, to medical tests and to any classifier deployed on imbalanced data: the precision of a prediction depends on the base rate, not only on the accuracy of the test.

Bayes’ theorem is also the foundation of a whole approach to inference. To fit a model to data, one treats the model’s parameters as the hypothesis, specifies a prior, and computes the posterior distribution of the parameters given the data. The posterior says not only which parameters fit best but how uncertain they are. Gaussian processes, Bayesian optimisation and Bayesian neural networks are all applications of this idea.

Estimating averages from data

Every simulation ends with the same task: estimate an average from a finite number of samples and say how uncertain the estimate is. The sample mean $\bar x = \frac1n\sum x_i$ estimates $\mu$, and the sample variance $s^2 = \frac{1}{n - 1}\sum (x_i - \bar x)^2$ estimates $\sigma^2$ (the $n - 1$ corrects for having used the data to locate the mean). If the samples are independent, the standard error of $\bar x$ is $s/\sqrt n$.

Simulation data, however, are almost never independent. Consecutive frames of a molecular dynamics trajectory are separated by a few femtoseconds and are nearly identical; consecutive configurations in a Markov chain Monte Carlo run differ by a single move. Treating them as independent grossly underestimates the error. The correct statement replaces $n$ by an effective number of independent samples, $n_{\text{eff}} = n/(2\tau)$, where $\tau$ is the integrated autocorrelation time in units of the sampling interval: roughly the number of steps it takes for the system to forget its past. A practical way of measuring the true error without computing $\tau$ explicitly is block averaging: divide the trajectory into blocks, average within each block, and treat the block averages as independent samples. As the blocks are made longer than the correlation time, the estimated error rises and then levels off at the correct value. The chapter on trajectory analysis develops this in detail.

Finally, the link between probability and fitting. If measurements scatter around a model’s predictions with independent Gaussian errors of equal width, the probability of the data given the parameters, the likelihood, is a product of Gaussians. Maximising it is the same as minimising the sum of squared residuals. Least squares is therefore maximum-likelihood estimation under Gaussian noise, and the loss functions of machine learning can be read as negative log-likelihoods, an idea developed in learning from data.