Matrix-Free Diagonal Hessian Approximations

Updated 29 June 2026

Matrix-free diagonal Hessian approximations are methods that estimate per-coordinate curvature using only vector operations, avoiding full Hessian formation.
They employ techniques such as curvature propagation, Hessian-vector products, finite differences, and quasi-Newton updates to achieve O(n) computational complexity.
These approaches are critical for large-scale optimization tasks in deep learning, stochastic and derivative-free optimization by providing efficient preconditioning and adaptive step sizing.

Matrix-free diagonal Hessian approximations are computational techniques that provide estimates of the diagonal of the Hessian matrix of a scalar-valued function, using only vector operations or backpropagation passes, without ever forming, storing, or manipulating the full $n \times n$ Hessian. These approximations are extensively used in large-scale optimization, including deep learning, stochastic optimization, derivative-free optimization, and scientific inverse problems, where explicit Hessian formation is computationally infeasible. A range of algorithmic frameworks, including Curvature Propagation, randomized estimators, specialized backpropagation recursions, secant/interpolation models, and diagonal quasi-Newton updates, provide matrix-free access to curvature information for preconditioning, variance reduction, and adaptive step size control.

1. Principles and Motivation

The central motivation is to obtain curvature information sufficient for second-order optimization or uncertainty quantification at a computational and storage cost comparable to a gradient evaluation. The full Hessian $H=\nabla^2 f(x)$ scales quadratically with parameter dimension $n$ , so explicit representation or even Hessian-vector product sweeps quickly become infeasible for large models. In many cases, however, only the diagonal entries $\operatorname{diag}(H)$ are needed. These entries encode per-coordinate curvature and can be exploited as preconditioners in stochastic optimization, adaptive step-size scaling in deep learning, or as a surrogate for trace estimators in statistics and machine learning (e.g., Score Matching).

Matrix-free diagonal Hessian approximations exploit the compositional structure of computational graphs, local block structures, or finite-difference/interpolation schemes to obtain either unbiased or highly accurate biased estimates, always using only vector storage and local operations. Compared to block-diagonal or low-rank approximations, they are strictly $O(n)$ in both complexity and memory, while still providing critical geometric information (Martens et al., 2012).

2. Algorithmic Techniques

A variety of algorithmic paradigms are used for matrix-free diagonal Hessian calculation:

2.1 Curvature Propagation (CP)

Curvature Propagation (Martens et al., 2012) is a general method built on the reverse-mode structure of computational graphs. It introduces random vectors $V$ sampled per node and computes a "curvature vector" through a modified backpropagation. For any scalar-valued function $f:\mathbb R^n \rightarrow \mathbb R$ with computation graph, a single CP sweep yields an unbiased rank-1 approximation $H \approx S(V)S(V)^\top$ , with the $i$ -th diagonal entry $\widehat{H}_{ii}=(S(V)_i)^2$ serving as an unbiased estimator for $H=\nabla^2 f(x)$ 0. Averaging $H=\nabla^2 f(x)$ 1 independent samples reduces variance by $H=\nabla^2 f(x)$ 2.

2.2 Diagonal Backpropagation (BL89, HesScale)

The Becker & LeCun 1989 (BL89) method, extended by HesScale (Elsayed et al., 2024, Elsayed et al., 2022), uses a coupled two-recursion backpropagation to propagate only the diagonal of the second derivative. At each layer, the diagonal Hessian is updated via element-wise operations, neglecting off-diagonal terms, resulting in a deterministic, linear-complexity approximation. HesScale further improves quality by computing the exact diagonal at the output layer, producing empirical errors an order of magnitude lower than stochastic Hutchinson or BL89 methods at the same cost.

2.3 Hessian-Vector-Product (HVP)-Based Extraction

Channel-wise parameters in neural networks (e.g., BatchNorm scale/shift) have exactly diagonal Hessians due to separability. Their diagonal entries can be extracted directly using Hessian-vector products computed via automatic differentiation, requiring only one additional backward pass per parameter block (Sun et al., 2024).

2.4 Secant, Interpolation, and Derivative-Free Models

Secant-based and interpolation approaches use only function and/or gradient values along coordinate-aligned or simplex-defined directions to estimate the Hessian diagonal. For example, the two-sided finite difference approximation

$H=\nabla^2 f(x)$ 3

provides an $H=\nabla^2 f(x)$ 4-accurate, matrix-free estimate for smooth functions (Coope et al., 2020, Jarry-Bolduc, 2021). Specialized positive basis or simplex sets allow for $H=\nabla^2 f(x)$ 5 cost diagonal estimation in derivative-free settings.

2.5 Diagonal Quasi-Newton and Barzilai–Borwein Updates

In large-scale nonconvex and multi-objective settings, diagonal-only quasi-Newton updates using Barzilai–Borwein or secant-type rules (elementwise $H=\nabla^2 f(x)$ 6) provide positive definite diagonal curvature estimates with low overhead and global convergence guarantees (Liu, 20 Dec 2025).

3. Properties, Complexity, and Statistical Guarantees

The matrix-free nature of these diagonal schemes is achieved via three key properties:

Storage: Only $H=\nabla^2 f(x)$ 7 vectors are maintained; no $H=\nabla^2 f(x)$ 8 matrices are stored or factored.
Computation: Each diagonal estimate costs at most the equivalent of one or two gradient or function evaluations. For CP and BL89/HesScale, one or two reverse-mode passes suffice (Martens et al., 2012, Elsayed et al., 2024).
Accuracy/Variance: Curvature Propagation yields unbiased diagonal estimates, and among all unbiased rank-1 schemes, uses noise that minimizes the variance on the diagonal (Martens et al., 2012). HesScale-type deterministic recursions achieve empirically superior bias at the cost of being slightly biased for nondiagonal Hessians, but exploit the fact that layer-block Hessians are nearly diagonal.

The choice of randomization, sample averaging, or exactness can be tailored to application:

Method	Bias/Variance	Cost
Curvature Propagation	Unbiased, $H=\nabla^2 f(x)$ 9 variance decay	$n$ 0–2 $n$ 1grad
BL89/HesScale	Biased, zero variance	$n$ 2–2 $n$ 3grad
HVP extraction	Exact for diagonal blocks	$n$ 4 $n$ 5grad for block
Finite differences	Biased, $n$ 6	$n$ 7– $n$ 8 function evals

Finite-difference and interpolation-based models have $n$ 9 error bounds and inherit the numerical stability and scaling of their underlying sampling patterns (Coope et al., 2020, Jarry-Bolduc, 2021).

4. Empirical Performance and Applications

Empirical results demonstrate that matrix-free diagonal Hessian schemes deliver:

Near-exact diagonal estimates on neural nets with tens to thousands of parameters using only tens to hundreds of samples (CP), or single backward passes (HesScale, BL89).
Superior optimization performance when used to precondition or scale Newton-like or adaptive updates, with wall-clock time speedups and faster convergence versus first-order or block methods (Martens et al., 2012, Elsayed et al., 2024, Elsayed et al., 2022).
In DNN optimization, channel-wise preconditioning via exact diagonal Hessian blocks (SGD-PH) yields systematic improvements on image classification tasks with only a %%%%30 $n$ 31%%%% time and $\operatorname{diag}(H)$ 21.2 $\operatorname{diag}(H)$ 3 memory cost over SGD, outperforming full-block or MC-based approximations in both performance and generalization (Sun et al., 2024).
In derivative-free optimization, diagonal conjugacy-based update mechanisms maintain regularized, well-conditioned scaling in high-noise regimes, outperforming limited-memory full-matrix methods in large benchmark suites (Kimiaei et al., 18 Jun 2026).

Typical performance metrics report a 2–3 $\operatorname{diag}(H)$ 4 reduction in data passes or wall-clock time to a given accuracy over first-order baselines, with essentially no failure cases due to breakdown of curvature estimation.

5. Variants in Structured and Derivative-Free Regimes

Matrix-free diagonal Hessian estimation extends to a variety of specialized settings:

Structured Least Squares: For $\operatorname{diag}(H)$ 5, block-diagonal or pure-diagonal secant approximations use only Jacobian-vector products and can be safeguarded for positive definiteness, yielding globally convergent matrix-free algorithms for $\operatorname{diag}(H)$ 6 up to $\operatorname{diag}(H)$ 7 (Awwal et al., 2020).
Multiobjective/Nonconvex Optimization: Diagonal Barzilai–Borwein quasi-Newton updates furnish low-cost preconditioners shared across all objectives, achieving R-linear convergence and zero failures in large-scale tests (Liu, 20 Dec 2025).
Derivative-Free Optimization: Positive basis or simplex-geometry-based finite-difference/interpolation exploits the O(n) decoupling of coordinate updates, providing accurate Hessian diagonal approximations even when only function queries are available (Coope et al., 2020, Jarry-Bolduc, 2021).

Sophisticated penalty-based diagonal updates using normalized displacement vectors and conjugacy constraints provide robust, noise-insensitive scaling for evolutionary strategies in noisy, high-dimensions (Kimiaei et al., 18 Jun 2026).

6. Limitations, Limiters, and Theoretical Considerations

All matrix-free diagonal Hessian approximations systematically ignore off-diagonal entries. This diagonal restriction is theoretically justified when the parameter coupling is weak, block structures are truly separable (as in channel-wise normalization), or the diagonal dominates (as for many standard network architectures after random initialization or at early optimization stages). In scenarios with strong nonlocal curvature coupling, neglecting off-diagonals can slow convergence or render the preconditioner ineffective. Some methods, such as robustified secant or conjugacy-penalty updates, employ explicit safeguarding, clamping, or regularization to ensure positive definiteness and maintain convergence.

Curvature Propagation is theoretically provable to yield unbiased diagonal (and full) Hessian estimates, with minimal variance among unbiased rank-1 schemes (Martens et al., 2012). HesScale is the leading deterministic approach for DNNs, with deterministic bias but empirically best $\operatorname{diag}(H)$ 8-diagonal error at trivial computational overhead (Elsayed et al., 2022, Elsayed et al., 2024). All methods are fundamentally limited by the quality of the local model, the suitability of their stochastic or coordinate sampling, and the scale of inter-parameter correlation.

7. Impact and Ongoing Developments

Matrix-free diagonal Hessian approximation is pervasive across large-scale optimization and inverse problems, including but not limited to stochastic training of deep nets (Martens et al., 2012, Elsayed et al., 2024, Elsayed et al., 2022), multi-objective optimization (Liu, 20 Dec 2025), nonlinear least squares (Awwal et al., 2020), derivative-free and evolutionary strategies (Coope et al., 2020, Kimiaei et al., 18 Jun 2026), and physical inverse problems with diagonalization via weighting (e.g., ultrasound tomography inversion) (Javaherian, 2022). Their impact is manifest in the scalability of second-order preconditioning, the democratization of curvature-based adaptation to large neural networks, and the feasibility of robust optimization in noise-dominated or model-free scenarios. Ongoing research addresses better interpolation in derivative-free settings, variance-minimized randomized schemes, coupling with block-diagonal or low-rank updates in mixed schemes, and further application to reinforcement learning and large-scale scientific inverse problems (Elsayed et al., 2024, Javaherian, 2022).

The principal advantage across all settings remains the capacity to extract curvature information at $\operatorname{diag}(H)$ 9 complexity, without matrix inversion or storage, delivering strong global convergence guarantees or empirical practical speedups in previously inaccessible problem scales.