🔢 mathematics

Fast and Stable Gradient Approximation for Bilinear Forms of Hermitian Matrix Functions

This paper introduces a fast and stable, forward-only gradient approximation method for bilinear forms of Hermitian matrix functions that reuses the Lanczos pass with minimal overhead, achieving unconditional stability and superior speed compared to existing approaches that rely on costly reorthogonalization or augmented block matrices.

Original authors: Navjot Singh, Kipton Barros, Xiaoye Sherry Li

Published 2026-05-14

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Navjot Singh, Kipton Barros, Xiaoye Sherry Li

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to navigate a massive, foggy mountain range (a giant mathematical matrix) to find the best path to a specific destination. The mountain is so huge that you can't see the whole map at once; you can only take small steps and feel the terrain under your feet. This is how scientists and machine learning models handle huge data problems: they use a clever shortcut called the Lanczos method to estimate the shape of the mountain without mapping every single rock.

However, there's a catch. While this shortcut is great for finding the path forward, figuring out how to change your path if the mountain shifts (calculating a "gradient" or derivative) is a nightmare. Existing methods are like trying to retrace your steps backward through the fog while holding a heavy backpack of extra equipment. They are slow, unstable, and often require you to constantly check your compass (a process called "reorthogonalization") to avoid getting lost, which wastes a lot of time and energy.

The Paper's Big Idea: The "Forward-Only" Shortcut

The authors propose a new way to do this that is like having a GPS that only looks forward. Instead of trying to retrace your steps backward through the complex, foggy mountain, their method says: "Let's just look at the small, simple map we already made while walking forward, and use that to guess how the path should change."

Here is the breakdown using everyday analogies:

1. The Problem: The Heavy Backpack

The Old Way: To figure out how to adjust your route, previous methods tried to run the entire mountain-climbing process in reverse. Because the mountain is so big, this reverse process is unstable. It's like trying to walk backward down a steep, slippery slope while carrying a heavy backpack. To stay safe, you have to stop every few seconds to check your footing (reorthogonalization), which makes the whole trip incredibly slow.
The Alternative: Another method tried to build a giant, double-sized model of the mountain just to calculate the slope. This is like building a massive scale model of the entire mountain range just to figure out which way is up. It works, but it's too heavy and expensive to carry.

2. The Solution: The "Small Map" Trick

The authors realized that the "forward" walk already gave them a tiny, simplified map (a small tridiagonal matrix) that represents the big mountain.

The Analogy: Imagine you walked up a mountain and sketched a tiny, 3-inch drawing of the peak on a napkin. You didn't map the whole mountain, but the napkin captures the essential shape.
The Innovation: Instead of walking backward through the fog, they simply do the math on that tiny napkin. They calculate how the small map changes, and then they "lift" that answer back up to the real mountain.
The Result: They don't need to carry the heavy backpack or build the giant model. They just use the napkin. This is fast and stable.

3. Why It's Safe (The "Ghost" Problem)

In the old methods, when you try to walk backward, you sometimes see "ghosts"—fake paths that look real but aren't. This happens because the math gets messy when you try to reverse the process.

The Paper's Claim: The authors prove that their "napkin method" is safe. The error in their answer is directly tied to how accurate the original forward walk was. If the forward walk was good, the backward guess is good. They show that their method doesn't need the constant "compass checks" (reorthogonalization) to stay stable. It just works, even on modern, fast computer chips.

4. Real-World Tests

The team tested this on three different "mountains":

Gaussian Processes (The Weather Forecast): They tested it on predicting weather patterns (log-determinants). Their method was just as accurate as the slow, heavy methods but much faster.
Network Sensitivity (The Social Web): They tested it on how a change in one friendship affects a whole social network. Their method was 1.4 to 2.6 times faster than the current best method, while being just as accurate.
Quantum Physics (The Atomic Dance): They tested it on learning the rules of how atoms move (Hamiltonian learning). They found that using their method allowed them to learn the rules much faster than using the heavy, dense calculations, without losing accuracy.

Summary

Think of the old methods as trying to solve a puzzle by taking it apart and putting it back together in reverse, which is slow and prone to dropping pieces. The new method is like looking at the puzzle from the front, realizing the pieces fit a certain pattern, and simply adjusting the pattern based on that view.

The Bottom Line: The authors have created a "forward-only" calculator that is fast, stable, and accurate. It allows computers to solve massive, complex math problems involving huge matrices without getting bogged down by the heavy, unstable work of reversing the process. It's a way to get the answer you need without carrying the heavy backpack.

Technical Summary: Fast and Stable Gradient Approximation for Bilinear Forms of Hermitian Matrix Functions

Problem Statement

In scientific computing and probabilistic machine learning, many objectives involve bilinear forms of the type $u^\top f(A(\theta))v$ , where $A(\theta)$ is a large, parameterized Hermitian matrix and $f$ is a scalar function (e.g., matrix exponential, logarithm). These arise in applications such as log-determinant estimation for Gaussian processes, network analysis, and quantum Hamiltonian learning.

While the Lanczos method efficiently approximates the forward evaluation of these forms by projecting the large matrix $A$ onto a small tridiagonal matrix $T$ , computing the gradient with respect to the parameters $\theta$ remains challenging. Existing approaches face significant trade-offs:

Reverse-mode differentiation (Backpropagation): Differentiating through the Lanczos recurrence requires storing all intermediate vectors and often necessitates expensive reorthogonalization to maintain numerical stability. Without reorthogonalization, the gradient calculation can become unstable in finite precision.
Augmented Block Methods: Approaches like applying Arnoldi to a $2n \times 2n$ block matrix (to compute Fréchet derivatives) double the problem size, increasing computational and memory costs.

The core problem is to derive a gradient approximation that is computationally efficient, numerically stable without full reorthogonalization, and avoids the overhead of backpropagating through the iterative Lanczos process.

Methodology

The authors propose a forward-only approximate gradient that reuses the forward Lanczos pass. The method relies on differentiating only the small projected problem rather than the full recurrence.

Core Derivation

Let the quadratic form be approximated by $\hat{\phi}(\theta; u) = \|u\|^2 e_1^\top f(T) e_1$ , where $T = V^\top A V$ is the $m \times m$ tridiagonal matrix from $m$ steps of Lanczos iteration with basis $V$ .

Projected Sensitivity: The sensitivity of the scalar objective with respect to the small matrix $T$ is computed using the Fréchet derivative of the matrix function. This yields a symmetric matrix $G$ such that the differential is $d\hat{\phi} = \text{tr}(G^\top dT)$ .
Ignoring Basis Variation: The total differential of $T$ is $dT = V^\top dA V + dV^\top A V + V^\top A dV$ . The authors approximate $dT \approx V^\top dA V$ , effectively ignoring the terms involving the variation of the Lanczos basis $V$ with respect to $\theta$ .
Error Analysis: The paper proves that the error introduced by ignoring the basis variation terms is proportional to the final Lanczos residual coefficient $\beta_m$ . Specifically, the error term is $2\beta_m e_m^\top G \eta$ . Consequently, the approximation is exact when the Krylov subspace is invariant (i.e., $\beta_m = 0$ ) and the error vanishes as the Lanczos approximation converges.
Gradient Computation: The final gradient is computed as:
$\nabla_\theta \hat{\phi} \approx \sum_{j=1}^m w_j^\top \left[ D_\theta A(\theta, v_j) \right]$
where $W = VG$. This requires only $m$ vector-Jacobian products (VJPs) of the matrix-vector map, reusing the same primitive used in the forward pass.

Algorithm

The proposed Algorithm 1 consists of:

Forward Pass: Run $m$ steps of Lanczos to obtain $V$ and $T$ .
Sensitivity Construction: Compute the eigendecomposition of $T$ , form the divided-difference matrix $F$ , and construct the sensitivity matrix $G$ .
Backward Pass: Compute $W = VG$ and sum the contractions $w_j^\top [D_\theta A(\theta, v_j)]$ .

Key Contributions

Forward-Only Gradient Derivation: The authors derive an approximate gradient that requires no backpropagation through the Lanczos recurrence, only the forward pass and operations on the small projected matrix.
Theoretical Error Bound: They prove that the error of the approximation is governed by the Lanczos residual $\beta_m$ , the same quantity controlling the forward approximation accuracy. This implies the gradient is accurate whenever the forward estimate is accurate.
Numerical Stability: Empirical results show the method inherits the stability of the forward Lanczos quadrature. Unlike adjoint-based methods, it remains reliable without full reorthogonalization, avoiding the instability that plagues naive backpropagation in finite precision.
Efficiency: The method avoids the $O(nk^2)$ cost of reorthogonalization at each step and the $2n \times 2n$ matrix overhead of block-Arnoldi methods, leading to significant speedups.

Experimental Results

The method was evaluated on three benchmarks:

Matrix Logarithm (Gaussian Processes):
- Task: Computing gradients for $u^\top \log(K(\theta))u$ .
- Comparison: Against an adjoint-based Lanczos method [22].
- Result: The proposed method matched the accuracy of the adjoint method (with reorthogonalization) once the forward estimate converged. Crucially, without reorthogonalization, the adjoint method became unstable, while the proposed method remained stable and accurate.
Network Sensitivity (Matrix Exponentials):
- Task: Computing sensitivities of graph communicability and subgraph centrality ( $u^\top \exp(A)v$ ).
- Comparison: Against a modified block-Arnoldi method [24].
- Result: The proposed method achieved speedups ranging from 1.42 $\times$ to 2.64 $\times$ over the block-Arnoldi baseline. At Krylov depth $m=32$ , both methods reached near machine-precision gradient error, but the proposed method was significantly faster.
Quantum Hamiltonian Learning:
- Task: Learning parameters of a transverse-field Ising Hamiltonian.
- Comparison: Against a dense Fréchet-derivative baseline.
- Result: The method achieved a 10.4 $\times$ speedup over the dense baseline while maintaining a parameter error of $1.05 \times 10^{-7}$ at $m=50$ . The optimization trajectory closely tracked the dense solution.

Significance and Claims

The paper claims that this approach offers a "fast and stable" alternative for differentiating bilinear forms of Hermitian matrix functions. Its primary significance lies in decoupling gradient stability from the expensive reorthogonalization procedures typically required for stable Lanczos differentiation. By proving that the gradient error is controlled by the same residual that governs the forward approximation, the authors demonstrate that high-accuracy gradients can be obtained simply by increasing the Krylov subspace dimension, without incurring the computational penalties of full reorthogonalization or augmented matrix methods.

The authors note limitations, specifically that the method is currently restricted to Hermitian matrices and scalar objectives, and that the storage of Lanczos vectors may present a bottleneck for very large-scale problems. Future work is suggested for non-Hermitian extensions and reducing storage requirements.