Fast Structure-Exploiting Eigendecomposition

• The paper demonstrates that warm-started Jacobi iterations reduce the per-step computational cost from O(d^3) to O(d^2) while preserving numerical accuracy.
• It introduces a dynamic, structure-aware framework that efficiently updates eigenpairs for smoothly evolving Hessians in hierarchical Gaussian process models.
• The method integrates cyclic Jacobi rotations, low-rank updates, and periodic Gram-Schmidt re-orthonormalization to ensure backward stability and scalability.

The infinite-dimensional QR algorithm (IQR) refers to techniques and algorithms that systematically generalize the classical finite-dimensional QR algorithm—central for matrix eigenproblems—to classes of operators or structured matrix sequences in infinite- or large-dimensional settings. Modern incarnations combine algebraic recursion, low-rank updates, dynamic basis-tracking, and iterative diagonalization to address computational bottlenecks in large-scale problems where direct application of dense QR is infeasible, such as Gaussian process models or structured random matrices. The IQR framework is exemplified by dynamically updated eigensolvers for sequences of large symmetric (often indefinite) matrices, as required in Riemannian Manifold Hamiltonian Monte Carlo (RMHMC) or hierarchical Gaussian process inference (Hayakawa et al., 9 Nov 2025).

1. Context and Motivating Applications

The need for scalable eigendecomposition arises in applications where a sequence of dense but structured symmetric matrices must be diagonalized at each step of an iterative algorithm. In hierarchical Gaussian process (GP) models, the Hessian of the negative log-posterior,

$H(q) = \nabla^2_q[-\ln P(q)],$

evolves smoothly with the model parameters $q$, which typically consist of concatenated GP feature coefficients and hyperparameters. These Hessians, despite being $d \times d$ and dense, have properties that enable efficient iterative diagonalization:

• Symmetry (and general indefiniteness).
• Slow temporal evolution: Over small time steps in an integrator (e.g., RMHMC leapfrog), $H(q(t+\epsilon)) \approx H(q(t))$. The IQR concept applies by using the approximate invariance of the eigenbasis across successive steps to drastically reduce computational cost. In practice, repeated full QR or divide-and-conquer eigendecomposition is prohibitively expensive for $d \gtrsim 500$.

2. Mathematical Foundations: Warm-Started Jacobi for Dynamic Structure

At the core is a dynamic-programming-style scheme combining basis recycling and rapid diagonalization by the cyclic Jacobi method. The workflow, instantiated for RMHMC integrators (Hayakawa et al., 9 Nov 2025), is:

At each time step:

1. Basis projection: Given eigenpairs $\Psi(t), \Lambda(t)$ of $H(q(t))$, form the projected matrix

$$A = \Psi(t)^\top H(q(t+\epsilon)) \Psi(t) \approx \Lambda(t) + O(\epsilon)$$

1. Few-sweep Jacobi: Apply 1–2 sweeps of cyclic Jacobi to $A$, accumulating orthogonal $Q$ so that $Q^\top A Q = \Lambda(t+\epsilon)$.
2. Update eigenvectors:

$\Psi(t+\epsilon) = \Psi(t) Q$

1. Orthonormalization: Every $T_0$ steps, perform Gram-Schmidt orthonormalization to mitigate drift.

Pseudocode (indices suppressed for clarity):

A = Psi_old.T @ H_new @ Psi_old Q = np.eye(d) for sweep in range(S): for i in range(d): for j in range(i+1, d): theta = 0.5 * atan2(2*A[i,j], A[i,i] - A[j,j]) c, s = cos(theta), sin(theta) # Apply Givens rotation in (i, j) plane G = ... # Build Givens rotation matrix A = G.T @ A @ G Q = Q @ G if max_offdiag(A) < zeta: break Psi_new = Psi_old @ Q

Only 1–2 sweeps (S ≈ 2) are necessary when warm starting from the previous eigenbasis.

3. Complexity and Performance Analysis

A direct QR or divide-and-conquer eigendecomposition of a $d \times d$ dense matrix costs $O(d^3)$ per step. For high-dimensional Bayesian models, this severely limits tractable system size. The IQR approach achieves:

• Warm-start Jacobi: Each sweep costs $O(d^2)$ since every Givens rotation is $O(1)$ and all pairs $(i < j)$ are updated, yielding $O(S d^2)$ total per step (with small $S$).
• Transformation of Hessian: The $$A = \Psi_\mathrm{old}^\top H_\mathrm{new} \Psi_\mathrm{old}$$ can be computed in $O(d^3)$ or better if block-structure is exploited or hardware accelerates (e.g., cuBLAS, $O(d^\omega)$ with $\omega < 3$).
• Update application: If accumulating $Q$ as explicit dense matrices, costs remain $O(d^3)$, but sequential Givens application keeps cost at $O(S d^2)$ per step.
• Periodic re-orthonormalization: Full Gram-Schmidt costs $O(d^3)$ but is amortized over $T_0$ steps.
• Net per-step cost:

$O(d^2) + O(d^3 / T_0)$

For $T_0 = 10$, the $O(d^2)$ term dominates, conferring an $O(d)$-fold speedup in the high-$d$ regime.

Empirical results confirm that this scheme achieves a $\sim 10 \times$ reduction in wall-clock time over naive eigendecomposition libraries for $d \sim 500$ and above, with equivalent numerical accuracy (Hayakawa et al., 9 Nov 2025).

4. Numerical Stability and Backward Error

The cyclic Jacobi method is well known for backward stability in symmetric diagonalization; quadratic convergence is observed once off-diagonal entries are sufficiently reduced. Critical algorithmic stabilization points in the IQR include:

• Tight halting threshold: An absolute tolerance $\zeta = 10^{-13}$ terminates sweep iterations, ensuring residual off-diagonality is negligible relative to working precision.
• Accumulative Givens rotations: Rather than recomputing from scratch, $Q$ is updated as a product of Givens rotations, avoiding catastrophic loss of orthogonality.
• Periodic Gram-Schmidt: Every $T_0$ steps, columns of $\Psi$ are orthonormalized to correct for small numerical drift due to accumulation of rounding errors.
• Backward error: The posterior eigenpairs $(\Psi, \Lambda)$ satisfy standard backward error criteria to machine precision.

5. RMHMC and Hierarchical GP Model Integration

The principal application motivating these developments is Riemannian-manifold Hamiltonian Monte Carlo for hierarchical Bayesian GP models. RMHMC integrates trajectories in a space whose local geometry is encoded by the Hessian metric, requiring eigenvalues and -vectors at each leapfrog step. Explicitly:

• Metric: The Hessian $H(q)$ supplies the local "mass matrix"—typically via Betancourt's soft-absolute transform.
• Gradient of the metric: Smooth evolution of $q \mapsto H(q)$ justifies the warm-starting protocol.
• Dynamic programming: The algorithm explicitly targets the RMHMC requirement of dynamic, structured eigendecomposition unavailable in generic autodiff frameworks.
• Real-world performance: RMHMC with IQR delivers rapid mixing and robust evidence estimation in logistic regression and propensity modeling with real data, outperforming libraries based on naive static QR routines (Hayakawa et al., 9 Nov 2025).

6. Scaling, Parallelization, and Broader Impact

From an implementation perspective:

• SIMD acceleration: The fixed-size Givens+Jacobi updates are highly compatible with GPU SIMD parallelism, further increasing throughput.
• Scalability: For $d \gg 10^3$, the $O(d^2)$ dominant cost makes leapfrog-based HMC feasible in models well beyond the reach of conventional dense solvers.
• Modular integration: The design encourages a customizable library abstraction, exposing dynamically programmed eigensolver objects that finely tune differentiation and update modes for arbitrary model structures.

The IQR paradigm thus reifies the general philosophy of "structure-aware," dynamic, and re-entrant linear algebra for complex stochastic simulation workflows, laying a practical foundation for large-scale inference in infinite- or high-dimensional Gaussian process problems and related structured probabilistic models (Hayakawa et al., 9 Nov 2025).