---
title: Hessian-Based Regularization
url: https://www.emergentmind.com/topics/hessian-based-regularizer
type: topic
---

# Hessian-Based Regularization

A Hessian-based regularizer is a penalty term constructed from second-order derivatives and added to an optimization objective to control curvature. In contemporary usage, the Hessian may be taken with respect to model parameters, network inputs, image intensities, or functions on manifolds, and the penalty may involve the trace, operator norm, Frobenius norm, Schatten norm, or a distributional total-variation analogue of the Hessian. Across these settings, the common aim is to bias optimization toward low-curvature solutions, flatter minima, improved adversarial robustness, or piecewise-linear and structure-preserving reconstructions [2208.05924, 2012.03801, 2009.06571, 1209.3318, 2302.05025, 2208.07787].

## 1. Conceptual and theoretical basis

In deep learning, one of the clearest motivations for Hessian regularization is an explicit generalization-error bound. Liu et al. quote a linear-model result attributed to Wei et al. in which both the average Jacobian norm and the average Hessian trace appear in the bound. Writing $\mu(W)=\hat{\mathbb E}[\|\mathrm{Jacobian}(f(x;W))\|_2]$ and $\nu(W)=\hat{\mathbb E}[\mathrm{tr}(\mathrm{Hessian}(f(x;W)))]$, the bound includes terms of the form
$$
(A\mu)^{2/3}(\theta B)^{1/3} n^{-1/3}
\;+\;
A\sqrt{B\nu\theta/n}
\;+\;
BA^2\theta/[n(\log^2(BA^2\theta/(\nu n))+1)]
\;+\;\zeta,
$$
so that keeping $\mu(W)$ and $\nu(W)$ small helps control the generalization gap [2208.05924].

A second motivation is geometric. Around a local minimizer $\omega^\*$, the second-order Taylor approximation
$$
\ell(\omega)\approx \ell(\omega^\*) + \tfrac12(\omega-\omega^\*)^T H(\omega^\*) (\omega-\omega^\*)
$$
makes the eigenvalues of $H(\omega^\*)$ a direct measure of sharpness. Penalizing $\sum_i \lambda_i=\mathrm{tr}(H)$ biases optimization toward regions with many small eigenvalues, i.e. flat minima. Liu et al. further connect this to the gradient-descent ODE
$$
\frac{d\omega}{dt}=-\nabla \ell(\omega),
$$
whose equilibrium stability is governed by $J=-\nabla^2\ell(\omega^\*)$: penalizing $\mathrm{tr}(H)$ reduces the sum of positive eigenvalues of $H$ and thus weakens Lyapunov stability, which is presented as making it easier for SGD to escape “easily-converged” but suboptimal minima [2208.05924].

A layerwise viewpoint yields a related but more localized interpretation. In “A Deeper Look at the Hessian Eigenspectrum of Deep Neural Networks and its Applications to Regularization,” the Hessian eigenspectrum of each layer is reported to be largely similar to that of the entire network, with the eigenspectrum of middle layers observed to be most similar to the overall Hessian eigenspectrum. The same work reports that both the maximum eigenvalue and the trace of the Hessian reduce as training progresses, and uses these observations to motivate penalizing layerwise Hessian traces, including a middle-layer-only variant [2012.03801].

## 2. Principal formulations

The simplest parameter-space formulation penalizes the trace of the empirical-loss Hessian with respect to parameters:
$$
R(\omega)=\mathrm{tr}[\nabla^2_\omega \ell_{\mathrm{emp}}(f(\cdot;\omega))], \qquad
L(\omega)=\ell_{\mathrm{emp}}(f)+\lambda R(\omega).
$$
This is the stochastic Hessian trace regularization framework denoted SEHT by Liu et al. [2208.05924].

A layerwise variant partitions parameters as $\theta=\{\theta_1,\dots,\theta_L\}$ and defines block Hessians
$$
H^{(l)}(\theta):=\frac{\partial^2\mathcal L_0}{\partial \theta_l\partial \theta_l}\in\mathbb R^{|\theta_l|\times |\theta_l|}.
$$
The regularizer is then
$$
R(\theta)=\sum_{l=1}^L \mathrm{Tr}[H^{(l)}(\theta)],
$$
with total objective
$$
\mathcal L_{\mathrm{total}}(\theta)=\mathcal L_0(\theta)+\gamma\sum_{l=1}^L \mathrm{Tr}[H^{(l)}(\theta)].
$$
Because middle layers were found to track the full-network spectrum most closely, the same paper proposes the restricted objective
$$
\mathcal L_{\mathrm{mid}}(\theta)=\mathcal L_0(\theta)+\gamma\sum_{l\in M}\mathrm{Tr}[H^{(l)}(\theta)],
\quad
M=\{l:\lfloor L/4\rfloor\le l\le \lfloor 3L/4\rfloor\},
$$
thereby penalizing only the middle half of layers [2012.03801].

A noise-smoothed formulation replaces the original loss $f(w)$ by
$$
F(w)=\mathbb E_{u\sim \mathcal N(0,\sigma^2 I)}[f(w+u)].
$$
Using a second-order Taylor expansion, the paper “Noise Stability Optimization for Finding Flat Minima: A Hessian-based Regularization Approach” states
$$
\mathbb E_u\!\Bigl[\tfrac12(f(w+u)+f(w-u))\Bigr]
=
f(w)+\tfrac12\sigma^2\,\mathrm{tr}(\nabla^2 f(w))+O(\|\sigma^2 I\|^{3/2}),
$$
so that, to first order in $\sigma^2$, the smoothed objective behaves like
$$
F(w)\approx f(w)+\frac{\sigma^2}{2}\,\mathrm{tr}(\nabla^2 f(w)).
$$
In this formulation, the Hessian penalty is induced by isotropic Gaussian weight perturbations rather than inserted explicitly as a symbolic trace term [2306.08553].

These formulations already indicate that “Hessian-based regularizer” is not a single construction. In the supplied literature it denotes at least three parameter-space families: explicit trace penalization, layerwise trace penalization, and stochastic smoothing schemes whose leading correction is proportional to the Hessian trace.

## 3. Matrix-free estimation and optimization

The dominant computational obstacle is that modern Hessians are too large to form explicitly. The standard matrix-free identity is Hutchinson’s estimator:
$$
\mathrm{tr}(H)=\mathbb E_v[v^T H v],
$$
for random $v$ with $\mathbb E[v_i]=0$ and $\mathbb E[v_i v_j]=\delta_{ij}$; the supplied works use Rademacher or Gaussian probes. The corresponding Hessian-vector product is obtained without forming $H$ by
$$
Hv=\frac{\partial}{\partial \omega}[\nabla \ell(\omega)\cdot v],
$$
which requires two automatic-differentiation calls. In the layerwise setting, repeating this for $n_{\mathrm{hutchinson}}$ probes yields an unbiased estimator of each layer trace, and the paper reports that $n_{\mathrm{hutchinson}}=1$ or $2$ together with an update frequency $f_r\approx 50$–$100$ suffices in practice, leading to $\approx 2\times$ slowdown compared to plain SGD [2208.05924, 2012.03801].

Liu et al. further accelerate trace estimation with a dropout-style sparse probe. They define a masked random vector $\sigma$ with
$$
\Pr(\sigma_i=\pm 1)=p,\qquad \Pr(\sigma_i=0)=1-2p,
$$
so that, conditional on the mask, $\sigma^T H \sigma$ is an unbiased estimator of the partial trace over selected coordinates. Their SEHT-D variant samples a layer-wise mask with keep probability $p_1$, then a within-layer keep probability $p_2$, forms $\sigma$ on the kept parameters, computes $v=\nabla \ell\cdot \sigma$, obtains $h=\partial_\omega v$, and accumulates $t=\sigma^T h$. Its per-iteration cost is stated as approximately $(1+K\cdot p)\times$ a standard update, where $K$ is the number of Hutchinson probes and $p$ the effective retained fraction [2208.05924].

Noise-based trace regularization uses a different estimator. The NSO algorithm draws perturbations $u_t^{(j)}\sim \mathcal N(0,\sigma^2 I)$ and averages symmetric stochastic gradients
$$
G_t^{(j)}=g(w_t+u_t^{(j)})+g(w_t-u_t^{(j)}),\qquad
\bar G_t=\frac1{2k}\sum_{j=1}^k G_t^{(j)},
$$
followed by $w_{t+1}=w_t-\eta_t \bar G_t$. The two-point symmetrization is designed to cancel the first-order Taylor noise term that affects naive one-point perturbation estimates. When $k=1$, the paper states that each step costs two back-propagations, the same as SAM [2306.08553].

When the target is a spectral norm rather than a trace, power-method and Lanczos iterations replace Hutchinson sampling. For input-Hessian operator-norm regularization, Mustafa et al. optimize
$$
\|H_f(x)\|_{p\to q}=\max_{\|v\|_p=1}\|H_f(x)v\|_q
$$
by projected gradient ascent over $v$, using Hessian-vector products computed through Pearlmutter’s algorithm. In their reported configuration, $T=10$ inner iterations produce a per-batch cost multiplier of about $(1+T)\simeq 11$ [2009.06571].

The paper “Generalizing and Improving Jacobian and Hessian Regularization” replaces power iteration by a batched, matrix-free Lanczos method for the extremal eigenvalue of
$$
M=(H-T)(H-T)^T.
$$
The stated objective is
$$
R_H(\theta,x)=\|H_f(\theta;x)-T\|_2
=
\sigma_{\max}(H_f-T)
=
\sqrt{\lambda_{\max}((H_f-T)(H_f-T)^T)}.
$$
Lanczos only requires efficient matrix-vector products and supports batched HVP/JVP/VJP implementations on GPU; the paper reports that the overhead relative to a single HVP is about $14.4$ s per Lanczos iteration on an A100, with total overhead ranging from $48$–$385\%$ depending on the number of iterations [2212.00311].

## 4. Decomposition of the Hessian and the sharpness debate

A major refinement of Hessian-based regularization is the observation that the parameter-space Hessian is not monolithic. The paper “Neglected Hessian component explains mysteries in Sharpness regularization” gives the exact decomposition
$$
H(\theta)=J(\theta)^T H_z J(\theta)+\nabla_z L(\theta)\cdot \nabla_\theta^2 z(\theta),
$$
and interprets the two terms differently. The Gauss–Newton term $J^T H_z J$ is positive semidefinite and is described as “feature-exploitation,” while the nonlinear modeling error matrix
$$
\mathrm{NME}(\theta)\eqdef \nabla_zL(\theta)\cdot \nabla_\theta^2 z(\theta)
$$
is generally indefinite and is interpreted as “feature-exploration.” The NME vanishes at interpolation points where $\nabla_z L=0$, which explains why many analyses neglect it, but the paper argues that modern networks do not remain in a regime where this term can be ignored [2401.10809].

This decomposition changes how several regularizers are interpreted. A gradient penalty update contains both $GN\cdot g$ and $NME\cdot g$. Because the NME depends on $\phi''(\cdot)$, the effectiveness of explicit Hessian-vector-product regularization becomes activation-dependent. The paper states that if $\phi''\equiv 0$, as with practical ReLU implementations, the NME is almost zero on the diagonal, and gradient penalties can become negligible or high-variance; smooth activations such as GELU yield a non-zero NME during training and make the penalty effective. By contrast, SAM is described as implicitly integrating both GN and NME over a neighborhood, making it substantially less sensitive to pointwise sparsity of $\phi''$ [2401.10809].

The same work directly challenges a common equivalence claim. It states that the long-held equivalence between weight noise and gradient penalties relies on ignoring the NME, and that this assumption does not hold for modern networks since they involve significant feature learning. Its ablations report that penalizing $\mathrm{tr}(GN)$ alone improves generalization, whereas penalizing the full $\mathrm{tr}(H)$, i.e. $GN+NME$, gives little or no benefit. A plausible implication is that “Hessian-based regularizer” must be specified at the level of Hessian components, not only at the level of a scalar summary such as trace or spectral norm [2401.10809].

## 5. Input-space curvature penalties and structured Hessian constraints

Another major branch of the literature regularizes the Hessian with respect to the network input rather than the parameters. For a classifier with scoring function $f:\mathbb R^d\to\mathbb R^k$, Mustafa et al. study margin functions $f^{(j)}(x)=f_t(x)-f_j(x)$ and define the Hessian of the margin as
$$
H_f(x)=\nabla_x^2 f(x).
$$
Their regularized objective adds both input-gradient and input-Hessian penalties:
$$
L(\theta)=\mathbb E_{(x_i,y_i)}\Bigl[
\ell(f(x_i),y_i)
+
\lambda_1 \sum_{j\ne y_i}\|\nabla_x f^{(j)}(x_i)\|_q
+
\lambda_2 \sum_{j\ne y_i}\|H_{f^{(j)}}(x_i)\|_{p\to q}
\Bigr].
$$
For $p=q=2$, the Hessian term is the spectral norm. The associated second-order robustness bound introduces
$$
K_R(f,x)=\max_{\|y\|_p\le R}\|H_f(x+y)\|_{p\to q},
$$
and states that a large adversarial perturbation requires both the local gradient and the local curvature to be small. This places Hessian regularization in direct continuity with gradient-norm defenses rather than as a separate paradigm [2009.06571].

The generalized-target framework of Lim et al. broadens the design space further. Instead of regularizing only toward the zero matrix, they consider
$$
R_H(\theta,x)=\|H_f(\theta;x)-T\|_\*,
$$
with $T$ any matrix admitting efficient left and right vector products. This yields explicit structural regularizers. For diagonality, Theorem 1 in the paper motivates
$$
R_{\mathrm{diag}}(\theta,x)=\|H_f-D(H_f\mathbf 1)\|_2,
$$
which drives off-diagonal Hessian entries toward zero. The same framework is also used for Jacobian symmetry. The paper emphasizes that Lanczos-based spectral-norm minimization is effective and stable for large Jacobian and Hessian matrices, including input Hessians of size $3072\times 3072$ on CIFAR-scale experiments [2212.00311].

Empirically, the input-Hessian and spectral-norm line of work is evaluated mainly through adversarial robustness rather than standard clean-error minimization. In Mustafa et al., Cross-Hölder regularization improves robustness over input-gradient regularization on MNIST and Fashion-MNIST, while in Lim et al. Lanczos-based Hessian spectral-norm regularization on ResNet-18 yields better robust accuracy than Hutchinson- or power-method baselines in the reported CIFAR-10 and CIFAR-100 experiments [2009.06571, 2212.00311].

## 6. Function-space, manifold, and imaging formulations

Outside neural-network parameter space, Hessian regularization appears as a function-space curvature penalty. On a flat Riemannian manifold $M$, the paper “Hessian Based Smoothing Splines for Manifold Learning” defines
$$
R[f]=H_M(f)=\int_M \|\mathrm{Hess}\,f(p)\|_F^2\,d\mu(p),
$$
with the underlying space taken as $W^{2,2}(M)$. The null space
$$
N=\{f\in W^{2,2}(M):\mathrm{Hess}\,f\equiv 0\}
$$
consists of functions that are locally affine. The corresponding smoothing spline solves
$$
f_{SS}=\arg\min_{f\in W^{2,2}(M)}
\sum_{i=1}^N (y_i-f(p_i))^2+\lambda\int_M \|\mathrm{Hess}\,f\|_F^2\,d\mu,
$$
and the minimizer is expressed through biharmonic Green’s functions and basis functions spanning the null space. When the manifold is unknown, the paper uses the Hessian-Eigenmaps estimator to build a discrete sparse matrix $H_M(X)$ so that
$$
g(X)^T H_M(X) g(X)\approx \sum_{i=1}^N \|\mathrm{Hess}\,f(p_i)\|_F^2,
$$
leading to the closed-form quadratic smoother
$$
\hat g=(I_N+\lambda H_M(X))^{-1} y
$$
or, with diagonal weights $W$, $\hat g=(W+\lambda H_M(X))^{-1}Wy$ [2302.05025].

In inverse imaging, Lefkimmiatis et al. propose per-pixel Hessian Schatten-norm regularization. For a discrete image $x\in\mathbb R^N$ with per-pixel Hessian $(Hx)_n\in\mathbb R^{2\times 2}$, the regularizer is
$$
R_p(x)=\sum_{n=1}^N \|(Hx)_n\|_{S_p}.
$$
This is a second-order analogue of total variation: it penalizes curvature rather than gradient magnitude and does not penalize linear ramps, so the paper presents it as a way to avoid the staircase effect associated with TV. The same work emphasizes translation invariance, rotation invariance, scaling invariance up to a multiplicative factor, and convexity for all Schatten $p$-norms with $p\ge 1$ [1209.3318].

Generalized Hessian–Schatten norm regularization combines Hessian–Schatten and TGV-style ideas. Its dual form is
$$
{\cal GHS}_p(g;\alpha_s,\alpha_f)=
\max_{\substack{
\|\overline{\mathbf N}\|_{\infty,S(q)}\le \alpha_s,\\
\|\mathrm{div}\,\overline{\mathbf N}\|_{\infty,2}\le \alpha_f}}
\langle \overline{\mathbf N}, \mathbf H*g\rangle,
$$
and its primal form is
$$
{\cal GHS}_p(g;\alpha_s,\alpha_f)=
\min_{\mathbf u}
\alpha_f\sum_{\mathbf r}\|\mathbf d*g(\mathbf r)-\mathbf u(\mathbf r)\|_2
+
\alpha_s\sum_{\mathbf r}
\Bigl\|
\tfrac12(\mathbf u*d^T+\mathbf d*\mathbf u^T)(\mathbf r)
\Bigr\|_{S(p)}.
$$
The authors note that as $\alpha_f\to \infty$ this reduces to the HS norm, and for $p=2$ to TGV-2 [2105.12099].

A non-convex variant replaces the convex $\ell_p$-type penalty on Hessian singular values by the $q$-shrinkage penalty $g_q$, defined indirectly through its proximal mapping
$$
\mathrm{prox}_{\rho g_q}(x)
=
s_q(x)
=
\max\{|x|-\rho^{2-q}|x|^{1-q},0\}\,\mathrm{sign}(x).
$$
The resulting regularizer is
$$
\Phi(u)=\sum_x [g_q(\sigma_1(x))+g_q(\sigma_2(x))],
$$
where $\sigma_1(x),\sigma_2(x)$ are the singular values of the Hessian at pixel $x$. The associated ADMM algorithm relies on per-pixel singular-value shrinkage, and convergence to a stationary point is derived through restricted proximal regularity, following the framework attributed in the paper to Wang–Yin–Zeng [2309.04593].

For continuous and piecewise-linear regression, Hessian total variation provides another function-space interpretation. Pourya et al. define, for $f\in C^2(\Omega)$,
$$
HTV(f)=\int_\Omega \|\nabla^2 f(x)\|_{S_1}\,dx.
$$
When $f$ is continuous piecewise-linear over a Delaunay triangulation, the generalized Hessian is concentrated on simplex facets, and the paper derives the closed form
$$
HTV(f)=\sum_{(A,B)\in N_T} |u_{A,B}^T(a_A-a_B)|\,Vol_{A,B},
$$
where $a_A$ and $a_B$ are simplex gradients. In the discrete CPWL model, this becomes
$$
HTV(f)=\|Lc\|_1,
$$
so that the learning problem reduces to a generalized LASSO
$$
c^\* \in \arg\min_{c\in\mathbb R^{N_g}}
\tfrac12\|y-Hc\|_2^2+\lambda\|Lc\|_1.
$$
Here the Hessian regularizer does not merely smooth; it enforces sparsity of second-order variations and therefore favors CPWL functions with few affine pieces [2208.07787].

## 7. Empirical behavior, advantages, and limitations

Parameter-space Hessian-trace penalties are reported to improve generalization on both vision and language tasks. On CIFAR-10 with ResNet-18, Liu et al. report $94.00\pm 0.24\%$ accuracy for the baseline with weight decay only, $95.37\pm 0.09\%$ for SEHT-D with $(K=1,p=0.01)$, $95.49\pm 0.06\%$ for SEHT-D with $(K=10,p=0.05)$, and $95.59\pm 0.06\%$ for full-net SEHT-H $(K=5)$. On CIFAR-100 with WRN-28-10, they report a baseline of $74.61\pm 0.52\%$ top-1 and $80.31\pm 0.33\%$ for SEHT-D $(K=1,p=0.05)$. On WikiText-2, SEHT-D reduces test perplexity from $95.65$ to $94.86$ for LSTM and from $111.64$ to $109.03$ for GRU [2208.05924].

Layerwise Hessian-trace regularization shows smaller but consistent gains in standard image classification. With $\gamma=1e{-2}$, the layerwise paper reports test error reductions such as VGG11 on CIFAR-10 from $18.20\pm 0.45$ to $15.11\pm 0.22$, and VGG11-BN on CIFAR-100 from $45.25\pm 0.44$ to $44.98\pm 0.05$. Penalizing only the middle half of layers is reported to yield nearly identical or slightly better test errors than penalizing all layers [2012.03801].

Noise-stability optimization reports improvements in both accuracy and explicit curvature statistics. The abstract states that the method delivers up to a $2.4\%$ test accuracy increase for fine-tuning ResNets on six image classification datasets, reduces the trace of the Hessian by $15.8\%$, and reduces the largest eigenvalue by $9.7\%$. The paper also reports that the method remains effective for improving generalization in pretraining multimodal CLIP models and chain-of-thought fine-tuning, and that it combines effectively with weight decay and data augmentation [2306.08553].

For robustness-oriented input-Hessian regularization, the main limitation is compute. Mustafa et al. explicitly state that direct Hessian computation is infeasible and that their Pearlmutter-based approximation still multiplies cost by about $(1+T)\simeq 11$ in the reported configuration. They nevertheless report that Cross-Hölder regularization increases robustness over input-gradient regularization on MNIST and Fashion-MNIST while maintaining clean accuracy of approximately $99\%$ on MNIST and approximately $88\%$ on Fashion-MNIST [2009.06571]. Lim et al. likewise emphasize the high computational complexity of Jacobian and Hessian regularization and motivate Lanczos specifically as a stable large-matrix alternative to Hutchinson- and power-method approximations [2212.00311].

The supplied literature also records important caveats. The layerwise trace paper notes that applicability to very large models, such as transformers, or other domains remains to be tested [2012.03801]. The sharpness-decomposition paper argues that explicit full-Hessian penalties can regularize the wrong component, and that penalizing feature exploitation without penalizing feature exploration is often the correct choice [2401.10809]. This suggests that empirical success depends not only on the magnitude of curvature being penalized, but on which curvature component is being suppressed.

In imaging and geometric learning, the empirical record is similarly application-specific. Hessian Schatten-$1$ and Hessian Schatten-$2$ regularizers are reported to outperform TV and Haar by approximately $0.5$ dB on average in inverse imaging, while preserving smooth ramps and removing staircase artifacts [1209.3318]. GHSN with $p=1$ is reported to give the highest PSNR among the compared higher-order image-reconstruction methods [2105.12099]. The non-convex shrinkage variant yields average SSIM values of $0.962/0.913$ under the two MRI masks reported in the paper, compared with $0.929/0.836$ for TV-1 and $0.938/0.862$ for TV-2 [2309.04593]. In low-dimensional regression, Delaunay-triangulation-based HTV regularization achieves test MSE $\simeq 13.14$ on the 4-D Power-Plant dataset and is described as using far fewer parameters than the compared neural networks [2208.07787].

Taken together, these results establish Hessian-based regularization as a broad second-order family rather than a unitary method. In the supplied work it serves at least four distinct roles: flattening parameter-space loss landscapes, constraining input-space curvature for robustness, enforcing matrix structure through spectral penalties, and controlling function-space bending or second-order sparsity in manifold learning and inverse problems. The main recurring trade-off is equally broad: stronger curvature control is repeatedly associated with improved generalization, robustness, or structural fidelity, but almost always at nontrivial computational cost and with sensitivity to the specific Hessian quantity being penalized.

Source: https://www.emergentmind.com/topics/hessian-based-regularizer