Technical Summary: Discrete Double-Bracket Flows for Isotropic-Noise Invariant Eigendecomposition

1. Problem Statement

The paper addresses the problem of full-basis matrix-free eigendecomposition for a symmetric operator under streaming observations. The observation model is given by:
$$C_k = C_{\text{sig}} + \sigma_k^2 I + E_k$$
where:

• $C_{\text{sig}} \in \mathbb{R}^{n \times n}$ is the signal covariance with distinct eigenvalues.
• $\sigma_k^2 I$ represents a time-varying, potentially unbounded isotropic background (noise floor).
• $E_k$ is an anisotropic perturbation.
• Access is restricted to matrix-vector products (MVPs) $v \mapsto C_k v$, meaning the full matrix $C_k$ is never explicitly stored or formed.

The Core Challenge: In high-background regimes where $\sigma_k^2 \gg \|C_{\text{sig}}\|_2$, standard algorithms (e.g., subspace iteration, Oja's rule, Euclidean gradient descent with QR retraction) suffer from severe degradation. Their stability, step sizes, and contraction rates become coupled to the noise floor $\sigma_k^2$. For instance, the contraction ratio in subspace iteration approaches unity as $\sigma_k^2 \to \infty$, and the effective step size in QR-Oja collapses to zero. This forces iteration counts to scale as $\Omega(\sigma_k^2 / \Delta)$, rendering them impractical for high-noise environments.

2. Methodology

The authors propose a Discrete Double-Bracket Flow that achieves algebraic invariance to isotropic shifts without explicitly estimating or subtracting the noise.

2.1 Algebraic Invariance via Lie Brackets

The central insight is that the isotropic component $\sigma_k^2 I$ lies in the center of the matrix algebra. Consequently, it vanishes under the Lie bracket operation (commutator). The authors construct a state-dependent generator $\Omega(M)$ defined on the Lie algebra $\mathfrak{so}(n)$:
$$\Omega(M) = [A(M), D(M)] = A(M)D(M) - D(M)A(M)$$
where:

• $A(M) = M^\top C_k M$
• $D(M) = \text{diag}(A(M))$

Because $[\alpha I, X] \equiv 0$ for any scalar $\alpha$ and matrix $X$, the generator satisfies:
$$\Omega(M; C_{\text{sig}} + \sigma_k^2 I + E_k) = \Omega(M; C_{\text{sig}} + E_k)$$
This ensures the update direction is pointwise identical whether the observation includes the isotropic noise or not. The noise is not filtered; it is structurally annihilated by the generator's definition.

2.2 Discrete Update via Cayley Retraction

The algorithm updates the rotation matrix $M_k \in SO(n)$ using the Cayley map to maintain orthogonality without expensive SVD or QR decompositions:
$$M_{k+1} = M_k \cdot \text{Cay}(\eta_k \Omega_k)$$
where $\text{Cay}(X) = (I - \frac{1}{2}X)^{-1}(I + \frac{1}{2}X)$.
In matrix-free settings, the inverse $(I - \frac{\eta}{2}\Omega)^{-1}$ is approximated via a truncated Neumann series of order $K$ (typically $K=3$ or $4$), requiring only $O(n)$ MVPs per iteration.

2.3 Extension to Top-$k$ Tracking

The framework is extended to the Stiefel manifold $\text{St}(k, n)$ for tracking the top-$k$ eigenvectors. By introducing a weighted Brockett objective $J(M) = \text{tr}(M^\top C M N)$ with a strictly ordered diagonal weight matrix $N$, the method isolates an ordered eigenframe while preserving the same $\sigma^2$-invariance property.

3. Key Contributions and Theoretical Results

3.1 Algebraic Invariance and Stability

• Pointwise Invariance: The discrete trajectory $\{M_k\}$ and Lyapunov values are identical for observations $\{C_{\text{sig}} + \sigma_k^2 I + E_k\}$ and $\{C_{\text{sig}} + E_k\}$.
• Stability Threshold: The maximal stable step size $\eta_{\max} = 1/L_C$ depends exclusively on the trace-free signal energy $\|C_e\|_2$ (where $C_e = \text{tf}(C_{\text{sig}})$), not on $\sigma_k^2$. Specifically, $L_C = (4\sqrt{n} + 8)\|C_e\|_2^2$.

3.2 Input-to-State Stability (ISS)

The authors establish that the algorithm is Input-to-State Stable with respect to the trace-free perturbation $E_e$.

• Noise Ball: The steady-state error bound (noise ball radius) is governed solely by $\|E_e\|$, $\|C_e\|$, and the spectral gap $\delta$. It is independent of the isotropic noise level $\{\sigma_k^2\}$.
• Sample Complexity: Under stochastic noise, the convergence rate is $O(1/k)$, and the sample complexity to reach accuracy $\epsilon$ is $O(\|C_e\|_2^2 \sigma_u^2 / (\delta^2 \epsilon))$, where $\sigma_u^2$ depends only on anisotropic noise moments.

3.3 Global Convergence

• Strict-Saddle Geometry: The objective function on $SO(n)$ has no spurious local minima. All non-optimal critical points are strict saddles.
• Haar-Almost-Sure Convergence: For the exact (noiseless) Cayley iteration, random initialization (Haar measure) converges to the global minimum set with probability 1.
• Finite-Time Entry: Under bounded trace-free perturbations, the perturbed orbit enters the spectrally separated domain (where local ISS analysis holds) in finite time from almost any initialization.

3.4 Computational Complexity

• **Full Basis ($SO(n)$):**$O(n)$ MVPs and $O(n^3)$ dense arithmetic per step.
• Top-$k$ ($\text{St}(k, n)$): $O(k)$ MVPs and $O(nk^2 + k^3)$ arithmetic per step, enabling scalability to $n \sim 10^4$ for small $k$.

4. Experimental Results

The paper validates theoretical predictions through numerical diagnostics rather than broad empirical benchmarks:

1. Iteration Complexity vs. Noise: The proposed Cayley method requires a constant number of iterations ($\approx 1600$) across four orders of magnitude of isotropic noise ($\sigma^2 \in [0, 1000]$). In contrast, "Raw Oja" fails completely for $\sigma^2 \ge 10$, and trace-removal baselines (TF-Oja) incur significant overhead.
2. Algebraic Correctness: The commutator implementation satisfies algebraic identities (e.g., Givens escape profiles, input bounds) to machine precision ($\sim 10^{-15}$).
3. Global Convergence: From 100 Haar-random initializations on $SO(15)$, the algorithm converges to the global minimum in 100% of trials.
4. Stiefel Invariance: Trajectories on $\text{St}(5, 50)$ with $\sigma^2 = 10^4$ are indistinguishable from those with $\sigma^2 = 0$ (difference $< 2.3 \times 10^{-12}$).

5. Significance and Claims

The paper claims to resolve a fundamental limitation in streaming eigendecomposition: the coupling of algorithmic stability to the noise floor. By leveraging the Lie bracket identity $[\alpha I, \cdot] \equiv 0$, the method renders isotropic noise "algebraically invisible."

• Structural vs. Estimation: Unlike trace estimation or denoising approaches, this method does not attempt to estimate or subtract the nuisance parameter. Instead, it constructs dynamics where the noise cannot exist.
• Robustness: The approach provides certified convergence guarantees even when the isotropic background is unbounded or impulsive.
• Scope: The authors explicitly limit their scope to isotropic shifts. Structured (non-isotropic) perturbations are outside the current framework. The full-basis variant remains $O(n^3)$, limiting it to $n \lesssim 10^3$, though the Stiefel extension offers scalability for low-rank tracking.

The work establishes a new paradigm for robust spectral estimation where the noise model is handled via geometric structure rather than statistical estimation.