---
title: 'O1NumHess: Constant-Gradient Hessian Reconstruction'
url: https://www.emergentmind.com/topics/o1numhess
type: topic
---

# O1NumHess: Constant-Gradient Hessian Reconstruction

O1NumHess is a seminumerical algorithm for reconstructing the full nuclear Hessian of a molecular system from a number of gradient evaluations that does not grow with system size, asymptotically \(O(1)\) in the number of gradients. It was introduced as a fast and accurate alternative to conventional seminumerical Hessians, which require gradients at \(O(N_{\mathrm{atom}})\) displaced geometries, and it relies on the off-diagonal low-rank (ODLR) structure of molecular Hessians: blocks associated with two distant groups of atoms have low rank rather than sparsity [2508.07544]. In this formulation, O1NumHess is a specific Hessian-reconstruction method for molecular quantum chemistry, not merely a generic label for first-order numerical Hessians.

## 1. Problem setting and motivation

For nuclear coordinates \(\boldsymbol{\xi}=(\xi_1,\dots,\xi_n)\), with \(n=3N_{\text{atom}}\), the target object is the nuclear Hessian
\[
H_{ij} \equiv \frac{\partial^2 E}{\partial \xi_i \partial \xi_j}.
\]
This matrix is required for harmonic vibrational frequencies and normal modes, IR/Raman spectra, thermochemical quantities, reaction-path following and transition-state theory, and vibronic spectra and nonadiabatic dynamics via multiple-state Hessians [2508.07544].

The standard computational routes are analytic Hessians, numerical Hessians from energy finite differences, and seminumerical Hessians from finite differences of gradients. In the standard seminumerical scheme, one displaces one Cartesian coordinate at a time and evaluates
\[
H_{ij} \approx \frac{g_i(\boldsymbol{\xi}+\eta \mathbf{e}_j) - g_i(\boldsymbol{\xi}-\eta \mathbf{e}_j)}{2\eta},
\]
which requires \(6N_{\text{atom}}\) gradient evaluations in the double-sided variant or \(3N_{\text{atom}}+1\) in the single-sided variant [2508.07544]. O1NumHess was proposed to remove this linear growth in gradient calls while still recovering the full Hessian.

The method is motivated by two practical bottlenecks. First, gradients can be expensive for TDDFT, post-HF methods, and large systems. Second, analytic Hessians, although typically fastest and most accurate when available, are difficult to implement, memory-hungry because of CP-SCF with \(3N_{\text{atom}}\) right-hand sides, and often unavailable for excited states and higher-level correlated methods [2508.07544]. O1NumHess targets precisely this regime: full Hessian recovery when gradients are available but conventional seminumerical or analytic Hessians are unattractive.

## 2. Structural basis: locality, low rank, and information compression

The mathematical starting point is the linearized relation between displaced gradients and the Hessian. For displaced geometries
\[
\boldsymbol{\xi}^{(j)} = \boldsymbol{\xi}^{(0)} + \Delta\boldsymbol{\xi}^{(j)},
\]
the gradients satisfy, to linear order,
\[
g_i^{(j)} \approx \sum_{k=1}^{n} H_{ik}\,\Delta\xi_k^{(j)}.
\]
After normalization by the Euclidean norm of each displacement, this becomes
\[
\mathbf{g}=\mathbf{H}\mathbf{\Delta\xi}.
\]
If \(\mathbf{\Delta\xi}\) has full rank \(n\), then \(\mathbf{H}=\mathbf{g}\mathbf{\Delta\xi}^{-1}\), which is the conventional picture and requires at least \(N_{\text{displ}}=n=3N_{\text{atom}}\) displacements [2508.07544].

O1NumHess departs from this dense-matrix viewpoint by using a structural decomposition. The paper argues that the molecular Hessian can be approximated as
\[
\mathbf{H} \approx \mathbf{H}^{\text{local}}+\mathbf{H}^{\text{lowrank}}.
\]
The local part arises mainly from second derivatives of nuclear–nuclear interactions and local electronic contributions; the long-range part is dominated by a small set of near-degenerate electronic states and is therefore low rank. In block form, if atoms are partitioned into distant groups, the corresponding off-diagonal Hessian blocks are not assumed to be small, but they are assumed to have rapidly decaying singular values [2508.07544].

This shifts the effective number of degrees of freedom from \(O(N_{\text{atom}}^2)\) for a generic dense symmetric matrix to \(O(N_{\text{atom}})\): a local block structure contributes \(O(N_{\text{atom}})\) parameters, and a rank-\(r\) symmetric correction \(\mathbf{U}\mathbf{S}\mathbf{U}^T\) contributes \(O(nr)=O(N_{\text{atom}})\) when \(r=O(1)\) [2508.07544]. Since each gradient is an \(n\)-vector, a constant number of gradients then provides \(O(N_{\text{atom}})\) scalar information, which is sufficient, in principle, to determine the independent parameters of the Hessian.

A common misconception is that O1NumHess assumes Hessian sparsity. The data support a narrower statement: the method assumes a local plus low-rank decomposition and, more specifically, ODLR structure. This distinction matters in systems such as conjugated polyenes or excited states, where the Hessian is not sparse but off-diagonal blocks remain numerically low rank [2508.07544].

## 3. Reconstruction algorithm

The algorithm has two principal reconstruction stages: a penalized local fit and an iterative low-rank correction.

The first stage estimates a local Hessian \(\bar{\mathbf{H}}^{\text{local}}\) by minimizing
\[
\text{cost}(\mathbf{H}) = \|\mathbf{g} - \mathbf{H}\mathbf{\Delta\xi}\|_F^2
+ \lambda \| \mathbf{W} \cdot \mathbf{H}\|_F^2,
\]
subject to \(\mathbf{H}=\mathbf{H}^T\) [2508.07544]. The Hadamard-weight matrix \(\mathbf{W}\) penalizes distant couplings according to
\[
W_{ij} = \max\bigl(0,\, r_{ij} - r_i^{\text{vdW}} - r_j^{\text{vdW}} - \Delta r^{(1)}\bigr)^\beta,
\]
with defaults \(\Delta r^{(1)}=1.0\) Bohr, \(\beta=3/2\), and \(\lambda=0.01\) a.u. [2508.07544]. Very distant pairs are set exactly to zero when
\[
r_{ij} > r_i^{\text{vdW}} + r_j^{\text{vdW}} + \Delta r^{(2)}, \qquad
\Delta r^{(2)} = \Delta r^{(1)} + 5\ \text{Bohr}.
\]
This produces a sparse, symmetric positive definite linear system in the unknown non-far Hessian entries. Its dimension is \(N_{\text{nz}}\times N_{\text{nz}}\), with \(N_{\text{nz}}=O(N_{\text{atom}})\), and it is solved iteratively using only matrix–vector products [2508.07544].

The second stage restores the missing long-range information by iterative low-rank updates. Gradients and displacements are first rescaled,
\[
g_{ij}^{\text{scaled}} = g_{ij}\, \frac{\epsilon}{\max\bigl(\epsilon, \sqrt{\sum_i g_{ij}^2}\bigr)},\qquad
(\Delta\xi^{\text{scaled}})_{ij} = (\Delta\xi_{ij})\, \frac{\epsilon}{\max\bigl(\epsilon, \sqrt{\sum_i g_{ij}^2}\bigr)},
\]
with \(\epsilon=10^{-3}\), რათა large forces are deemphasized and low-frequency information is emphasized [2508.07544]. Starting from \(\mathbf{H}_0=\bar{\mathbf{H}}^{\text{local}}\), the update is
\[
\mathbf{H}_{n+1}^{\text{unsym}}=
\mathbf{H}_n + (\mathbf{g}^{\text{scaled}}-\mathbf{H}_n\mathbf{\Delta\xi}^{\text{scaled}})
(\mathbf{\Delta\xi}^{\text{scaled}})^T,
\]
followed by symmetrization,
\[
\mathbf{H}_{n+1}=\frac{1}{2}\bigl(\mathbf{H}_{n+1}^{\text{unsym}}+(\mathbf{H}_{n+1}^{\text{unsym}})^T\bigr).
\]
Each correction has rank at most \(N_{\text{displ}}\), so the iteration explicitly enforces a local-plus-low-rank form [2508.07544].

The paper states that convergence is checked on the residual \(\|\mathbf{g}^{\text{scaled}}-\mathbf{H}_n\mathbf{\Delta\xi}^{\text{scaled}}\|\), with stopping when the norm is below \(10^{-8}\) or its relative change per iteration is below \(10^{-8}\) [2508.07544]. Because the true Hessian is only approximately ODLR and finite-difference or SCF noise is present, the residual is interpreted as a measure of numerical error rather than expected to vanish identically.

## 4. Displacement design and the meaning of “\(O(1)\) gradients”

The constant-gradient claim depends on a specific deterministic construction of collective displacement directions. O1NumHess begins from a model Hessian, \(\mathbf{H}^{\text{Swart}}\), built cheaply in internal coordinates using a modified Swart model and then transformed to Cartesian form [2508.07544]. This model is used only to construct informative displacements.

For non-linear molecules, the initial displacement set contains seven directions: three translations, three rotations about principal axes of inertia, and one symmetric breathing mode [2508.07544]. The translations require no gradients because the forces are zero in the absence of an external field. Rotational gradients at a non-equilibrium geometry are obtained analytically from the gradient at the reference geometry, so they do not require additional electronic-structure calls. The breathing mode is treated with double-sided finite differences because bond-stretch anharmonicity is strong; the remaining directions use single-sided differentiation [2508.07544].

Additional directions are generated through local completeness. For each atom \(A\), the method defines a neighborhood \(\mathcal{N}(A)\) using the near-range cutoff. It projects out the subspace already spanned by current displacements restricted to that neighborhood, diagonalizes the residual local model Hessian, and takes the leading eigenvector as a raw local mode. Local modes from different neighborhoods are then combined into a global mode by choosing signs \(\sigma_A=\pm 1\) greedily to maximize the norm of the global displacement; the resulting vector is orthonormalized and scaled so that its maximum component equals the finite-difference step \(\eta=0.005\) Bohr [2508.07544].

The critical asymptotic observation is that each neighborhood has bounded size for fixed cutoffs, so each local neighborhood requires only a bounded number of independent displacement directions. Once all neighborhoods are locally complete, adding more atoms does not require new global directions. The paper therefore states that, for sufficiently large systems,
\[
N_{\text{displ}} \to N_{\text{displ}}^\infty=\text{constant},
\]
and the total number of gradient evaluations, including equilibrium, breathing-mode doubling, and any post-processing directions, converges to a constant of order 100–120 [2508.07544].

This clarifies the phrase “\(O(1)\) gradients.” It does not mean a single gradient evaluation. It means that the number of gradients saturates at a system-size-independent constant for sufficiently large systems. The paper also states that O1NumHess “requires only about 100 gradients for sufficiently large systems” [2508.07544].

A further robustness step addresses spurious imaginary modes. After a first Hessian is built, the method diagonalizes it, selects a small number of negative-eigenvalue eigenvectors, orthogonalizes them against the displacement set, adds them as extra displacements, and recomputes the Hessian. The number of such additional directions is described as small, so the overall gradient count remains \(O(1)\) [2508.07544].

## 5. Accuracy, scaling, and implementations

The reported numerical behavior is benchmarked on long alkanes and polyenes, WCCR10 transition-metal reactions, and S30L-CI non-covalent complexes [2508.07544]. The paper states that frequency, zero-point energy, enthalpy, and Gibbs free energy errors are only about two times those of conventional double-sided seminumerical Hessians. For large systems, O1NumHess is described as always faster than the conventional numerical Hessian algorithm, frequently even faster than the analytic Hessian, and it requires only about 100 gradients for sufficiently large systems [2508.07544].

For model systems, the reported vibrational-frequency mean absolute deviations are a few \(\text{cm}^{-1}\) for well-behaved ground-state chains, with \(n\)-C\(_{32}\)H\(_{66}\) giving \(0.8\)–\(2\ \text{cm}^{-1}\) depending on \(\Delta r^{(1)}\), and tens of \(\text{cm}^{-1}\) for challenging conjugated excited states such as C\(_{32}\)H\(_{34}\) S\(_1\) and T\(_1\) [2508.07544]. For WCCR10, the MADs of reaction enthalpies and Gibbs free energies are reported as about twice those of conventional double-sided seminumerical Hessians, with reaction \(\Delta G\) MAD \(\approx 1.18\ \text{kcal/mol}\). For S30L-CI, the complexation \(\Delta G\) MAD is reported as \(\approx 0.92\ \text{kcal/mol}\) [2508.07544].

In computational complexity terms, conventional double-sided seminumerical Hessians require \(6N_{\text{atom}}\) gradients, whereas O1NumHess saturates at roughly \(100\)–\(120\) gradients for large systems [2508.07544]. The reconstruction overhead is described as small relative to the gradient computations. The paper also reports that memory use is typically \(1\)–\(10\times\) lower than for analytic Hessians on the tested systems, with smaller savings of \(2\)–\(3\times\) in PCM cases [2508.07544].

The method is implemented in an open-source Python 3 core library, `O1NumHess`, together with a quantum-chemistry wrapper, `O1NumHess_QC`. The wrapper provides nuclear distance matrices, model Hessians, and interfaces to electronic-structure programs; the method is fully interfaced with BDF through the keyword `O1NumHess` in `$bdfopt`, and also interfaced with ORCA through Python [2508.07544]. The core is described as method-agnostic: it assumes only a gradient function for arbitrary coordinate vectors and a distance matrix between coordinates, and it does not assume that coordinates are grouped in threes or that the target function is an electronic energy [2508.07544].

## 6. Interpretation, limitations, and relation to other Hessian methodologies

The defining assumptions are spatial locality, low-rank off-diagonal blocks, and the availability of reliable gradients [2508.07544]. The paper explicitly notes several problematic regimes: very small systems, where constant overhead can dominate; highly delocalized systems such as metallic clusters or extended \(\pi\)-systems with dense low-lying excited-state manifolds, where ODLR ranks may be higher; pathological regions of the potential-energy surface such as conical intersections or very flat surfaces; and noisy gradients from poor SCF convergence or inadequate grids, whose errors may be amplified by the reconstruction [2508.07544].

The article’s terminology can be confused with another usage of “O1 numerical Hessian.” In matrix-based derivative-free optimization, an “order-1 numerical Hessian” denotes a Hessian approximation with \(O(\Delta_u)\) error in the sampling radius. The generalized simplex Hessian in “A matrix algebra approach to approximate Hessians” is such an order-1 numerical Hessian, constructed from function values alone, while its centered variant is order-2 [2304.03222]. This is conceptually distinct from O1NumHess, which is a seminumerical molecular Hessian algorithm using gradient evaluations and ODLR reconstruction rather than a generic order-of-accuracy label.

This distinction suggests a broader classification of numerical Hessian methods. One family, exemplified by generalized simplex Hessians, uses only function evaluations and pseudoinverse-based matrix formulas, yielding order-\(1\) or order-\(2\) approximation guarantees for full or partial Hessians under rank and smoothness assumptions [2304.03222]. O1NumHess belongs to a different family: it uses gradients at collective displacements, exploits problem-specific Hessian structure, and aims to reconstruct the full molecular Hessian with a constant number of gradient calls asymptotically [2508.07544]. A plausible implication is that the name “O1NumHess” should be read historically as a project-specific algorithm name rather than as a statement about first-order truncation accuracy.

Within computational chemistry and beyond, the method is presented as generic for twice-differentiable multivariate functions whenever a distance between coordinates can be defined, far-coordinate Hessian blocks are numerically low rank, and gradients are available at arbitrary points [2508.07544]. The paper specifically mentions classical atomistic potentials, PDE discretizations with integral operators of Calderón–Zygmund type, and machine-learning models with localized basis functions plus global pooling as plausible settings where analogous structure may arise [2508.07544]. This suggests that O1NumHess is best understood not only as a molecular Hessian algorithm, but also as a structured inverse problem: recovering a full dense second-derivative operator from few directional gradient measurements under local-plus-low-rank priors.

Source: https://www.emergentmind.com/topics/o1numhess