Generalization of Zeroth-Order Method for Quotients of Quadratic Functions

1. Problem Statement

The paper addresses the optimization of the generalized operator norm and the extremal generalized Rayleigh quotient under zeroth-order (derivative-free) constraints. Specifically, the goal is to solve:
$$\left\| \frac{A}{B} \right\| := \max_{v \in \mathbb{R}^d \setminus \{0\}} \frac{\|Av\|}{\|Bv\|}$$
where $A \in \mathbb{R}^{m \times d}$ and $B \in \mathbb{R}^{\ell \times d}$ are matrices accessible only via forward evaluations (black-box access).

Key Constraints & Challenges:

• Adjoint-Free/Matrix-Free: The adjoints $A^T$ and $B^T$ are unavailable or inconsistent with the forward operators (common in legacy codes, PDE solvers, and tomography).
• Well-posedness: $B$ must have full column rank ($\text{rank}(B)=d$) so that $B^TB \succ 0$.
• Limitations of Existing Methods: Recent state-of-the-art zeroth-order methods (e.g., Bresch et al. [7–9]) rely on tangent-space sampling on a generalized sphere ($S^{d-1}_{B^TB}$). This requires knowledge of the metric $B^TB$ or its square root to project search directions, which is infeasible in strict black-box settings.

2. Methodology

The paper proposes a fundamentally different approach based on unconstrained sampling on the entire unit sphere rather than restricted tangent spaces.

A. Unconstrained Sampling Strategy

Instead of sampling search directions $x_k$ from the tangent space of the generalized sphere (which requires $B^TB$), the method samples $x_k$ uniformly from the standard unit sphere $S^{d-1}$:
$$x_k = \frac{x}{\|x\|}, \quad x \sim \mathcal{N}(0, I_d)$$
This removes the need for metric-dependent projections or adjoint information.

B. Optimal Step Size via Closed-Form Solution

For a current iterate $v_k$ and search direction $x_k$, the algorithm seeks an optimal step size $\tau_k$ to maximize the objective along the retraction path. The objective function along the path becomes a quotient of quadratic functions:
$$s_k(\tau) = \frac{\|A(v_k + \tau x_k)\|^2}{\|B(v_k + \tau x_k)\|^2} = \frac{a_k + 2b_k\tau + c_k\tau^2}{d_k + 2e_k\tau + f_k\tau^2}$$
The coefficients ($a_k, \dots, f_k$) are computed using only forward evaluations of $A$ and $B$.

• Theoretical Breakthrough: The paper proves that finding the optimal $\tau_k$ is equivalent to solving a 2D generalized eigenvalue problem involving $2 \times 2$ Gram matrices $M_k$ and $N_k$.
• Result: A closed-form solution for $\tau_k$ is derived (Theorem 3.1), ensuring the step size is optimal for the specific random direction without line search heuristics.

C. Zeroth-Order Riemannian Optimization Framework

The method constructs estimators for Riemannian optimization quantities using only function evaluations:

1. Gradient Estimator: An unbiased estimator of the Riemannian gradient is constructed using the random direction $x_k$ and the optimal step size logic.
2. Hessian Estimator (Quasi-Newton): The paper derives an unbiased estimator for the Riemannian Hessian (Corollary 3.3) using fourth-order moments of the uniform distribution on the sphere. This enables a Zeroth-Order Quasi-Newton method (Algorithm 3).
3. Variance Reduction: By averaging $m$ independent samples (Algorithm 2), the variance of the gradient estimator is reduced, accelerating convergence.

3. Key Contributions

1. Unconstrained Sampling on $S^{d-1}$:
   The paper introduces a sampling scheme that operates on the full unit sphere rather than the tangent space of a generalized sphere. This eliminates the requirement for $B^TB$ or its adjoint, making the method applicable to strict black-box scenarios where previous methods failed.

2. Closed-Form Optimal Step Size:
   Unlike standard zeroth-order methods that use fixed or heuristic step sizes, this work derives an exact, closed-form solution for the optimal step size $\tau_k$ for the specific quadratic quotient objective. This transforms the subproblem into a solvable 2D generalized eigenvalue problem.

3. Unbiased Hessian Estimation:
   The authors provide a principled derivation of an unbiased estimator for the Riemannian Hessian in a derivative-free setting. This allows for the construction of a second-order zeroth-order algorithm (Quasi-Newton) that captures curvature information without explicit derivatives.

4. Theoretical Convergence Guarantees:

   • Almost Sure Convergence: The sequence of iterates converges almost surely to the leading generalized eigenspace (Theorem 3.7).
   • Critical Point Analysis: The paper proves that the sequence of gradient proxies ($\alpha_k$) converges to zero almost surely, ensuring the algorithm does not get stuck at non-optimal critical points (under mild spectral assumptions).
   • Well-Definiteness: It is proven that the denominator in the step size calculation remains non-zero almost surely, ensuring the algorithm is well-defined.

4. Numerical Results

Experiments were conducted on synthetic random Gaussian matrices with dimensions up to $d=1000$.

• Performance vs. State-of-the-Art:
  • In low dimensions, the proposed method is comparable to recent tangent-space sampling methods [9].
  • In high dimensions ($d \gg 10$), the unconstrained sphere sampling significantly outperforms tangent-space methods. This is attributed to the larger search space and the avoidance of expensive projection operations (matrix-matrix multiplications) required by tangent-space methods.
• Convergence Rates:
  • The method exhibits linear convergence for the Relative Quotient Error (RQE).
  • The minimal squared residual (MSQR) for the generalized eigen-equation shows a convergence rate of $O(1/k)$, which is superior to the sublinear rates often seen in derivative-free methods.
• Second-Order Advantage:
  The Quasi-Newton variant (Algorithm 3) demonstrates stricter convergence behavior and faster error reduction near critical points compared to the first-order variant (Algorithm 2), effectively closing the performance gap near optima.

5. Significance

This work represents a significant advancement in derivative-free optimization for spectral problems.

• Practical Impact: It enables the computation of operator norms and Rayleigh quotients in scenarios where adjoint operators are unavailable, inconsistent, or too costly to compute (e.g., legacy simulation codes, complex PDE solvers).
• Theoretical Impact: It bridges the gap between zeroth-order random search and Riemannian optimization by providing unbiased estimators for both gradients and Hessians on manifolds without explicit manifold projections.
• Algorithmic Efficiency: The closed-form step size and the ability to operate without tangent space projections make the algorithm computationally efficient and scalable to high-dimensional problems.

In summary, Bresch generalizes zeroth-order methods to handle quotients of quadratic functions by replacing constrained tangent-space sampling with unconstrained sphere sampling, deriving exact step sizes, and establishing a rigorous second-order framework that outperforms existing state-of-the-art techniques in high-dimensional, black-box settings.