Hessian Smoothing Techniques

• Hessian smoothing is a collection of techniques that use second derivative information to control curvature and suppress noise across various applications.
• In finite element methods, it employs polynomial-preserving recovery to achieve superconvergence and accurate Hessian estimates, improving mesh adaptation and error estimation.
• In neural networks and optimization, adaptive smoothing and Hessian-regularization enhance convergence rates, generalization, and interpretability even in noisy and high-dimensional settings.

Hessian smoothing encompasses a collection of methodologies for improving the regularity, stability, accuracy, or interpretability of functions, solutions, or estimates by manipulating or exploiting the properties of the Hessian matrix (the matrix of second derivatives). Its applications span numerical PDEs, mesh adaptation in finite element methods, neural network interpretability, geometry processing, inverse problems, optimization, and manifold learning. The treatment of Hessian smoothing varies considerably across disciplines, with each area adapting the notion to its analytical and computational requirements.

1. Theory and Mathematical Frameworks

The core mathematical structures for Hessian smoothing capitalize on the Hessian operator’s sensitivity to function curvature. In PDE discretization and geometry, smoothness energies such as

$$E_{H^2}(u) = \frac{1}{2} \int_\Omega \|H_u\|_F^2 dA$$

penalize roughness, favoring solutions whose second derivatives (curvature) are small in an averaged sense. In manifold learning, generalizations appear as

$H_M(f) = \int_{M} |\mathrm{Hess} f(p)|_F^2 du(p)$

for flat or curved manifolds, where the Hessian is appropriately generalized to Riemannian settings. These functionals lead to minimization problems resulting in splines, denoised data, and regularized estimates.

In the context of optimization, Hessian smoothing can refer to regularization techniques (trace or norm penalization of the Hessian), strategies for adapting smoothing parameters based on Hessian approximation (as in dynamic kernel shaping), or specific averaging schemes to suppress noise in stochastic Hessian estimates.

2. Hessian Smoothing in Numerical Discretization and the Finite Element Method

A principal setting for Hessian smoothing is the post-processing of finite element solutions. In "Hessian Recovery for Finite Element Methods" (Guo et al., 2014), superconvergent recovery of the Hessian is accomplished by the repeated application of a polynomial-preserving recovery (PPR) operator. The process consists of:

1. Applying PPR to the finite element solution $u_h$ to estimate the gradient $\nabla u_h$,
2. Applying PPR again, componentwise, to the recovered gradient to obtain the Hessian $H_h u_h$.

This PPR-PPR method yields several key properties:

• Polynomial preservation: The recovered Hessian exactly reproduces the Hessian for polynomials up to degree $k+1$ on arbitrary meshes, and up to degree $k+2$ or $k+3$ on translation-invariant meshes (depending on the parity of $k$).
• Superconvergence: For sufficiently regular meshes and smooth solutions,

$$\|Hu - H_h u_h\|_{L^2(\Omega)} \leq C h^k \|u\|_{k+2,\infty,\Omega}$$

where $h$ is the mesh size and $k$ the finite element order.

• Symmetry: Under symmetric sampling, the recovered Hessian is symmetric at mesh nodes.
• Empirical validation: Numerical results confirm higher-order accuracy even for low-order elements on structured meshes and competitive or superior performance compared to weighted average or quadratic fit approaches.

Hessian smoothing here is essential for mesh adaptation, error estimation, and for solving fully nonlinear PDEs where the principal part involves the Hessian.

3. Smoothing, Regularization, and Hessian-Based Energies in Geometry and Manifold Learning

In geometric processing and data interpolation, Hessian-based smoothness energies (such as the squared Frobenius norm of the Hessian) offer critical advantages over Laplacian-based energies. Papers such as (Stein et al., 2017) and (Stein et al., 2019) demonstrate that minimizing the Hessian energy subject to data constraints and free (natural) boundary conditions yields solutions that are unbiased by the geometric shape or placement of the boundary. The natural boundary conditions for the Hessian energy ensure that, in the absence of prescribed data on the boundary, solutions default to globally affine (linear) functions, avoiding strong artificial alignment present in Laplacian-based smoothing.

On curved surfaces, as in (Stein et al., 2019), the Hessian is generalized to a covariant one-form Dirichlet energy,

$$E(u)=\frac{1}{2}\int_\Omega (\nabla d u):(\nabla d u) + \kappa|d u|^2\, dA$$

where the curvature term $\kappa|d u|^2$ ensures consistency with intrinsic geometry. Discretizations using mixed finite elements (such as Crouzeix-Raviart) yield practical, convergent algorithms on meshes with complex topology.

In manifold learning, the use of a Hessian-based curvature penalty generalized from thin-plate splines leads to spline estimators represented as combinations of biharmonic Green’s functions and linear null-space components (Kim, 2023). Practical algorithms employ local PCA and Hessian Eigenmaps to directly estimate the penalty from ambient high-dimensional data clouds, sidestepping explicit manifold recovery.

4. Hessian Smoothing in Neural Networks and Black-Box Optimization

For piecewise linear networks, such as those employing ReLU activations, the Hessian is exactly zero almost everywhere. Nevertheless, feature interaction and interpretability analyses benefit from second-order information. The SmoothHess framework (Torop et al., 2023) addresses this by defining feature interactions via the Hessian of a smoothed (Gaussian-convolved) version of the function: $$\nabla^2 h_{f, \Sigma}(x_0) = \int_{\mathbb{R}^d} f(z) \nabla_x^2 q_\Sigma(z - x_0) dz$$ Efficient estimation is realized through Stein's lemma, using stochastic gradient sampling: $$\hat{H}(x_0) = \frac{1}{2n}\sum_{i=1}^n \left[ \Sigma^{-1} \delta_i [\nabla_x f(x_0 + \delta_i)]^T + [\nabla_x f(x_0 + \delta_i)] \delta_i^T \Sigma^{-1} \right]$$ The smoothing parameter (covariance matrix $\Sigma$) explicitly controls the scale of interactions analyzed. SmoothHess empirically outperforms surrogate (SoftPlus, Swish) Hessian approaches and provides superior fit to local Taylor expansions, with rigorously proven sample complexity bounds on the spectral norm error.

Related, Hessian smoothing is leveraged for regularization of deep neural networks. Penalizing the trace of the Hessian (using stochastic estimators such as Hutchinson sampling) biases optimization towards flatter minima, which are empirically associated with improved generalization (Liu et al., 2022).

In derivative-free optimization, smoothing functionals through convolution enable robust gradient and Hessian estimation from noisy function values. Dynamic anisotropic smoothing (Reifenstein et al., 2 May 2024) adapts the covariance structure of the smoothing kernel to align with the (estimated) local Hessian, yielding optimal gradient error scaling and efficient tuning even in heterogeneous curvature landscapes.

5. Hessian Smoothing and Regularization in Optimization Algorithms

In large-scale convex and non-convex optimization, the reliability of Hessian estimates directly influences convergence rates. Innovations such as Hessian averaging (Na et al., 2022) mitigate stochastic noise in Hessian oracles by maintaining (weighted) averages of Hessian estimates across iterations, reducing noise at rate $1/\sqrt{t}$ and enabling local $Q$-superlinear convergence rates: $$\|x_{t+1} - x^\star\| \leq C(\Upsilon \sqrt{\log t / t})^t$$ where $\Upsilon$ represents the normalized stochastic noise level. Weighted averaging schemes optimize transition time to the superlinear regime, essentially matching the best rates achievable with deterministic second-order methods.

Universal frameworks for approximate Newton-type methods, such as those built on Gradient-Normalized Smoothness (Semenov et al., 16 Jun 2025), analytically link Hessian approximation error and function smoothness to global convergence rates. The central iteration

$$x_{k+1} = x_k - (H_k + \frac{\|\nabla f(x_k)\|}{\gamma_k} I)^{-1}\nabla f(x_k)$$

with $\gamma_k$ set according to local properties, achieves optimal rates for H\"older and self-concordant objectives with only inexact (e.g., Fisher or Gauss-Newton) Hessian information, and encompasses both convex and non-convex settings.

6. Hessian Smoothing in the Numerical Solution of Nonlinear and Degenerate PDEs

In the context of fully nonlinear PDEs, such as $k$-Hessian equations, numerical smoothing via Hessian-preserving iterative schemes is critical for stability. The use of subharmonicity-preserving iterations (Awanou, 2014) ensures admissibility, convergence, and selection of viscosity solutions even in the absence of regularity. The enforcement of discrete subharmonicity or $k$-convexity through explicit Hessian smoothing steps regularizes iterates and enhances robustness against singularities or data degeneracy.

In semismooth Newton frameworks for optimization with smoothing operators (Draganescu, 2011), multigrid preconditioners are constructed to efficiently invert principal minors of the Hessian corresponding to inactive constraint sets. These preconditioners exploit the spectral smoothing introduced by operators such as $\mathcal{K}^{*}\mathcal{K} + \beta I$, leading to mesh-size dependent improvements in iteration count, despite suboptimal (but diminishing) spectral distances versus the unconstrained case.

7. Hessian Smoothing in Pluripotential Theory and Nonlinear Complex PDE

In complex analysis, smoothing of $m$-subharmonic functions utilizes solutions of higher-order (complex) Hessian equations. Richberg-type theorems (Pliś, 2013) guarantee the existence of smooth strictly $m$-subharmonic approximations to continuous or strictly $m$-subharmonic functions by solving the complex Hessian Dirichlet problem: $$(dd^c u)^m \wedge \beta^{n-m} = dV,\quad u|_{\partial\Omega} = \varphi$$ for strictly pseudoconvex domains. The smoothing ensures

$u < v < u + h$

for any continuous $m$-subharmonic $u$ and $h > 0$, providing regularization essential for potential theory and the paper of complex Monge-Ampère and Hessian equations.

---

Summary Table: Principal Forms of Hessian Smoothing

Domain	Technique/Formulation	Primary Application
Finite Element Methods	PPR-PPR (polynomial-preserving recovery)	Mesh adaptation, error estimation
Geometry Processing	Minimize $\int \|H_u\|_F^2 dA$ with natural BCs	Smoothing/interpolation on irregular domains
Neural Networks (ReLU)	Gaussian-convolved Hessian via Stein’s Lemma	Feature interaction, interpretability
Optimization (Convex/Nonconvex)	Hessian averaging, GNS-based regularized Newton	Superlinear rate, Hessian inexactness
Derivative-Free Optimization	Adaptive anisotropic smoothing kernel (DAS)	Gradient estimation from noisy values
Pluripotential Theory	Smoothing via solutions to $m$-Hessian equations	Regularization of $m$-subharmonic functions

---

Hessian smoothing unifies a range of analytical and computational approaches for controlling curvature, achieving regularity, and enhancing both the reliability and interpretability of solutions across a broad spectrum of problems in numerical analysis, optimization, machine learning, and applied geometry. Its modern forms employ polynomial-preserving recovery, variational principles, kernel convolution, stochastic estimation, operator-symbol analysis, and deep connections to the qualitative theory of PDEs.