hessian-vector products

In machine learning, Hessian and other curvature matrices suffer from high computational complexity, lack of unified representation, and poor scalability, severely hindering applications such as second-order optimization, uncertainty quantification, and model compression. To address this, we propose a unified abstraction paradigm based on linear operators—specifically, matrix-vector multiplication interfaces—and systematically develop the first scalable, user-friendly curvature matrix representation framework. Building upon this paradigm, we introduce *curvlinops*, an open-source PyTorch library that supports automatic differentiation, incorporates structural priors (e.g., sparsity, Kronecker factorization), and maintains compatibility with major deep learning frameworks. Experiments demonstrate that our approach drastically simplifies implementation across curvature-driven tasks, enabling efficient, memory-bounded, and scalable computation even on models with tens of billions of parameters. This advances the practical adoption and democratization of curvature-based methods in large-scale deep learning.

Accurate estimation of the Hessian spectrum for large-scale foundation models has long been hindered by prohibitive computational costs, with existing approaches relying on small models or strong structural approximations. This work proposes a sharded local finite difference method compatible with Fully Sharded Data Parallel (FSDP), combined with stochastic Lanczos quadrature, enabling the first efficient and accurate computation of the true Hessian spectrum on hundred-billion-parameter language models. The method incurs only constant-level training overhead and reveals that widely used block-diagonal approximations severely break down even in moderately sized LLMs. Furthermore, it establishes a scalable paradigm for Hessian spectral analysis, supporting numerical stability studies under fp32/bf16 precision and end-to-end scaling law modeling.

Second-order optimization is hindered by the prohibitive computational cost of Hessian evaluation. This paper proposes an asynchronous, client-separated second-order optimization framework that, for the first time, decouples gradient and curvature (Hessian) computations across distinct clients and updates them asynchronously—eliminating manual hyperparameter tuning required by conventional lazy-Hessian approaches. Integrating cubic regularization with an inexact Hessian update mechanism, the framework effectively alleviates computational latency and communication bottlenecks in high-dimensional settings. We establish rigorous theoretical guarantees of strong convergence. Empirical evaluations on both synthetic and real-world datasets demonstrate consistent superiority over standard and lazy-Hessian baselines, achieving up to √τ real-time speedup while preserving the rapid convergence rate of second-order methods and enhancing practical training efficiency.

For large-scale optimization problems lacking second-order derivative information—such as large-scale data fitting and curve modeling—this paper proposes a parameterized compact Hessian approximation method that replaces the dense Hessian with a low-rank structure, significantly improving scalability. The method introduces the first unified parameterized framework that encompasses and generalizes classical compact quasi-Newton formulas, including BFGS and SR1, and incorporates a flexible vector selection mechanism to better adapt to heterogeneous problem structures. It integrates limited-memory techniques, low-rank matrix decomposition, and quasi-Newton update theory, with task-specific designs for eigenvalue computation, tensor decomposition, and nonlinear regression. Experimental results demonstrate that the proposed method maintains convergence accuracy while substantially reducing memory consumption and computational complexity compared to standard dense or conventional compact Hessian approximations.

Multi-shape Laplacian eigenbases are inherently incompatible due to spectral misalignment across non-isometric shapes, hindering consistent cross-shape functional representation and transfer. Method: This paper proposes a correspondence-guided coupled quasi-harmonic basis construction method. By approximately jointly diagonalizing Laplacian operators and incorporating geometrically consistent correspondence priors—such as stable region indicator functions—it constructs shape-shared, frequency-aligned quasi-harmonic bases. Contribution/Results: The method overcomes the limitation of conventional per-shape independent spectral bases and establishes, for the first time, a geometrically consistent multi-shape frequency-domain representation. It significantly improves accuracy and consistency in cross-shape function transfer across applications including shape editing, pose transfer, non-rigid matching, and shape similarity assessment. The framework provides a scalable theoretical foundation and practical toolset for multi-shape spectral analysis and collaborative processing.

Computing pointwise mean curvature on high-dimensional data manifolds typically incurs a prohibitive $O(m^4)$ complexity due to the explicit construction of the shape operator, hindering scalability. This work proposes an efficient estimation method that circumvents the explicit formation of the Hessian matrix by leveraging an exact algebraic identity. The approach combines covariance eigendecomposition, the cyclic property of the trace operator, and truncated singular value decomposition, while analytically approximating the expectation of nullspace eigenvectors under the Haar measure. This strategy substantially reduces computational complexity, achieving 50–300× speedup on real-world datasets with negligible loss in accuracy, thereby rendering mean curvature a practical geometric feature for machine learning applications.

This work addresses the instability and lack of convergence guarantees in training neural networks on Riemannian manifolds, which arise from neglecting manifold curvature. To this end, the authors propose a curvature-aware training framework that replaces Euclidean derivative supervision with a covariant Sobolev loss, aligns gradients via parallel transport, and incorporates a Laplace–Beltrami regularizer to enhance stability. For the first time, the framework explicitly integrates manifold curvature and the Jacobian of parallel transport into Sobolev-based training, enabling the derivation of geometry-dependent smoothness constants. Building on this, the authors establish curvature-aware linear and quadratic convergence theories. Empirical validation on Lie groups such as SO(3) and SE(3), as well as tasks involving surface imaging and physics-informed learning, demonstrates stable training dynamics and confirms the predicted theoretical convergence rates.

This work investigates the generalization ability of nonlinear least-squares models with ridge regularization at local minima attained after training. Leveraging average algorithmic stability analysis, it characterizes the local geometry of parameter space through the empirical Jacobian Gram matrix and a residual-curvature term. For the first time, an effective dimension is defined with respect to the trained model rather than the initialization point, and combined with curvature and covering complexity to derive a generalization error bound that depends on the geometric structure of learned features rather than the number of parameters. The framework explicitly links this bound to intrinsic data manifold dimensions or the number of stable activation regions in ReLU networks. Experiments confirm Jacobian contraction during training, the tightness of residual-curvature linearization, and a strong alignment between the proposed bound and observed generalization gaps.

This work addresses Riemannian optimization problems with Euclidean bound constraints—common in covariance estimation, neuroimaging, and signal processing—and proposes the first method extending the generalized Cauchy point strategy from L-BFGS-B to Riemannian manifolds, enabling joint optimization over bounded Euclidean variables and manifold variables. The approach integrates limited-memory quasi-Newton updates in the tangent space with a boundary-handling mechanism and is implemented within the Manopt.jl framework. Experimental results demonstrate that the method matches the performance of classical L-BFGS-B, substantially outperforms interior-point methods, and achieves speedups of several orders of magnitude in blind source separation and principal component analysis tasks.

To address the computational inefficiency arising from nested optimization in Fréchet mean estimation on Riemannian and Finsler manifolds, this paper proposes a joint iterative optimization framework based on local coordinate charts: it simultaneously updates the mean location and geodesic distances, thereby eliminating the inner-loop optimization. This work introduces the first gradient-based joint update scheme for both geodesic parameters and the Fréchet mean; extends the formulation to Finsler manifolds; and incorporates an adaptive sampling strategy to enable scalability to large-scale datasets. We establish theoretical guarantees of global convergence and local quadratic convergence rate. Experiments demonstrate that the proposed method significantly outperforms existing baselines in both accuracy and runtime, while comprehensive evaluations across diverse manifold geometries and data scales validate its convergence behavior, robustness, and scalability.