Intro to Optimization - Mathified

by Jonathon Vandezande · Nov 12, 2026

Our previous gamified blog on optimization ignored the math behind it. In this blog, we hope to make optimization theory more concrete while keeping it easy to understand. As such, we have provided a modicum of simple equations and lots of interactive graphics.

We will (without justification) assume that our potential-energy surface is smooth, so that near any point $\mathbf{q}$ it is well approximated by a Taylor-series expansion:

Here,

• $E_q$ is the energy at point $\mathbf{q}$,
• $\delta\mathbf{q}$ is a small step away from $\mathbf{q}$,
• $\mathbf{g}_q$ is the gradient (vector of first derivatives) at $\mathbf{q}$,
• $\mathbf{H}_q$ is the Hessian (matrix of second derivatives) at $\mathbf{q}$,
• and so on with higher-order derivative tensors at $\mathbf{q}$.

We can examine how each of these terms varies across the Müller–Brown PES.

The gradient vector ($\mathbf{g}_q$) shows us the instantaneous slope at the point, and its negative ($-\mathbf{g}_q$) points the way a ball would start rolling on the surface. The Hessian matrix ($\mathbf{H}_q$) describes the curvature. Its diagonal elements give how quickly the slope along each coordinate changes as you move along that coordinate, and its off-diagonal elements give how moving along one coordinate changes the slope along another. Higher-order tensors (generalizations of matrices to higher orders) show higher-order couplings between coordinates.

If we choose a small step $\delta\mathbf{q}$ and contract it with the derivative tensors, then the actual contribution of each term drops off significantly on this PES, as the $n$th-order term is contracted $n$ times with very small values and divided by $n!$.

For small steps, higher-order terms are often inconsequential, and truncating the Taylor-series expansion is a viable option. If we keep only the terms through second order, we get a quadratic model of the surface around $\mathbf{q}$. Setting the model's derivative with respect to $\delta\mathbf{q}$ to zero and solving gives the Newton step:

This step jumps to the stationary point of the quadratic model. That point is a minimum only when the Hessian is positive definite (all of its eigenvalues are positive); otherwise, the Newton step heads toward a saddle point or maximum.

In computational chemistry, gradients are usually only marginally more expensive than energies (roughly one extra SCF cycle for DFT and about twice the cost of an energy for NNPs), but Hessians are often 3N times more expensive than gradients. For minimization, the extra information gained from computing the Hessian is rarely worth the cost. Instead of computing the true Hessian at the start of an optimization, we often start with a guess (e.g., Schlegel, Lindh, or the identity matrix) and iteratively update it as we optimize. The most common update, BFGS (Broyden–Fletcher–Goldfarb–Shanno), uses the change in the gradient between consecutive steps to estimate the curvature along the step direction. It preserves positive definiteness if the initial Hessian is positive definite and each update satisfies the curvature condition.

Quasi-Newton minimization with BFGS works well in the approximately quadratic region of a potential well, where the true Hessian is positive definite. Near a transition state or in a non-quadratic region, two problems arise. First, a BFGS approximation initialized positive definite cannot learn negative curvature. The explorer skips updates that fail the curvature condition. An exact initial Hessian can already contain negative curvature, but skipped updates can still limit what BFGS learns. Second, higher-order terms can become significant.

The first problem is fixed by using Bofill updates, which can represent negative curvature, and rational function optimization (RFO) steps. RFO shifts the Hessian by a scalar, $\delta\mathbf{q} = -(\mathbf{H}_q - \mu\mathbf{I})^{-1} \mathbf{g}_q$, with $\mu$ chosen below both zero and the lowest eigenvalue. The shifted Hessian is then positive definite, so the step goes downhill along every mode. The shift does not guarantee a short step: when the gradient has only a small component along a mode with negative curvature, the RFO step can be very long, so optimizers also cap it at a trust radius. A partitioned variant (P-RFO) instead goes uphill along the lowest mode to find a transition state. The second problem is often mitigated by taking smaller steps or using nonlinear coordinate transforms.

The explorer below draws all three steps from the exact Hessian at a point you choose. Turn on the quadratic model to see its contours: ellipses in a well and hyperbolae where the Hessian is indefinite. In a well, the Newton and RFO steps point in nearly the same direction, while P-RFO climbs along the lowest mode. Where the Hessian is indefinite, as it is near a saddle point, Newton and P-RFO both step towards the model's saddle point, while RFO steps downhill.

Further improvements can be made by taking into account the unique nature of molecular systems. Modern computational-chemistry optimizers generate initial guess Hessians that use empirical estimates for the second derivatives of bond stretches, angle bends, and dihedral twists (Schlegel, 1984). Nonlinear transformations to internal coordinates based on these stretches, bends, and twists can reduce anharmonic coupling between molecular motions, helping the quadratic approximation remain accurate over larger displacements and allowing optimizers to take larger steps (Fogarasi et al., 1992). For nonbonded species, whole-body translational and rotational coordinates are also added (Wang and Song, 2016). Further refinements can be made to how these transformations are applied or how steps are taken. These tricks have been added to the Rowan optimizer but are beyond the scope of this blog.

For fun minigames showcasing optimization theory, including a transition-state search, see our companion blog Intro to Optimization - Gamified.