Stochastic Rayleigh-Ritz Ascent: Global Linear Convergence under a Spectral Gap

AMS Classification: 65F35, 15A60, 68W20

Introduction

Estimating the operator norm of a matrix $A$ amounts to maximizing the Rayleigh quotient of $A^{*}A$ over the unit sphere

Although the objective can be evaluated using only forward applications of $A$, classical algorithms for finding the maximizer, such as the power method [20, 10, 24] or Riemannian gradient ascent [2], require access to the adjoint $A^{*}$. An exact adjoint implementation, however, may be unavailable or computationally expensive in certain scenarios. For instance, in imaging applications such as computed tomography [9], the adjoint is often replaced by a computationally tractable approximation and discretization scheme to reduce memory use and improve computational speed on GPUs [1, 14, 16, 15]. The resulting adjoint mismatch can adversely affect reconstruction quality [28, 30, 25, 7, 6]. For operator-norm estimation, replacing $A^*$ by an approximation or a discretization, say $V^* \approx A^*$, also changes the operator used in the iteration, cf. [7 Sec. 1.1] and [6 Sec. 3.3], and therefore may be misleading for a successful application of power-schemes [27 Sec. 28 & 48]. More precisely, if $V^*A$ is non-normal, power-schemes may fail to converge or may return estimates far from the leading singular value. A full understanding of the impact of $V^*$ requires further analysis [27 p. 446], which should be bypassed by avoiding the use of either $V^*$ or $A^*$. These considerations motivate adjoint-free methods that estimate the norm of $A$ using only forward evaluations so that they can be applied to black-box implementations while retaining provable convergence guarantees.

Random-plane optimization provides a simple approach for solving (1) without involving the adjoint $A^*$. In [7], the authors proposed a method that at each iteration samples a direction in a tangent plane to the sphere at the current iterate. Then, the objective is exactly maximized over the unit circle in the plane spanned by the current iterate and the sampled direction. The resulting subproblem can be solved using only applications of $A$ and inner products. The detailed description of the algorithm can be found in § 1.1. This approach was generalized to the estimation of adjoint mismatch between two linear operators in [6] and extended to estimate generalized Rayleigh quotients in [8]. Together, these works establish sublinear convergence guarantees and residual estimates, while the numerical trials therein indicate much faster linear convergence rates. Results in [8] show that the random-plane iteration can be viewed as a Riemannian zeroth-order method with an optimized step size, and propose a variance reduction through multi-sample gradient estimators. Given this connection to zeroth-order optimization methods, we explore the connections to other algorithms in the literature, indicate their linear convergence guarantees, and discuss their differences.

Riemannian zeroth-order gradient methods [21] with constant step size are shown to converge sublinearly to the critical points of non-convex functions. Their Euclidean counterparts, e.g., in [17, 19], exhibit stronger (local) linear rates under the (local) Polyak–Łojasiewicz (PL) condition. In the context of maximizing the nonconstant Rayleigh quotient (1), nonleading eigenvectors are critical points that are not globally optimal, limiting the performance of the mentioned general guarantees. Inexact Riemannian gradient descent [31] provides another possible framework through controlled gradient errors and Kurdyka–Łojasiewicz (KL) inequalities. A single random tangent direction does not, in general, satisfy the uniform relative-accuracy requirement of that analysis, while its alternative absolute-error conditions would require separate verification.

In contrast, first-order methods analyzed specifically for eigenvalue problem (1) provide more insight into the local geometry. Riemannian steepest ascent [2, 4, 3] and the randomized submanifold method [13] converge linearly under the Riemannian PL condition. In particular, in [2] the rate of convergence is directly linked to the spectral gap of $A^*A$. However, these methods use the adjoint to perform updates and are not available in the setting of this paper.

Motivated by these observations, we establish explicit local and global linear rates for the expected objective error under a spectral-gap assumption. Our local rate estimate, which depends on the spectral gap, spectral width, and dimension, is near-optimal up to a constant multiplicative factor that corresponds to the size of the convergence neighborhood.

The paper is structured in the following way. The preliminaries in § 1.1 cover necessary background on the random-plane algorithm for estimating the operator norm [7]. § 1.2 states the main result. § 2.1 establishes the spectral inequalities and their connection to local expected linear convergence. We discuss results of § 2.1 in the context of the PL condition in § 2.2. § 2.3 derives the global linear convergence in expectation. In § 2.4, we extend the analysis to the zeroth-order estimation for the generalized Rayleigh quotient considered in [8]. We compare our theoretical findings with the empirical performance of the random-plane algorithm in § 3, followed by a conclusion.

Preliminaries

Algorithm. We consider the problem of finding the eigenvector corresponding to the leading eigenvalue of $A^{*}A$ for a given matrix $A \in\mathbb{R}^{m \times d}$ without involving the adjoint $A^{*}$. This can be achieved through the stochastic random-plane algorithm [7], Alg. 1. For our analysis, we focus on its Rayleigh-Ritz reformulation discussed in [8], Sec. 8.6. That is, let $v^{0} \in\mathbb{S}^{d-1} := \{x \in\mathbb{R}^{d} \mid\lVert x \rVert= 1\}$ be the initial guess sampled uniformly at random, with $\lVert\cdot\rVert$ denoting the Euclidean norm in $\mathbb{R}^{d}$. Then, the sequence of iterates $(v^{k})_{k \ge0}$ is constructed by iteratively solving

where the direction $x^{k}$ is sampled uniformly at random on the unit sphere in the tangent space $T_{v^{k}} = \{x \in\mathbb{R}^{d} \mid\langle x,v^{k}\rangle= 0\}$ to the sphere $\mathbb{S}^{d-1}$ at the point $v^{k}$. This is achieved by sampling $y \sim\mathcal{N}(0,I_{d})$, the standard Gaussian random vector, and setting

where $P_{T_{v^{k}}} := I_{d} - v^{k}(v^{k})^{*}$ is the orthogonal projector onto the tangent space and $I_{d}$ is the $d \times d$ identity matrix.

Probability Space and Distributions. We can associate the samples $(v^{k})_{k \ge0}$ and $(x^{k})_{k \ge0}$ to the sequences of random variables. Following [8], Rem. 4.1, we consider a probability space $(\Omega,\mathcal{A},\mathbb{P})$ and a sequence of random variables $V^{k}:\Omega\to\mathbb{S}^{d-1}$ with law $\mathbb{P}_{V^{k}} := V^{k}_{\#}\mathbb{P} = \mathbb{P}(V^{k})^{-1}$. By construction, $\mathbb{P}_{V^{0}} = \sigma_{\mathbb{S}^{d-1}}$, the Haar measure on the $(d-1)$-dimensional unit sphere and $v^{0} \sim V^{0}$. Moreover, for each sample $v^{k}$, we consider random variables $X^{k}:\Omega\to T_{v^{k}} \cap\mathbb{S}^{d-1} \simeq\mathbb{S}^{d-2}$ almost surely. Since they depend on $v^{k}$, by [8], Lem. 4.1, their conditional distribution given $V^{k}=v^{k}$ is uniform, i.e., $\mathbb{P}_{X^{k}\mid V^{k}=v^{k}} = \mathcal{U}(T_{v^{k}} \cap\mathbb{S}^{d-1}) = \sigma_{\mathbb{S}^{d-2}}$. As a result, we can view $(V^{k})_{k\in\mathbb{N}}$ as a Markov chain with values on $\mathbb{S}^{d-1}$. Furthermore, we say that something holds almost surely (a.s.) if its exceptional event has probability zero.

Solutions of (2). Since $v^k$ and $x^k$ are orthogonal by construction, any $z \in\operatorname{span}\{v^k,x^k\}$ with $\lVert z\rVert=1$ can be written as $z=u_1v^k+u_2x^k$ with $\lVert u\rVert=1$. The subproblem (2) then corresponds to finding the leading eigenvalue $\rho_{k+1}$ and eigenvector $u\in\mathbb{R}^2$ of the system matrix

see [5 Rem. 3.4] for further details.

If $a_k=0$, $H_k$ is a diagonal matrix and the maximizer is determined by the larger of the values $\rho_k$ and $s_k$,

Otherwise, if $a_k\ne0$, according to [7 Prop. 2.4], the unique maximizer of (2) is given by

In all cases, we have

and the algorithm produces a nondecreasing sequence $(\rho_k)_{k\ge0}$.

Sublinear Convergence. To put the generated sequence $(v^k)_{k\ge0}$ in perspective of finding the leading eigenvector of $A^*A$ for $A\in\mathbb{R}^{m\times d}$, let us denote the eigenvalues of $A^*A\in\mathbb{R}^{d\times d}$ by $\lambda_1\ge\lambda_2\ge\cdots\ge\lambda_d\ge0$. Throughout the paper, the eigenvectors $u$ are always unit vectors, i.e., $\lVert u\rVert=1$. We define the spectral error as

and the residual at $k$th iteration as

The main results of [8] provide the following convergence guarantees for $\varepsilon_k$ and $r^k$.

Main Contributions

While Thm. 1.1 establishes almost-sure convergence of the spectral error $\varepsilon_k$ and the iterates (2) and a sublinear bound for the expected squared residual, it does not quantify the decay of $\mathbb{E}[\varepsilon_k]$. In contrast, other algorithms for finding the leading eigenvector, such as the (randomized) power method [23, 20, 12] or Riemannian gradient ascent [4, 29], converge at a linear rate when the spectral gap is positive. Our main contribution lies in closing this gap and establishing global linear convergence guarantees in expectation under a positive spectral gap, summarized in the following main theorem.

The proof combines two estimates: First, for any $\eta\in(0,1)$, we establish local expected contraction

in the neighborhood of global maximum described by the set $\{\varepsilon_k\le\eta\delta\}$ with $0<\eta<1$. Second, we derive a geometric tail bound for the first entrance time into this region. Combining these estimates yields the global rate from iteration zero, see Thm. 2.3, 2.7, and 2.8. Then, we derive convergence rates for the squared residual norm $r^k$ and the squared distance to the leading eigenspace by bounding these values in terms of $\varepsilon_k$ in Cor. 2.9.

Linear Convergence under a Spectral Gap

Throughout the remainder of the paper, we rely on assumptions of Thm. 1.2. The case $d\le2$ is excluded, since by Theorem 1.1, the algorithm will terminate after at most one iteration a.s.

Since we are interested in a linear rate of

by (5) and (6), we first consider the local linear convergence in the neighborhood of the solution in § 2.1 and then study the probability of entering it in § 2.3.

Ascent Bounds and Local Contraction

In this subsection we first derive qualitative bounds on the expected error. We start with an analogue of sufficient decrease condition of $\varepsilon_k$ in terms of $\lVert r^k\rVert^2$.

Proof. Utilizing (4) and characteristic equation gives

Both $\rho_{k+1}$ and $s_k$ belong to $[\lambda_d,\lambda_1]$, and $\rho_{k+1}\ge\max\{\rho_k,s_k\}$. Consequently, by (6), we get

Now, we transform $a_k^2$. By (4), we have

since $x^k\perp v^k$ by construction in (3). Next, [8 Lem. 4.1] gives

which yields

Finally, we notice that (7) implies

Consequently, taking the conditional expectation in (10) and combining it with (13) and (14) proves the assertion. □

Next, we relate $\lVert r^k\rVert^2$ to the errors $\varepsilon_k$.

Proof. Let $(u_i)_{i=1}^{d}$ be an orthonormal eigenbasis of $A^*A$ and write $v^k=\sum_{i=1}^{d}\gamma_i^k u_i$ for some $\{\gamma_i^k\}_{i=1}^{d}\subset\mathbb{R}$. Then

Set $\delta_i:=\lambda_1-\lambda_i$, so that $\delta_1=0$, $\delta_2=\delta$, and $\delta_i\ge\delta$ for $i\ge2$. We have $\varepsilon_k$ as

On the other hand, we have

using that $\{u_i\}_{i=1}^{d}$ is a basis. Therefore,

Finally, by using $\delta_1^2=0=\delta\delta_1$ and $\delta_i\ge\delta>0$ for $i\ge2$, we obtain $\delta_i^2\ge\delta_i\delta\ge0$, and the assertion follows by bounding

The right hand side, follows by using $\delta_i\le\Delta$,

In the following we consider neighborhood of the global maximum, where errors satisfies $\varepsilon_k\le\eta\delta$ for some $\eta\in(0,1)$. This condition is defined by the set

Next, we show that within this set we obtain a linear convergence rate.

Proof. By Lem. 2.1 and 2.2 we obtain

For $v^k\in\mathcal{E}_{\eta,\delta}$ we have $\delta-\varepsilon_k\ge(1-\eta)\delta$, which proves (19). Since $(\varepsilon_k)_{k\ge0}$ is nonincreasing, this event also holds for subsequent iterations. Equation (20) follows by induction and the tower property,

Choosing $\eta=\frac{1}{2}$ gives the conditional contraction factor

once the iterates $v^k$ satisfy $\lVert Av^k\rVert^2\ge(\lambda_1+\lambda_2)/2$. The relevant gap is therefore the gap of $A^*A$, normalized by its spectral width. Since $\Delta\le\lambda_1$, a simpler, possibly weaker factor is

These are upper bounds on the expected error, and do not assert a deterministic contraction at each realization of a random step. Moreover, the obtained rate is tight up to the multiplier $1-\eta$, as indicated by the next example.

Its leading eigenvalue is $\rho_{k+1} = \lambda_1 - \delta+ \delta(\alpha_k^2 + \beta_k^2)$. Hence the exact ascent is $\epsilon_k - \epsilon_{k+1} = \rho_{k+1} - \rho_k = \delta\beta_k^2$. By utilizing (12), we obtain

and, consequently, by iterating we get

Compared with Thm. 2.3, this rate matches the rate without the $1-\eta$ multiplier, i.e. using $\eta=0$. Therefore, our local linear convergence rate is near-optimal and we refer to the optimal local linear convergence rate as $q_0 := q(0,\delta,\Delta)=1-\frac{\delta}{\Delta(d-1)}$.

The dimension dependence of our convergence guarantee contrasts with that of the randomized power method applied to $A^*A$, where the initialization is uniformly sampled on the unit sphere [20, 10]. Under the assumption of a positive spectral gap, the contraction factor becomes $\lambda_2/\lambda_1=1-\delta/\lambda_1$, without an extra $d-1$ factor in the denominator. In view of Ex. 2.4, the random plane algorithm encounters the curse of dimensionality. Whether a different adjoint-free algorithm may avoid it is an open problem.

Interpretation through the PL Condition

Thm. 2.3 admits an interpretation through the Polyak–Lojasiewicz (PL) condition [26, 18]. A differentiable function $H$ on a Riemannian manifold satisfies a PL inequality on $\mathcal{U}$, relative to a minimum value $H_*$, if

where $H(v)\ge H_*$ on $\mathcal{U}$ and the gradient norm is induced by the Riemannian metric; see also [13 Definition 5.4]. This gradient-domination condition is implied by Euclidean strong convexity, but requires neither convexity nor a unique minimizer [18].

For $\epsilon(v):=\lambda_1-\|Av\|^2$ on $\mathbb{S}^{d-1}$ as in (6), we have $\min_{\|v\|=1}\epsilon(v)=0$ and $\operatorname{grad}\epsilon(v)=-2r(v)$, where $r(v):=A^*Av-\|Av\|^2v$ is defined in (7). Lem. 2.2 therefore gives, for $\eta\in(0,1)$ and $v\in\mathcal{E}_{\eta,\delta}$,

Thus, $\epsilon$ satisfies a PL inequality on $\mathcal{E}_{\eta,\delta}$ with $\mu=\mu_{\eta,\delta}=2(1-\eta)\delta$, without requiring geodesic strong convexity. Combining this inequality with Lem. 2.1 yields

where $q$ is defined in (18).

Gradient domination was also used to establish local linear convergence rates for Riemannian gradient ascent for finding the leading eigenvector in [2 Prop. 4] under the assumption of the spectral gap and with the exact logarithmic map on the sphere. However, a global PL inequality relative to $\min_{\lVert v\rVert=1}\varepsilon(v)=0$ cannot hold since it implies that the function only has global minima. This condition fails since nonleading eigenvectors are critical points with $\varepsilon(v)>0$. Moreover, the local PL inequality does not quantify the time needed to reach $\mathcal{E}_{\eta,\delta}$. Sec. 2.3 addresses this remaining step.

Global Linear Convergence Rates

Now, we focus on entering the neighborhood $\mathcal{E}_{\eta,\delta}$. The almost sure convergence rate established in Thm. 1.1 implies that the number of iterations needed

is almost surely finite. We first study the tails of its distribution.

Proof. By Bayes rule, we have

Since $V^k\in\mathcal{E}_{\eta,\delta}$ is equivalent to $T_{\eta,\delta}\le k$, the first integral is zero. Next, we focus on the complement

for $v^k\notin\mathcal{E}_{\eta,\delta}$ characterized by the random variable $X^k$. Let $u_1\in\mathbb{S}^{d-1}$ be the eigenvector of $A^*A$ corresponding to $\lambda_1$. We write $\alpha_k:=\langle v^k,u_1\rangle$, and define

Notice that the denominator is nonzero, since $|\alpha_k|=1$ would imply $\epsilon_k=0$ and contradict $v^k\notin\mathcal{E}_{\eta,\delta}$. Furthermore, let $P_{\operatorname{span}\{v^k,x^k\}}$ be the orthogonal projection onto $\operatorname{span}\{v^k,x^k\}$ such that

Orthogonality gives

and by $u_1=\alpha_kv^k+\sqrt{1-\alpha_k^2}w^k$ we obtain $\langle u_1,x^k\rangle=\sqrt{1-\alpha_k^2}\langle w^k,x^k\rangle$ such that

since $|\langle x^k,w^k\rangle|\leq1$ by the Cauchy-Schwarz inequality. On the event

the vector $z:=P_{\operatorname{span}\{v^k,x^k\}}u_1/\|P_{\operatorname{span}\{v^k,x^k\}}u_1\|$ is well defined and satisfies

Since $\|z\|=1$, expanding $z$ in an orthonormal eigenbasis $\{u_i\}_{i=1}^d$ of $A^*A$ gives

Hence, by optimality from (2), on $G_k$ it holds

Therefore,

and it remains to bound the conditional probability of $X^k\in G_k$. Recall that given $V^k=v^k$ the random vector $X^k$ is uniformly distributed on the sphere $T_{v^k}\cap\mathbb{S}^{d-1}\cong\mathbb{S}^{d-2}$. Since $w^k\in T_{v^k}\cap\mathbb{S}^{d-1}$, the event $G_k=\{|\langle X^k,w^k\rangle|^2\geq1-\frac{\eta\delta}{\Delta}\}$ consists of two antipodal spherical caps with polar angle $\theta=\arccos((1-\frac{\eta\delta}{\Delta})^{1/2})$. According to [11, 22], their normalized surface area is given by the regularized incomplete Beta function $I_x(a,b)$ evaluated at $x=\sin^2(\theta)=1-(1-\frac{\eta\delta}{\Delta})=\frac{\eta\delta}{\Delta}$, $a=\frac{d-2}{2}$ and $b=\frac{1}{2}$, yielding

By [7 Lem. A.1], we obtain the more restrictive bound

which proves (21). Now, we return to (23) and substitute (24) and (27) to get

Repeating the bound $k$ times gives

where we used the fact that the probability is at most one.

Remark 2.6. The initial survival probability can be bounded more sharply. In the last step of Prop. 2.5, we derive the stronger upper bound on $\mathbb{P}(T_{\eta,\delta}>0)$. For $V^{0}\sim U(\mathbb{S}^{d-1})$, we have by (26) that

The event on the left consists of two antipodal spherical caps, now, on the $(d-1)$-dimensional sphere. Hence, by (24), we obtain

Now, we are ready to establish the linear convergence rates from the main theorem, cf. Thm. 1.2.

Proof. For $0\le j\le k$, we split

By $\varepsilon_k \le\Delta$ and Prop. 2.5, the first term is at most

The event $\{T_{\eta,\delta}\le j\}$ implies $\varepsilon_j\le\eta\delta$. Therefore, taking the expectation over $X_j,\ldots,X_k$ first by Thm. 2.3 and (20) in particular gives

This proves (28). Now, we would ideally choose $j$ such that $(1-p)^j=q^{k-j}$. This corresponds to $(q(1-p))^j=q^k$ and taking the logarithm gives optimal $j_\ast$

Since this value is not necessarily an integer, we pick $j=\lceil j_\ast\rceil$, the rounded-up value. As a result, we have

where we used that $\lceil j_\ast\rceil-1\le j_\ast\le\lceil j_\ast\rceil$. Finally, we observe that

which proves (30). $\square$

We also establish an alternative linear rate in the case when $q\ne1-p$.

Proof. In contrast to Thm. 2.7, we perform a more detailed split on $T_{\eta,\delta}$, more precisely

The first term is bounded as in (31). The second term can be bounded analogously to (32) by $\mathbb{E}[\varepsilon_k\mathbf{1}_{\{T_{\eta,\delta}=t\}}]\le\eta\delta q^{k-t}\mathbb{P}(T_{\eta,\delta}=t)$. Using

we obtain

If $q\ne(1-p)$, we get a geometric sum (35), which gives

Note that

which gives

Now, if $q=(1-p)=\varrho_2$, the sum is exactly $k$, yielding

Splitting $\varrho_2^k=\varrho_2^{k(1-\alpha)}\varrho_2^{k\alpha}$ for any $\alpha\in(0,1)$ gives

Consider the function $g(t)=t\varrho_2^{\alpha t}$. Since $g(t)\to0$, as $t\to\infty$, $g(0)=0$ and $g(t)\ge0$ on $[0,+\infty)$, the function is bounded on $[0,+\infty)$ and its maximum value is a critical point

Consequently, we bound

Substituting this bound into (36) concludes the proof.

The global bound accounts for the time needed to enter the superlevel set. Its cap probability $p(\eta,\delta,\Delta)$ can be small for large $d$, as we see an explicit dependence on the dimension in (21). Our numerical trials indicate that the lower bound on $p(\eta,\delta,\Delta)$ in (21) is significantly suboptimal and the rate is much faster.

Figure 1 shows the global rates $\varrho_1$ from Thm. 2.7 and $\varrho_2$ from Thm. 2.8 over the local convergence rate $q$ for various $p$. We observe that rates $\varrho_1$ and $\varrho_2$ are close if $\lvert1-p-q\rvert$ is large and $\varrho_2$ is better than $\varrho_1$ for small $\lvert1-p-q\rvert$.

Next, we establish bounds for the expected error and the distance to the leading eigenvector in terms of errors, which allows us to extend linear convergence rates to these quantities.

Therefore, combining Cor. 2.9 with Thm. 2.7 or Thm. 2.8 gives a linear convergence rate for the average norm of the residual and the average distance to the leading eigenvector.

Linear Convergence for Estimating the Generalized Rayleigh Quotient

In this section, we briefly discuss the extension of the derived linear convergence rates to the generalized-Rayleigh-quotient method proposed in [8]. We consider for the computation of the extremal generalized Rayleigh quotient

of $A\in\mathbb{R}^{d\times d}$ and positive definite $B=B^{*}\in\mathbb{R}^{d\times d}$. To solve (41), the authors of [8] consider the exact random-plane iteration

with

similar to (2). Hence, we obtain a random sequence $(v^{k})_{k\in\mathbb{N}}$ with $\langle v^{k},Bv^{k}\rangle=1$, where the latter maximization problem can be solved uniquely [8 Thm. 3.1], almost surely. Transforming $\widetilde{v}^{k}:=B^{1/2}v^{k}$, and constructing

the transformed search direction fulfills

and is generally not uniformly distributed. Moreover, it holds

Consequently, the generalized Rayleigh quotient (41) becomes the ordinary Rayleigh quotient of $C:=B^{-1/2}A^{\mathrm{H}}B^{-1/2}$, and the transformed iteration satisfies

Moreover, the eigenvalues of $C$ and the generalized eigenvalues of $(A^{\mathrm{H}},B)$ given by

coincide. For that reason, we redefine in the manner of (4), (6), and (7)

and

Thereby, the system matrix $H_k$ from (4) is defined by $a_k := \langle C\tilde{v}^k,\tilde{x}^k\rangle$, and $s_k := \langle\tilde{x}^k,C\tilde{x}^k\rangle$. The proof of Lem. 2.1 carries over through (11). Using the orthogonal projection with respect to the changed Riemannian metric yields

Then, the conditional second moment can be, instead of (13), bounded by

To bound the numerator, we have

Together, we obtain

Combining this with (10) proves the generalized version of Lem. 2.1

The change of variables applies directly to the proof technique of Lem. 2.2, hence we obtain

Together with (42) we obtain for $\varepsilon_k\le\eta\delta$ that

which provides the generalized local linear rate constant (18).

Similar to Prop. 2.5, let $\tilde{u}_1$ be a leading unit eigenvector of $C$. The eigenbasis argument in (26) gives

and

Consequently, with $\varepsilon:= \sqrt{\eta\delta/\Delta}$ it holds

By [8 Lem. 5.5] the entrance probability becomes

Then, Theorems 2.7 and 2.8 apply unchanged with $q$ and $p$ replaced by $q_{B}$ and $p_{B}$, respectively.

Numerical Results

In this section, we substantiate the theoretical results by exploring the actual speed of the random plane iteration (2–3) from [7], compared with the derived convergence rates in Thm. 2.7 and 2.8.

Firstly, we study the local convergence rate in the neighborhood of the solution. Secondly, we analyze the global convergence behavior with respect to the entrance time and subsequent contraction. Finally, we study the influence of the spectrum $\{\lambda_i\}_{i=1}^{d}$ of $A^{*}A$ on the convergence.

All algorithms are implemented in Python and the code is publicly available⁴. The experiments are performed on an off-the-shelf MacBook Pro 2020 with Apple M1 (8-Core GPU, 3.2 GHz) and 16 GB RAM.

Controlled Spectra and Initialization

In this subsection, we analyze how the spectral gap, dimension, and initialization affect the convergence behavior and compare empirical and theoretical local rates.

Given a desired spectral gap $\delta\in[0,1/2]$, we generate matrices $A\in\mathbb{R}^{d\times d}$ via

with $Q$ being sampled uniformly on the orthogonal group $O(d)=\{V\in\mathbb{R}^{d\times d}\mid V^{*}V=VV^{*}=I_d\}$, denoted by $Q\sim\mathcal{U}(O(d))$. The spectrum is set to $\lambda_1=1$, $\lambda_2=1-\delta=\lambda_1-\delta$ and the lower eigenvalues evenly spaced in $[0,1/2]$, including $\lambda_d=0$. Consequently, $A^{*}A$ has spectral gap $\delta$, spectral width $\Delta=1$ and $\|A\|=1$.

Furthermore, for given $\eta\in[0,1]$ we initialize the algorithm at

where $u_1$ and $u_2$ are the eigenvectors of $A^*A$ corresponding to $\lambda_1$ and, respectively, $\lambda_2$. Hence, $\epsilon_0=\eta\delta/2\leq\eta\delta$, so the local convergence rate (18) applies from the first iteration,

Thus, the guaranteed contraction becomes weaker as the dimension increases, the gap $\delta$ decreases, or $\eta$ increases.

In our experiments, we fix one matrix for each spectral configuration

and perform 96 independent runs of the random-plane algorithm for each $\eta\in\{0.01,0.1,0.5\}$. All runs start from the deterministic initial vector (43)—randomness arises from the sampled search directions.

For each configuration, we use at most 2 500 iterations and otherwise four theoretical decay times, where one decay time corresponds to $1/(-\log q)$ iterations:

Figure 2 reports the mean of the relative errors $\epsilon_k/\epsilon_0$ over 96 runs as a function of iteration number and compare it with the optimal theoretical local convergence rate $q_0=q(0,\delta,\Delta)$ from Ex. 2.4.

The results show approximately straight mean-error lines on logarithmic vertical axes over the observed contraction ranges, consistent with linear decay rates. Larger spectral gaps $\delta$ produce faster convergence, while larger $d$ slows it. These observations support the theorem’s qualitative prediction of local linear convergence and its dependence on

the spectral gap and the dimension. The empirical rates, however, are faster than the actual theoretical bounds. Across all tested configurations, the observed rates exceed the theoretical rates by factors of approximately 1.9 – 6.3. In contrast to theoretical results, we did not observe a strong empirical dependence on $\eta$. For instance, at $d = 20$ and $\delta= 0.1$, the observed rate $-\log\hat{q}$ changes only from approximately 0.0163 to 0.0162 as $\eta$ increases from 0.01 to 0.5, while the theoretical rate $-\log q$ falls from 0.00522 to 0.00264.

Global convergence

In this section, we aim to study the time and computational cost required to enter the neighborhood of local convergence.

Therefore, we use the same setting as in Sec. 3.1, but now sample the initialization $v^0$ uniformly on the unit sphere. For each spectral configuration

we fix one matrix and perform 128 independent runs with an iteration limit $K = 3500$. The same initial vectors and random draws are reused across $\eta$, hence, only the threshold defining the local contraction region changes.

The first row of Figure 3 shows the empirical mean error over $N = 128$ independent runs, together with pointwise 95% bootstrap confidence intervals (shaded), the limiting local reference $q_0$, and the theoretical global bounds. The confidence intervals are computed

from 150 bootstrap resamples, each consisting of 128 complete trajectories sampled with replacement.

In the second row of Figure 3, we compare the empirical survival probability, estimated by the fraction of runs that have not yet entered the local contraction region, with this theoretical bound. By Prop. 2.5, the entrance-time survival probability is bounded by $(1-p(\eta,\delta,\Delta))^k$. The horizontal red line marks $u_N=1-0.051^{1/N}\approx0.0231$, the one-sided 95% binomial upper confidence limit for the survival probability when no survivors are observed among $N$ independent runs. For each spectral configuration, the red shaded band in the first row of Figure 3 spans the earliest and latest iterations, across $\eta$, at which the empirical survival curves fall below this reference level. These bands summarize the observed entrance behavior across the chosen thresholds; they are neither confidence intervals for the entrance time nor ranges at which 95% of runs have entered. The mean-error curves also illustrate changes in the observed decay around these entrance ranges.

For the tested configurations, entrance takes longer for larger dimensions, smaller spectral gaps, and smaller values of $\eta$, which define stricter local regions. tbl:1 summarizes the sample median entrance times. The global bounds established in Thm. 2.7 and Thm. 2.8 are substantially more conservative than the observed behavior. The local contraction factor $q(\eta,\delta,\Delta)$ provides a more informative description after entrance, indicating room for improvement in the global entrance-time analysis.

Effect of the Lower spectrum

Lastly, we study in this section the impact of the remaining spectrum beyond the leading and second leading singular values, i.e., the spectral gap, on the achieved convergence behavior in practice. Therefore, let us consider $\{A_{(i)}\}_{i=1}^{3}$ as in Sec. 3.1, with $\lambda_1=1$, $\delta=0.1$, i.e. $\lambda_2=\lambda_1-\delta$. For $A_1$, the remaining singular values spread over $[\sqrt{\lambda_2},\ldots,0]$ equidistantly, for $A_2$ they are clustered on the second-largest $\lambda_2=\ldots=\lambda_{d-1}=\lambda_1-\delta$, with $\lambda_d=0$ and for $A_3$ we cluster them on the smallest, i.e. $\lambda_3=\ldots=\lambda_d=0$. Therefore, all matrices have operator norm $\lVert A_{(i)}\rVert=1$, spectral gap $\delta$, and spectral width $\Delta=1$.

We compare two initialization settings for each of the three spectral configurations for $d=100$. In the first setting, we sample $v^0$ uniformly on the unit sphere $\mathbb{S}^{d-1}$ and perform 256 independent runs with 20 000 iterations per run. The mean spectral error relative to the initial error with the pointwise 95% percentile bootstrap confidence intervals, computed using 200 samples of the recorded trajectories, is visualized in cf. Figure 4 (left). In the second setting, for each $\eta\in\{0.01,0.1,0.5\}$, we sample random initial vectors inside the corresponding local contraction region $\mathcal{E}_{\eta,\delta}$ and perform 256 independent runs with 5 000 iterations per run. Specifically, we independently sample

and set

This construction gives $\epsilon_0=\upsilon\eta\delta\leq\eta\delta$ and $\mathbb{E}[\epsilon_0]=\eta\delta/2$. It defines a random distribution inside $\mathcal{E}_{\eta,\delta}$, rather than uniform surface measure on that region. For this setting we report the same as in the first setting in Figure 4 (middle). The observed ratio is depicted in Figure 4 (right), here, we consider $\eta\in\{0.01,0.1,0.5\}$.

The main observation in Figure 4 is that a zero eigenvalue with multiplicity $d-2$ produces the slowest convergence. For the spread spectrum, the rate is intermediate, whereas the configuration with a repeated second-largest eigenvalue yields the fastest observed convergence. Changing $\eta$ has a comparatively modest effect on the fitted empirical rates. Since the theoretical rate $q(\eta,\delta,\Delta)$ increases with $\eta$, the empirical-to-theoretical rate ratios increase with $\eta$, as illustrated in Figure 4. Table 2 reports the fitted rates $e^{-k\log\hat{q}}$, their ratios to the theoretical rate, and the final relative mean errors for $\eta=0$, i.e. compared with $q_0=q(0,\delta,\Delta)$.

Spectral configuration	Observed rate $-\log\hat{q}$	Ratio $\frac{\log\hat{q}}{\log q_0}$	Final mean / initial error
Spread spectrum	$5.905\times10^{-3}$	$5.84$	$5.597\times10^{-11}$
Clustered at $\lambda_2$	$9.731\times10^{-3}$	$9.63$	$5.464\times10^{-22}$
Clustered at zero	$5.657\times10^{-3}$	$5.60$	$7.441\times10^{-8}$

Conclusions

In this paper, we considered a problem of finding the largest singular value of a matrix without involving the adjoint matrix using the random-plane iteration algorithm. Complementing existing results on global almost sure convergence, we established local and global linear convergence rates in the presence of the spectral gap, with the local rate being near-optimal. The global rate depends strongly on the dimension of the underlying problem, and numerical experiments indicate that, in practice, a much faster rate applies and a better global rate could be a subject for future research. We also observe differences between local and global convergence behavior, showing that considering both scenarios is necessary. Interestingly, in the last numerical experiment, we saw how other eigenvalues affect the speed of convergence, which the theoretical analysis does not cover.