Stochastic Zeroth-Order Method for Computing Generalized Rayleigh Quotients

Technical Summary: Stochastic Zeroth-Order Method for Computing Generalized Rayleigh Quotients

Problem Statement
The paper addresses the maximization of the generalized Rayleigh quotient, defined as:
$$R(A, B) = \max_{v \in \mathbb{R}^d \setminus \{0\}} \frac{\langle v, Av \rangle}{\langle v, Bv \rangle}$$
where $A \in \mathbb{R}^{d \times d}$ and $B \in \text{Sym}^d_{\succ 0}$ (symmetric positive definite). This problem is fundamental in numerical linear algebra, appearing in spectral equivalence, generalized singular value decompositions, and stability analysis of partial differential equations (via the numerical abscissa).

Conventional methods for solving this problem typically rely on matrix-adjoint products ($A^T$) and the explicit computation of $B^{-1}$. The authors identify two primary limitations with these approaches:

1. Computational Cost and Stability: Computing $B^{-1}$ or performing matrix-vector products involving it is expensive and prone to numerical inaccuracies.
2. Adjoint Mismatch: In applications like computed tomography, the transpose $A^T$ is often replaced by an approximate, computationally tractable operator. This "adjoint mismatch" leads to significant reconstruction errors, rendering methods requiring exact adjoints unsuitable.

Methodology
The authors propose a stochastic zeroth-order Riemannian algorithm that operates directly on the manifold $S^{d-1}_B = \{v \in \mathbb{R}^d : \langle v, Bv \rangle = 1\}$. The method avoids computing $A^T$ or $B^{-1}$ entirely, relying solely on matrix-vector products with $A$ and $B$.

The core algorithm (Algorithm 1 for single-sample, Algorithm 2 for multi-sample) proceeds as follows:

1. Initialization: Start with a random vector $v_0$ on the manifold.
2. Sampling: At each iteration $k$, sample a random direction $x_k$ from the tangent space $T_{v_k} S^{d-1}_B$. This is achieved by projecting a standard Gaussian vector onto the tangent space and normalizing.
3. Optimal Step Size: Instead of using a fixed or backtracking step size, the algorithm analytically solves the one-dimensional maximization problem along the sampled direction:
   $$\tau_k = \arg\max_{\tau \in \mathbb{R}} r(A, B, v_k + \tau x_k)$$
   Theorem 3.1 provides a closed-form solution for $\tau_k$ involving only inner products of $v_k$ and $x_k$ with $A$ and $B$.
4. Update: The iterate is updated via the retraction map: $v_{k+1} = R_{v_k}(\tau_k x_k)$.

The multi-sample variant (Algorithm 2) averages multiple sampled directions to reduce the variance of the gradient estimate before computing the step size.

Theoretical Contributions
The paper establishes rigorous convergence guarantees under the condition that the dimension of the eigenspace corresponding to the maximal eigenvalue, $\dim(G_{\max})$, is less than $d-1$.

• Non-Termination: The authors prove that if $\dim(G_{\max}) < d-1$, the algorithm does not terminate prematurely (i.e., the gradient proxy $b_k$ does not vanish) with probability one.
• Convergence to Global Maximizers: The sequence of objective values $r(A, B, v_k)$ converges almost surely to the global maximum $R(A, B)$. Furthermore, the sequence of iterates converges to the set of global maximizers $G_{\max}$.
• Convergence Rates:
  • The norm of the Riemannian gradient vanishes at a sublinear rate of $O(1/n)$ with probability one.
  • Specifically, $\min_{k=0,\dots,n} \| \text{grad } f(v_k) \|^2 \leq \frac{C}{n+1}$, where $C$ depends on the condition number of $B$ and the norm of $A$.
• Zeroth-Order Interpretation: The method is reinterpreted as a stochastic zeroth-order Riemannian gradient ascent. The authors show that the update direction $b_k x_k$ serves as an unbiased estimator of the Riemannian gradient (scaled by $1/(d-1)$), linking the proposed method to existing stochastic optimization frameworks on manifolds.

Numerical Results
Extensive numerical experiments validate the theoretical findings and compare the proposed method against state-of-the-art alternatives:

• Comparison with ZO-RGA: The proposed method significantly outperforms stochastic zeroth-order Riemannian gradient ascent (ZO-RGA) from [32], particularly in high dimensions ($d=100, 500$). The analytical selection of the step size $\tau_k$ proves superior to constant step sizes or Armijo backtracking used in ZO-RGA.
• Comparison with Gen-Oja: In a Karhunen-Loève problem, the proposed method converges faster and yields more accurate eigenfunction approximations than both deterministic and stochastic Gen-Oja methods.
• Robustness: The algorithm demonstrates robustness against ill-conditioned matrices $B$, though convergence time increases with the condition number $\kappa(B)$.
• Sample Size ($m$): Increasing the number of samples $m$ (from 1 to 100) significantly improves convergence speed and reduces variance, with minimal additional runtime cost due to vectorization.
• Rayleigh-Ritz Subproblem: The authors compare their method to a Rayleigh-Ritz approach on the subspace spanned by the current iterate and $m$ samples. While the Rayleigh-Ritz method can be faster for specific large $m$ and $d$ combinations, the proposed method offers a better trade-off between computational complexity ($O(dm)$ vs. $O(d^2m)$) and convergence speed.

Significance and Claims
The paper claims to introduce a "simple and efficient" algorithm that solves the generalized Rayleigh quotient maximization without requiring adjoint or inverse computations. Its primary significance lies in:

1. Adjoint-Free Operation: It provides a provably convergent solution for problems where $A^T$ is unavailable or mismatched, a common scenario in imaging and inverse problems.
2. Theoretical Rigor: It provides the first convergence analysis for a stochastic zeroth-order Riemannian method applied specifically to the generalized Rayleigh quotient, establishing almost sure convergence to global maximizers and sublinear rates.
3. Performance: Empirical results demonstrate that the method outperforms existing stochastic and deterministic baselines, particularly in high-dimensional settings and when dealing with ill-conditioned matrices.

The work extends a previous preprint [41] by generalizing the framework to Riemannian geometry and generalized Rayleigh quotients, incorporating additional clarifications and substantial improvements in convergence analysis.