A Comprehensively Tight Analysis of Gradient Descent for PCA
Year: 2021

Abstract:
We study the Riemannian gradient method for PCA on which a crucial fact is that despite the simplicity of the considered setting, i.e., deterministic version of Krasulina’s method, the convergence rate has not been well-understood yet. In this work, we provide a general tight analysis for the gap-dependent rate at O( 1 ∆ log 1 ε ) that holds for any real symmetric matrix. More importantly, when the gap ∆ is significantly smaller than the target accuracy ε on the objective suboptimality of the final solution, the rate of this type is actually not tight any more, which calls for a worst-case rate. We further give the first worst-case analysis that achieves a rate of convergence at O( 1ε log 1 ε ). The two analyses naturally roll out a comprehensively tight convergence rate at O( 1 max{∆,ε}log 1 ε ). Particularly, our gap-dependent analysis suggests a new promising learning rate for stochastic variance reduced PCA algorithms. Experiments are conducted to confirm our findings as well.

text: 1 ϵ ) that holds for any real symmetric matrix. More importantly, when
the gap ∆ is significantly smaller than the target accuracy ϵ on the objective suboptimality of the final solution, the rate of this type is actually not tight any more, which calls for a worst-case rate. We further give the first worst-case analysis that achieves a rate of convergence at O( 1ϵ log 1 ϵ ). The two analyses naturally roll out a comprehensively tight convergence rate at O( 1max{∆,ϵ}log 1 ϵ ). Particularly, our gap-dependent analysis suggests a new promising learning rate for stochastic variance reduced PCA algorithms. Experiments are conducted to confirm our findings as well.

heading: 1 Introduction
text: Gradient descent is a basic, well-established, and celebrated optimization method for minimizing convex functions [11, 3], but remains far less understood theoretically for non-convex functions though it has been proven useful in practice for many non-convex problems [22, 17]. In the past decade, the research on gradient methods for non-convex problems has been gaining increasing attention in optimization and especially machine learning [6, 8, 4]. In this work, we focus on a notable example, namely the Riemannian gradient descent for the leading eigenvector computation of a real symmetric matrix A ∈ Rn×n, which lays the foundation of Krasulina’s method for PCA. This problem can be solved by either a traditional projection based method from the numerical algebra, such as the well-known power method or Lanczos algorithm [13], or a gradient search based method for the associated optimization problem such as the Riemannian gradient descent [18, 21]. The projected gradient descent (with constant step-sizes), as another closely related solver, can be interpreted either as a projection method or as a search method in this case (see Section 6 for interpretations). Compared to projection methods which have arguably been well-understood already [13, 9, 10], the search method, i.e., the Riemannian gradient descent, is lacking a clear and deep theoretical understanding on their convergence behaviors. This may largely be because the underlying optimization problem is geodesically non-convex:
min x∈Sn−1 f(x) = −1 2 x⊤Ax, (1)
where the objective function f(x) is constrained on the (n−1)-dimensional unit sphere manifold Sn−1 = {x ∈ Rn×1 : ∥x∥2 = 1} [1].
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
Previous analyses of the Riemannian gradient descent for Problem (1) achieve a rate of (global1) convergence at O( 1∆2 log 1 ϵ ) which depends quadratically on the relative gap ∆ between A’s two (distinct) largest eigenvalues [20, 19]2. The poor dependence of the rate on the gap characterizes a loose upper bound on the iteration complexity of the method, due to the artifacts of the used analysis techniques. Ding et al. [5] gave a tight convergence rate at O( 1∆ log 1 ϵ ). However, the analysis is limited to positive definite matrices. This means that the general tight convergence analysis for real symmetric matrices, in fact, has not been settled. Our first goal in this work is to give a general analysis for the tight rate of this type, i.e., dependent on the gap, which is termed as the (tight) gap-dependent analysis hereafter. More importantly, when the gap is significantly smaller than the target accuracy ϵ on the objective sub-optimality of the final solution rather than its point distance to the optimal solutions, the rate of this type is actually not tight any more, which calls for a worst-case rate independent of any gap. Our second goal is to further give the first worst-case or gap-independent analysis being able to achieve a rate of convergence at O( 1ϵ log 1 ϵ ). The two analyses naturally roll out a comprehensively tight rate of convergence at O( 1max{∆,ϵ} log 1 ϵ ). Our analyses are based on novel adaptations of the analyses for projection methods working for constant matrices [7, 10] to varying matrices during iterations. We also make a comparison between projection and search methods to help understand their subtle differences. Particularly, the gap-dependent analysis suggests a new step-size for search methods which gives rise to a novel learning rate for stochastic variance reduced PCA algorithms. Experiments verify our theoretical findings and demonstrate the effectiveness of the new learning rate. The main contributions are summarized as follows:
• We prove the first comprehensively tight convergence rate of the Riemannian gradient descent for Problem (1) at O( 1max{∆,ϵ} log 1 ϵ ).
• We propose an adaptive learning rate scheme for stochastic variance reduced PCA algorithms which can dramatically improve their convergence.
• We experimentally verify the established theory and demonstrate the effectiveness of the proposed adaptive learning rate.

heading: 2 Related Work
text: We focus on the gradient methods for the leading eigenvector computation. The classic Oja’s algorithm and Krasulina’s method, as stochastic counterparts of the projected gradient descent and Riemannian gradient descent with diminishing step-sizes, respectively, converge at a rate O( 1∆2ϵ ) [2] (Theorem 1.1 where E[ΨT ] = O( c 2 T ) and c = O( 1 ∆ ) using their notations except for the iteration number T ). The quadratic dependence on the gap was originally omitted there, probably because the gap was regarded as a constant. However, as argued in Musco et al. [10], this is insufficient. Interestingly, there are much fewer theoretical studies of the two methods in the deterministic setting. Absil et al. [1] introduced the Riemannian gradient descent for Problem (1) (Algorithm 2 there) with the Armijo line search for step-sizes, and implicitly stated a rate at O( 1∆ log 1 ϵ ) (Theorem 4.6.3-iii) which, however, is local (Theorem 4.6.3-ii). Wen et al. [17] proposed a special retraction for the Riemannian gradient descent based on the Cayley transform (Algorithm 1). It was empirically found comparable and more stable than the Lanczos algorithm, however, only with the guarantee of convergence to a critical point (Theorem 2). Since Problem (1) is geodesically non-convex, theories established for geodesically convex functions [23] can’t be applied unless locally. Note that we only consider global convergence in this work. Xu et al. [19] specifically analyzed the Riemannian gradient descent for Problem (1), and originally stated a rate O(min{ 1∆2 log 1ϵ , 1ϵ }) (Theorem 1). However, the constituent part O( 1ϵ ) resulted from the same mistake
3 as made in Balsubramani et al. [2] and mentioned at the beginning of this section. This means that their analyses are all gap-dependent and the better one of the rates is O( 1∆2 log 1 ϵ ). It is also the case for the work of Xu and Gao 2018 [20] (Theorem 4.1). Particularly, Ding et al. [5] improved the dependence of these rates on the gap from quadratic to linear at O( 1∆ log 1 ϵ ). As explained in the previous section, however, it is not general. In
1Only global convergence is considered throughout the paper. 2It needs be noted there that the rate O( 1
ϵ ) was mistakenly stated and actually should be O( 1 ∆2ϵ ) according
to their analyses. See the next section for explanations. 3Precisely, ψ(xT ) = O( c 2
T ), where c = O( 1 ∆p ) implied by a’s expression, using their notations.
addition, to the best of our knowledge, before our work, there has been no worst-case or gap-free analysis of the Riemannian gradient descent for Problem (1). Comparisons are summarized in Table 1.
Shamir [14] proposed a learning rate scheme for stochastic variance reduced PCA algorithm, namely VR-PCA, based on Oja’s algorithm. It can be used for VR-PCA’s variant as well that is based on Krasulina’s method [2, 16]. In this work, we put forward yet another learning rate scheme which we find much more effective than that suggested by [14] for these stochastic variance reduced PCA algorithms.

heading: 3 Notions and Notations
text: The only assumption on the given matrix A, made in our analyses, is real symmetry, i.e., A⊤ = A ∈ Rn×n. Let A’s orthonormal eigenvectors be v1,v2, · · · ,vn ∈ Sn−1, corresponding to eigenvalues λ1, λ2, · · · , λn ∈ R, i.e., Avi = λivi, in descending order. If λ1 > λn, then let ∆p (1 ≤ p ≤ n− 1) be the first nonzero eigenvalue gap, i.e., ∆p = λp − λp+1 > 0, such that λ1 = · · · = λp > λp+1. Particularly, the first nonzero relative eigenvalue gap is defined as ∆ = ∆pλ1−λn if λ1 > λn, otherwise as ∆ = +∞. The update equations of the two gradient methods for Problem (1) with constant step-sizes are as follows:
Projected gradient descent (PGD) : xt+1 = xt − η∇f(xt)
∥xt − η∇f(xt)∥2 , (2)
Riemannian gradient descent (RGD) : xt+1 = xt − η∇̃f(xt)
∥xt − η∇̃f(xt)∥2 , (3)
where ∇f(x)=−Ax and ∇̃f(x)=−(I−xx⊤)Ax (throughout the paper I stands for an appropriately sized identity matrix) represent the Euclidean and Riemannian gradients of f(x), respectively, η>0 constant step-size, and denominators for normalization, i.e., xt+1 ∈ Sn−1. Consider Algorithm that the Riemannian gradient descent for Problem (1) with constant step-size η starts from a random initial point x0 = x̂0∥x̂0∥2 ∈ S
n−1 (where x̂0 is entry-wise standard Gaussian), and then repeatedly runs Equation (3) for T iterations. We will focus on the analysis of this algorithm next. The results for the projected gradient descent are immediate by a much simpler path through the same analysis.

heading: 4 Gap-Dependent Analysis
text: Let Vp=[v1 · · · vp]. The solution space to Problem (1) then is {Vpα: α ∈ Sp−1} ≜ V , where S0 = {±1}. To study the convergence of iterates to V , it suffices to analyze the principal angle between xt and Vp, i.e., θ(xt,Vp) = cos−1(∥V⊤p xt∥2). Let V−p = [vp+1 · · · vn] and □† represent the pseudo-inverse of a matrix. Then tan θ(x0,Vp) = ∥(V⊤−px̂0)(V⊤p x̂0)†∥2. Suppose that T iterations are sufficient for sin2 θ(xt,Vp) = 1− ∥V⊤p xt∥22 to reach the target accuracy ϵ > 0, i.e., sin2 θ(xT ,Vp) < ϵ. The goal of our job here is to upper bound T with the gap information and show that the bound is tight, namely
Theorem 1 The Riemannian gradient descent for Problem (1) with step-size η = Θ( 1λ1−λn ) ≤ 1 λ1−λn converges in T = Θ( 1 ∆ log n ϵ ) iterations, i.e., sin
2 θ(xT ,Vp) < ϵ, w.h.p. (with high probability).
Proof First, it is trivial and T = 0 for the case that λ1 = λn since f(x) is a constant. Thus, it is always assumed that λ1 > λn in what follows. Let ht(x) = 1 + η(x − x⊤t Axt) and write ht(A) = ∑n i=1 ht(λi)viv ⊤ i = Vnht(Σn)V ⊤ n , where ht(Σn) = diag(ht(λ1), · · · , ht(λn)). We can write that
xT = hT−1(A)xT−1
∥hT−1(A)xT−1∥2 = · · · =
∏0 t=T−1 ht(A)x0
∥∏0t=T−1 ht(A)x0∥2 =
∏T−1 t=0 ht(A)x0
∥∏T−1t=0 ht(A)x0∥2 , (4)
where the second last equality is equivalently renormalized. Noting that the initial point x0 can be expressed as x0= ∑n i=1(v
⊤ i x0)vi, the last numerator above can be written as∏T−1
t=0 ht(A)x0 = ∑n i=1(v ⊤ i x0) ∏T−1 t=0 ht(A)vi = ∑n i=1(v ⊤ i x0) ∏T−1 t=0 ht(λi)vi.
Accordingly, it holds that
∥V⊤p ∏T−1 t=0 ht(A)x0∥22 = ∑p i=1(v ⊤ i x0) 2 ∏T−1 t=0 h 2 t (λi), ∥∏T−1t=0 ht(A)x0∥22 = ∑n i=1(v ⊤ i x0) 2 ∏T−1 t=0 h 2 t (λi).
By the fact that x⊤Ax ∈ [λn, λ1] for x ∈ Sn−1 (Corollary 4.7 in [15]), we have that ht(λi) ≥ 1 + η(λn − λ1) for all i and t. If 1 + η(λn − λ1) ≥ 0, i.e., η ≤ 1λ1−λn , it then holds that ht(λi) ≥ 0 for all i and t. When η ≤ 1λ1−λn , we can get that sin2 θ(xT ,Vp) = 1− ∥V⊤p xT ∥22
= 1− ∑p i=1(v ⊤ i x0) 2 ∏T−1
t=0 h 2 t (λi)∑n
i=1(v ⊤ i x0)
2 ∏T−1
t=0 h 2 t (λi)
=
∑n i=p+1(v ⊤ i x0) 2 ∏T−1
t=0 h 2 t (λi)∑n
i=1(v ⊤ i x0)
2 ∏T−1
t=0 h 2 t (λi)
(5)
≤ ∏T−1 t=0 h 2 t (λp+1) ∑n i=p+1(v ⊤ i x0) 2
∏T−1 t=0 h 2 t (λ1) ∑p i=1(v ⊤ i x0) 2 =
T−1∏
t=0
( ht(λp+1)
ht(λ1) )2 tan2 θ(x0,Vp)
≤ (1 + η(λp+1 − λn) 1 + η(λ1 − λn) )2T tan2 θ(x0,Vp)
≤ exp{− 2η∆pT 1 + η(λ1 − λn) } tan2 θ(x0,Vp) ≜ ϵ,
where the second last inequality have used the fact that 1+η(λp+1−x)1+η(λ1−x) is nonnegative and monotonically decreasing for x ∈ [λn, λ1]. Solving the last equality about ϵ for T yields that T = 1+η(λ1−λn)2η∆p log tan2 θ(x0,Vp) ϵ . Note that V ⊤ −px̂0 and (V ⊤ p x̂0)
† are entry-wise standard Gaussian since V−p and Vp are orthonormal. By standard Gaussian matrix concentration theory [12], w.h.p. ∥V⊤−px̂0∥22 ≤ c1n and ∥(V⊤p x̂0)†∥22 ≤ c2p for some fixed constants c1, c2. We then have w.h.p. that tan2 θ(x0,Vp) ≤ cnp for some fixed constant c and log tan2 θ(x0,Vp) = O(log(n)). Hence, T = O( 1∆ log n ϵ ) for η = Θ( 1 λ1−λn ).
To justify the tightness of this upper bound, we now show that the Riemannian gradient descent for Problem (1) requires Ω( 1∆ log n ϵ ) iterations to converge. To this end, we resume the analysis from Eq. (5) as follows:
sin2 θ(xT ,Vp) = 1− ∥V⊤p xT ∥22 = ∑n i=p+1(v ⊤ i x0) 2 ∏T−1 t=0 h 2 t (λi)∑n
i=1(v ⊤ i x0)
2 ∏T−1
t=0 h 2 t (λi)
≥ (v ⊤ p+1x0)
2 ∏T−1
t=0 h 2 t (λp+1)∑p+1
i=1 (v ⊤ i x0)
2 ∏T−1
t=0 h 2 t (λi)
= 1
1 + ∑p i=1(v ⊤ i x0) 2
(v⊤p+1x0) 2 ·
∏T−1 t=0 ( ht(λ1) ht(λp+1) )2
≥ 1 1 + ∑p i=1(v ⊤ i x0) 2
(v⊤p+1x0) 2 ·
∏T−1 t=0 ( 1 1−η∆p ) 2 =
(1− η∆p)2T tan2 θ̂ 1 + (1− η∆p)2T tan2 θ̂
≥ exp{−2η∆pT/(1− η∆p)} tan 2 θ̂
1 + exp{−2η∆pT/(1− η∆p)} tan2 θ̂ ,
where tan2 θ̂ = (v⊤p+1x0) 2
∑p i=1(v ⊤ i x0) 2 and we have used for the last inequality that (1 − η∆p)2T =
exp{2T log(1 − η∆p)} ≥ exp{− 2η∆pT1−η∆p } by log(1 + x) ≥ x 1+x for x > −1. Setting exp{− 2η∆pT1−η∆p } tan 2 θ̂ = 2ϵ yields that T = 1−η∆p2η∆p log tan2 θ̂ 2ϵ and sin
2 θ(xT ,Vp) ≥ 2ϵ1+2ϵ ≥ ϵ, assuming ϵ < 12 . Hence, for η = Θ( 1 λ1−λn ), we need T = Ω( 1 ∆ log n ϵ ) iterations to get sin2 θ(xT ,Vp) < ϵ, which completes the proof. □
Remark 1 In the gap-dependent analysis above, from the optimization point of view, we care more about the distance of the iterate xt to the optimal solution set V , i.e., sin2 θ(xt,Vp) = 1 2 minv∈V ∥xt − v∥22, with emphasis on the eigenvector computation. In fact, it is also common to consider the objective sub-optimality gap, i.e., (f(xt)−minx∈Sn−1 f(x)) = 12 (λ1 − x⊤t Axt), with emphasis on the eigenvalue computation. From the proof of Lemma 4.6 in Xu and Gao 2018 [20], we know that λ1 − x⊤TAxT ≤ (λ1 − λn) sin2 θ(xT ,Vp). By Theorem 1, with ϵ-rescaling, we then have that λ1 − x⊤TAxT < ϵ for T = O( 1∆ log nϵ ).

heading: 5 Gap-Independent Analysis
text: The issue of the gap-dependent rate in Theorem 1 is that when the gap ∆ is significantly smaller than the target accuracy ϵ, the convergence in terms of sin2 θ(xT ,Vp) is too slow, e.g., T = Θ( 1ϵ2 log n ϵ ) for ∆ = Θ(ϵ2). In this case, it is often acceptable to swap the true leading eigenvector for any ϵ-accurate solution xT ∈ Sn−1 in terms of the objective sub-optimality gap, i.e., λ1 − x⊤TAxT < ϵ. In this section, by extending the analysis of Musco et al. 2015 [10] for constant matrices (i.e., ht(x) = x), we show without the need of the eigenvalue gap information that even if the gap ∆ is significantly smaller than the target accuracy, we only need O( 1ϵ log n ϵ ) (far less than O( 1 ∆ log n ϵ ) for, e.g., ∆ = Θ(ϵ3/2), see Remark 1) iterations to achieve that λ1 − x⊤TAxT < ϵ, as stated in the following theorem on the worst-case performance.
Theorem 2 The Riemannian gradient descent for Problem (1) with step-size η = O(1) ≤ 1λ1−λn converges in T = O( 1ϵ log n ϵ ) iterations, i.e., λ1 − x⊤TAxT < ϵ.
Proof Sketch Due to space limit, we only give proof sketch here. The complete proof and all the missing proofs can be found in the supplementary material. We assume again that λ1 > λn, and η ≤ 1λ1−λn such that ht(λi) = 1 + η(λi − x ⊤ t Axt) ≥ 0 for all i and t. It suffices to show that λ1 − x⊤TAxT < 2η ϵ always holds no matter whether λ1 is significantly larger than λ2 in the sense that h0(λ1) ≥ (1 + δ/2)h0(λ2) for 0 < δ ≤ 2. We take T = ⌈ 2δ log n(1+tan2 θ0) ϵ ⌉+ 1, where θ0 = θ(x0,v1). Case 1 that h0(λ1)≥(1 + δ/2)h0(λ2). Consider the polynomial
pT (x) = √ (1 + δ/2)h0(λ2) ∏T−1 t=0 ht(x) (1+δ/2)ht(λ2) ≥ 0, x ∈ [λn, λ1]
and matrix form pT (A) = ∑n i=1 pT (λi)viv ⊤ i =VnpT (Σn)V ⊤ n , where pT (Σn)=diag(pT (λ1), · · · , pT (λn)). We have Fact 1. For x ∈ [λn, λ2], h0(λ1) ≥ (1 + δ/2)h0(x) implies ht(λ1) ≥ (1 + δ/2)ht(x) for all t, by
Lemma 3 If η ≤ 2λ1−λn then x ⊤ t+1Axt+1 ≥ x⊤t Axt. Hence, pT (λ1) ≥ √ h0(λ1). Also, pT (λi) ≤ √ h0(λ2)(1 + δ/2)
−T+ 12 , i = 2, · · · , n. Eq. (4) can be rewritten as xT =
pT (A)x0 ∥pT (A)x0∥2 , where the rank-1 approximation inequality for pT (A) holds:
∥pT (A)− xTx⊤T pT (A)∥2F ≤ (1 + tan2 θ0)∥pT (A)− pT (λ1)v1v⊤1 ∥2F = (1 + tan2 θ0) ∑n i=2 p 2 T (λi)
≤ (1 + tan2 θ0)(n− 1)h0(λ2)(1 + δ/2)−2T+1 < h0(λ2)ϵ.
At the same time,
∥pT (A)− xTx⊤T pT (A)∥2F = ∥pT (A)∥2F − ∥xTx⊤T pT (A)∥2F = ∑n
i=1(1− (x⊤T vi)2)p2T (λi) ≥ (1− (x⊤T v1)2)p2T (λ1) ≥ (1− (x⊤T v1)2)h0(λ1).
Thus, λ1 − x⊤TAxT < 2η ϵ holds by noting that h0(λ1)− x⊤T h0(A)xT = h0(λ1)− ∑n i=1(x ⊤ T vi) 2h0(λi) ≤ (1− (x⊤T v1)2)h0(λ1) < ϵh0(λ2).
Case 2 that h0(λ1) < (1+ δ/2)h0(λ2). Consider the polynomial qT (x) = √ h0(λ1) ∏T−1
t=0 ht(x) ht(λ1) .
and write xT = qT (A)x0 ∥qT (A)x0∥2 . Define the index set α = {i : 1 1+δ/2h0(λ1) ≤ h0(λi) < h0(λ1)} and rewrite
xT =
∑ i∈α qT (λi)viv ⊤ i x0
∥qT (A)x0∥2 +
∑ i/∈α qT (λi)viv ⊤ i x0
∥qT (A)x0∥2 ≜ x̃(α)T + x̃ (−α) T .
We then have ∥h 1 2 0 (A)x (α) T ∥22 ≥ 11+δ/2h0(λ1) ≥ (1 − ϵ)h0(λ1) for δ = 2ϵ, where the first inequality is by the definition of the index set. For x(−α)T ≜ x̃ (−α) T
∥x̃(−α)T ∥2 = qT (A−α)x0∥qT (A−α)x0∥2 and
A−α ≜ ∑ i/∈α λiviv ⊤ i , a similar argument to Case 1 can be applied to get that h0(λ1) − (x (−α) T ) ⊤h0(A−α)x (−α) T < h0(λ1)ϵ. Then ∥h 1 2 0 (A)x (−α) T ∥22 = (x (−α) T ) ⊤h0(A−α)x (−α) T > (1− ϵ)h0(λ1). Thus, λ1 − x⊤TAxT < 2η ϵ holds by noting that
x⊤T h0(A)xT = ∥h 1 2 0 (A)x̃ (α) T ∥22 + ∥h 1 2 0 (A)x̃ (−α) T ∥22
= ∥x̃(α)T ∥22 ∥h 1 2 0 (A)x (α) T ∥22 + ∥x̃ (−α) T ∥22 ∥h 1 2 0 (A)x (−α) T ∥22 > ( ∥x̃(α)T ∥22 + ∥x̃ (−α) T ∥22 )(1− ϵ)h0(λ1) = (1− ϵ)h0(λ1).
□
Remark 2 When the gap is significantly larger than the target accuracy ϵ, e.g., ∆ = O(1), and sin2 θ(xT ,Vp) ≥ ϵ(λ1−λn)∆ for T = Θ( 1 ∆ log n(λ1−λn)∆ ϵ ) (see the proof of Theorem 1 on the tightness), we have that λ1 − x⊤TAxT ≥ (λ1 − λn)∆ sin2 θ(xT ,Vp) ≥ ϵ by Lemma 2 in Xu et al. [19]. Thus, we need T = Ω( 1∆ log n ϵ ) iterations to get λ1 − x⊤TAxT < ϵ. Otherwise, it needs O( 1ϵ log n ϵ ) iterations to make it. Thus, a comprehensively tight rate can be stated as O( 1max{∆,ϵ} log n ϵ ).
Remark 3 The analysis of Case 1 in the proof of Theorem 2 can give us a gap-dependent rate:
h0(λ1) = 1 + η(λ1 − x⊤0 Ax0) = (1 + η∆1 1 + η(λ2 − x⊤0 Ax0) )(1 + η(λ2 − x⊤0 Ax0))
≥ (1 + η∆1 1 + η(λ2 − λn) )h0(λ2) ≥ (1 + η∆1 2 )h0(λ2).
Taking δ = η∆1 in T = ⌈ 2δ log n(1+tan2 θ0) ϵ ⌉+ 1 gives us T = O( 1η∆1 log n ϵ ). It is straightforward to extend to the setting λ1 = · · · = λp > λp+1 of the previous section to get T = O( 1∆ log nϵ ) with ∆p and η = Θ( 1λ1−λn ). But it is worth mentioning this is about the convergence in terms of λ1 − x⊤TAxT < ϵ, rather than sin2 θ(xT ,Vp) < ϵ.
Remark 4 Compared to ht(A) ≡ A which leads to simple monomials in [7, 10], ht(A) = I+ η(A− x⊤t Axt · I) varies with iterations and requires complicated polynomials in our analysis. In addition, from Remark 1-2, it holds for ∥x∥2 = 1 that 12∆p sin2 θ(x,Vp) ≤ f(x) − f(v1) ≤ 1 2 (λ1 − λn) sin2 θ(x,Vp). Note that either these inequalities or Wedin’s theorem can’t be used to translate results between Theorem 1 and Theorem 2, because the translated results by these inequalities may possibly be worse and Wedin’s theorem (which upper bounds perturbation of eigenspace in
subspace distance by those of both the given matrix and corresponding gap) is inapplicable here. For example, by the second inequality and Theorem 1, we can get λ1−xTAxT < ϵ for T = O( 1∆ log nϵ ). But when ∆ is much less than ϵ, T will clearly be worse than the results in Theorem 2. On the other hand, by the first inequality and Theorem 2, we can get sin2 θ(xT ,Vp) < ϵ for T = O( 1ϵ∆p log n ϵ∆p
), which is worse than the results in Theorem 1.

heading: 6 Projection Versus Search
text: To be specific, in this section, we refer to the power iteration (PI) and Riemannian gradient descent (RGD) as the projection and search methods for Problem (1), respectively, and the projected gradient descent (PGD) belongs to both classes.
The above analyses for the search method apply to the PGD by simply setting ht(x) = 1 + ηx. In this case, η represents the step-size and needs satisfy 0 < η ≤ 1ρ to make sure ht(λi) ≥ 0 for all i, where ρ = maxi |λi| is A’s spectral radius. Assuming that ∆1 > 0, it will converge to v1. On the other hand, the analyses [7, 10] of the PI apply to the PGD as well. Then η can be any nonzero real value and the PGD converges to vσ(1) assuming that the largest two values of {|ht(λj)| : j = 1, · · · , n} are distinct, where σ maps i to j such that |ht(λj)| is the i-th largest value of {|ht(λj)| : j = 1, · · · , n}. In addition, the PI can be seen as the PGD with constant step-size η = ±∞, and can also be interpreted as the RGD with varying step-sizes ηt = 1
x⊤t Axt when A is
positive definite.
Table 2 summarizes the comparisons between these methods for Problem (1), where ϵ represents the target accuracy for λ1 − x⊤TAxT , and (x⊤TAxT ,xT ) describes which eigenpair it approximates. Table 3 shows that the PGD even with varying step-sizes 0 < ηt ≤ 1ρ always gets a larger contraction factor ht(λ2)ht(λ1) (see the equation right below Eq. (5)) than the PI, but the RGD can achieve a smaller one when step-sizes satisfy that ηt > 1x⊤t Axt . In practice, it may be worth a try to set ηt = τ x⊤t Axt with τ > 1. However, the Riemannian gradient descent update requires more computations than the power iteration in each step, possibly leading to an even worse wall-clock running time. This can be easily seen from the equivalence between the two methods with ηt = 1x⊤t Axt .
In the next section, we will investigate the effect of this step-size scheme ηt = τx⊤t Axt (τ > 1) on both projected and Riemannian gradient descent for eigenvector computation as well as the stochastic variance reduced algorithms for PCA (VR-PCA) based on Oja’s algorithm or Krasulina’s method.

heading: 7 Experiments
text: The purpose of the experimental study for corroborating our findings in above sections is twofold. One is to verify the rate O((1/max{∆, ϵ}) log(n/ϵ)) of the Riemannian gradient descent on real symmetric matrices in the form of A ∈ Rn×n. The other is to experiment on the step-size scheme for both projected and Riemannian gradient descent methods on PCA data in the form of C ∈ Rd×n such that A = 1nCC ⊤. Experiments were done on a laptop (dual-core 2.30GHZ CPU and 8GB RAM).
The constituent part O((1/∆) log(n/ϵ)) of the rate has been experimentally verified in Ding et al. [5]. Hence, only the worst-case rate O((1/ϵ) log(n/ϵ)) needs verification, for which we generate synthetic data A=VΣV⊤ as follows. Given n=1000, s=100,∆=0.01 and entry-wise i.i.d. standard normal U ∈ Rn×n,a ∈ R(n−s−1)×1, we consider orthogonal V = U(U⊤U)−1/2 and diagonal Σ = diag(2I−∆diag(0, 1, · · · , s), 1ndiag(|a1|, · · · , |an−s−1|)). In this case, the convergence for each method considered in this work should behave differently in terms of different measures, i.e., O((1/ϵ) log(n/ϵ)) for (λ1 − x⊤TAxT ) versus O((1/∆) log(n/ϵ)) for sin2 θ(xT ,v1).
We implemented the PGD with η=1ρ , ηt = 0.6 x⊤t Axt , 1 x⊤t Axt , 1.6 x⊤t Axt , and RGD with step-size schemes η = 1λ1−λn , ηt = 0.6 x⊤t Axt , 1 x⊤t Axt , 1.6 x⊤t Axt
, in MATLAB. Each is compared to the PI. All the methods start from the same random initial point x0 and run for T = 100 iterations. Figure 1 reports their performance across iterations (or the average performance over time on 10 runs) in two measures, where the y-axis of the figures in the first and third columns is in log scale of (λ1 −x⊤t Axt)/λ1. It is clear for each method that (λ1 − x⊤t Axt) converges sub-linearly while sin2 θ(xt,v1) is decreasing even more slowly, matching their respective theoretical rates. Consider the case that ∆ = 0.01 and accuracy ϵ = 0.1 here. First, from the figures on the left in Figure 1, it is clear that the (relative) target in objective converges sub-linearly (note that y-axis is log scale here), roughly matching O((1/ϵ) log(n/ϵ)) iterations. From the figures on the right, the target sin2 θ(xt,v1) converges about one order of magnitude slower than the target in objective, roughly matching O((1/ϵ)2 log(n/ϵ)) iterations which is basically consistent with O((1/∆) log(n/ϵ)) for ϵ = 0.1.
We also experiment on two real datasets4 Schenk that was used in [19, 5] and GHS_indef. Schenk is a 10728× 10728 real symmetric matrix with 85000 nonzero entries. The performance is depicted in Figure 2, where all the methods converge linearly due to a large gap ∆ and their differences are sharper. Particularly, we add the plot of the average performance over time on 10 runs. Although the RGD with ηt = 1x⊤t Axt is theoretically the same as the PI (thus performing exactly the same over iterations in Figures 1-2), it needs more computational costs reflected by the wall-clock time in Figure 2. Nonetheless, Figures 1-2 showcase the great potential of the RGD with aggressive step-size schemes to outperform the PI, e.g., ηt = 1.6x⊤t Axt here, in practice. In the meanwhile, the performance of the PGD even with adaptive step-sizes is worse than that of the PI. We tested significantly larger
4https://sparse.tamu.edu/
constant τ in τ x⊤t Axt . PGD’s performance is still worse than that of the PI. Similar results appear on GHS_indef which is a 76638× 76638 real symmetric matrix with 859520 nonzero entries.
Algorithm 1 VR-PCA (Oja) 1: Input: C = (c1, · · · , cn) ∈ Rd×n, x̃0,
τ , b = |αt|, m = ⌈nb ⌉, S, where [n] ={1, 2, · · · , n} 2: Output: xS 3: for s = 0, 1, · · · , S − 1 do 4: gs = 1 nCC
⊤x̃s, x0 = x̃s, ηs = nτ∥C⊤x̃s∥22 5: for t = 0, 1, · · · ,m− 1 do 6: Pick αt ⊂ [n] uniformly at random 7: gt = 1 b ∑ i∈αt cic ⊤ i (xt − x̃s) + gs 8: xt+1 = xt+ηsgt ∥xt+ηsgt∥2 9: end for
10: x̃s = xm 11: end for
Algorithm 2 VR-PCA (Krasulina) 1: Input: C = (c1, · · · , cn) ∈ Rd×n, x̃0, τ ,
b = |αt|, m = ⌈nb ⌉, S 2: Output: xS 3: for s = 0, 1, · · · , S − 1 do 4: g̃s = 1 n (CC
⊤x̃s−∥C⊤x̃s∥22x̃s), x0 = x̃s, ηs =
nτ ∥C⊤x̃s∥22
5: for t = 0, 1, · · · ,m− 1 do 6: Pick αt ⊂ [n] uniformly at random 7: g̃t = 1 b ∑ i∈αt(cic ⊤ i (xt − x̃s) −
((c⊤i xt) 2xt − (c⊤i x̃s)2x̃s) + g̃s
8: xt+1 = xt+ηsg̃ ∥xt+ηsg̃t∥2 9: end for
10: x̃s = xm 11: end for
We now check the performance of two versions of the VR-PCA based on Oja’s algorithm and Krasulina’s method, corresponding to Algorithm 1 in the work of Shamir 2015 [14] (Algorithm 1 above) and variance reduced version (see the supplementary material) of Algorithm 1 in the work of Tang 2019 [16] (Algorithm 2 above), respectively, where the constant learning rate is replaced with the proposed adaptive learning rate scheme. The common PCA datasets are used and summarized in Table 4. We use b = 100. Note that we only update the learning rate at the epoch level and keep it unchanged within each epoch, similar to the case of the computation of the full gradient. Figure 3 reports the average performance of the VR-PCA with our adaptive learning rate on three real datasets, following the previous experimental setting on Schenk except for 5 runs used for average. It is easy to see that the proposed learning rate scheme (i.e., ηs = τn∥Cx̃s∥22 ) in the PCA setting consistently
works far better than that proposed in Shamir 2015 [14] (i.e., η = √ n
∥C∥2F ) across the datasets, for both
versions of the VR-PCA. Interestingly, an opposite case (see Figure 2 for the deterministic setting) happens that the performance of the Oja-based VR-PCA performs significantly better than that based on Krasulina’s method. This shows that the proposed learning rate scheme is very promising for the practitioners of the VR-PCA.

heading: 8 Conclusions
text: We presented two analyses of the Riemannian gradient descent for eigenvector computation. Specifically, the gap-dependent one is the first general analysis that achieves a tight rate O( 1∆ log n ϵ ) in terms of the point distance to the optimal solutions. However, if we are only concerned with the objective distance to the optimal solutions, the rate of this type is not worst-case optimal. We succeeded in establishing a worst-case rate O( 1ϵ log n ϵ ) independent of any gap, by adapting an existing analysis for projection methods to our context of which a prominent characteristic lies in varying matrices ht(A) during iterations, and ended up with a comprehensively tight rate O( 1max{∆,ϵ} log n ϵ ). In addition, we discussed the equivalence and subtle differences between projection and search methods that give rise to new insights on the step-size schemes, and conducted experimental studies to verify the established theories. Particularly, the proposed adaptive learning rate scheme for the VR-PCA works far better than those proposed previously. One limitation of our analyses lies in the great difficulties of their extensions to the block setting where X ∈ Rn×k for k ≥ 1 (this work corresponds to k = 1), because X⊤t AXt is no longer a scalar but a matrix and the key matrix function ht(A) then is not well-defined for the analysis. We consider overcoming the underlying difficulties in this setting as an important future work.

Review summaries:
In this paper, the authors considered Riemanninan gradient descent for solving the top eigenvector/value problem for real-symmetric matrices. Rates in two regimes are given, depending on the gap of the top eigenvalues compared to the desired accuracy. The analysis is based on polynomial recursion derived from the algorithm updates. Additionally, for positive-definite matrices, a stepsize choice is proposed based on the upper bound on the contraction factor as suggested by the recursion. Numerical experiments are conducted for the power iteration, projected GD and Riemanninan GD.

The paper studies Riemannian gradient descent for largest eigen-vector computation real symmetric matrices. For matrices with eigen-gap $\Delta$ and (global) convergence error $\epsilon$, the paper gives a gap-dependent convergence rate of $\frac{1}{\Delta} \log(1/\epsilon)$ and gap-independent analysis which gives a convergence rate of $\frac{1}{\epsilon} \log(1/\epsilon)$. This improves upon previous work that achieves a $\frac{1}{\Delta^2} \log(1/\epsilon)$ rate and generalizes related work that that achieves the same rate for projected gradient descent, power iteration and Riemannian gradient descent restricted to PSD matrices.  The paper highlights an omission in previous work that mistakenly stated a rate of $\frac{1}{\epsilon} \log(1/\epsilon)$ which is in fact $\frac{1}{\Delta^2\epsilon} \log(1/\epsilon)$.   Experiments comparing the power iteration, projected gradient descent and RGD are done for proposed adaptively chosen step-sizes. 

This paper considers Riemannian gradient method for PCA. Contributions of the paper are following: 1) proved convergence rate of $O(\frac{1}{\max{\Delta, \epsilon}}\log\frac{1}{\varepsilon})$ for any real symmetric matrix; 2) showed without the need of the eigenvalue gap information that even if the gap $\Delta$ is significantly smaller than the target accuracy, it is possible to achieve a speed of convergence $O(\frac{1}{\epsilon}\log\frac{1}{\varepsilon})$. Numerical results verify obtained rate of convergence.

In this work authors presented two analyses of the Riemannian gradient descent for eigenvector computation.  First authors provided a general tight analysis that holds for any real symmetric matrix. In addition, authors further give the first worst-case analysis that achieves a rate of convergence.  For real-world application, this gap-dependent analysis suggests a new promising learning rate for stochastic variance reduced PCA algorithms. Authors provided experiments to confirm the findings.  Overall this is a solid paper with both theoretical analysis and empirical studies.