Noisy Accelerated Power Method for Eigenproblems with Applications

TL;DR

This paper presents a scalable, noise-tolerant algorithm for efficiently computing dominant eigenvectors of symmetric matrix pairs, combining approximation theory and convex optimization to outperform traditional methods in speed and accuracy.

Contribution

The paper introduces a novel noisy accelerated power method that decomposes eigenproblem solutions into approximate subproblems, achieving faster convergence with theoretical guarantees.

Findings

Linear running time in input size

Significant acceleration in canonical correlation analysis

Empirical results show improved accuracy and speed

Abstract

This paper introduces an efficient algorithm for finding the dominant generalized eigenvectors of a pair of symmetric matrices. Combining tools from approximation theory and convex optimization, we develop a simple scalable algorithm with strong theoretical performance guarantees. More precisely, the algorithm retains the simplicity of the well-known power method but enjoys the asymptotic iteration complexity of the powerful Lanczos method. Unlike these classic techniques, our algorithm is designed to decompose the overall problem into a series of subproblems that only need to be solved approximately. The combination of good initializations, fast iterative solvers, and appropriate error control in solving the subproblems lead to a linear running time in the input sizes compared to the superlinear time for the traditional methods. The improved running time immediately offers acceleration for several applications. As an example, we demonstrate how the proposed algorithm can be used to accelerate canonical correlation analysis, which is a fundamental statistical tool for learning of a low-dimensional representation of high-dimensional objects. Numerical experiments on real-world data sets confirm that our approach yields significant improvements over the current state-of-the-art.

Figures

Tables

Table 1. Table 1: Summary of different algorithms solving CCA

Problem	Algorithm	LS solver	Time complexity \bigstrut
CCA-1	Appgrad[17]	GD	$𝒪\left(\frac{nd\kappa }{\Delta }\log \frac{1}{ϵ}\right)$*
	CCALin[9]	SVRG	$𝒪\left(\frac{d\left(n+\sqrt{n\tilde{\kappa }}\right)}{\Delta }\log \frac{1}{ϵ}\right)$
	ALS-VR[31]	SVRG	$𝒪\left(\frac{d\left(n+\tilde{\kappa }\right)}{{\Delta }^2}{\log }^2\frac{1}{ϵ}\right)$
	SI[31]	ASVRG	$𝒪\left(\frac{d{n}^{3/4}\sqrt{\tilde{\kappa }}}{\sqrt{\Delta \left|{\lambda }_1\right|}}{\log }^2\frac{1}{ϵ}\right)$
	LazyCCA[2]	ASVRG	$𝒪\left(\frac{d\left(n+\sqrt{n\tilde{\kappa }}\right)}{\sqrt{\Delta }}\log \frac{1}{ϵ}\right)$
	NAPI	ASVRG	$𝒪\left(\frac{d\left(n+\sqrt{n\tilde{\kappa }}\right)}{\sqrt{\Delta }}\log \frac{1}{ϵ}\right)$
CCA-$k$	Appgrad[17]	GD	$𝒪\left(\frac{ndk\kappa }{{\Delta }_k}\log \frac{1}{ϵ}\right)$*
	CCALin[9]	SVRG	$𝒪\left(\frac{dk\left(n+\sqrt{n\tilde{\kappa }}\right)}{{\Delta }_k}\log \frac{1}{ϵ}\right)$
	LazyCCA[2]	ASVRG	$𝒪\left(\frac{dk\left(n+\sqrt{n\tilde{\kappa }}\right)}{\sqrt{{\Delta }_k}}\log \frac{1}{ϵ}\right)$
	NAPI	ASVRG	$𝒪\left(\frac{dk\left(n+\sqrt{n\tilde{\kappa }}\right)}{\sqrt{{\Delta }_k}}\log \frac{1}{ϵ}\right)$

Table 2. Table 2: Brief Summary of Data sets

Data set	Description	$d_1$	$d_2$	$n$ \bigstrut
Mediamill	Image and its labels	101	120	30,993 \bigstrut
XRMB	Acoustic and articulation measurements	273	112	1,429,236
MNIST	Left and right halves of images	392	392	60,000

Equations

$$\mathbf{A}\mathbf{w}=\lambda \mathbf{B}\mathbf{w}$$

$$\begin{array}{rlrl}& \underset{\mathbf{x},\mathbf{y}}{\mathrm{maximize}}& & {\mathbf{x}}^{\top }{\mathbf{\Sigma }}_{12}\mathbf{y}\\ & \text{subject to}& & {\mathbf{x}}^{\top }{\mathbf{\Sigma }}_{11}\mathbf{x}={\mathbf{y}}^{\top }{\mathbf{\Sigma }}_{22}\mathbf{y}=1.\end{array}$$

$$\underset{\mathbf{x}}{\max }\frac{{\mathbf{x}}^{\top }{\mathbf{\Sigma }}_{11}\mathbf{x}}{{\mathbf{x}}^{\top }{\mathbf{\Sigma }}_{22}\mathbf{x}}\text{and}\underset{\mathbf{y}}{\max }\frac{{\mathbf{y}}^{\top }{\mathbf{\Sigma }}_{22}\mathbf{y}}{{\mathbf{y}}^{\top }{\mathbf{\Sigma }}_{11}\mathbf{y}}.$$

$$\cos {\theta }^{(i)}\left(\mathcal{X},\mathcal{Y}\right)=\underset{\mathbf{x}\in \mathcal{X}}{\max }\underset{\mathbf{y}\in \mathcal{Y}}{\max }{\mathbf{x}}^{\top }\mathbf{B}\mathbf{y}=\left|{\mathbf{x}}_i^{\top }\mathbf{B}{\mathbf{y}}_i\right|,$$

$$\cos \theta \left(\mathbf{X},\mathbf{Y}\right)={\sigma }_k\left({\mathbf{X}}^{\top }\mathbf{BY}\right)\text{and}\sin \theta \left(\mathbf{X},\mathbf{Y}\right)=\Vert {\mathbf{X}}_{\perp }^{\top }\mathbf{BY}\Vert .$$

$$\tan \theta \left(\mathbf{X},\mathbf{Y}\right)=\Vert {\mathbf{X}}_{\perp }^{\top }\mathbf{BY}{\left({\mathbf{X}}^{\top }\mathbf{BY}\right)}^{-1}\Vert .$$

$$\Vert \mathbf{X}{\mathbf{X}}^{\top }\mathbf{B}-\mathbf{Y}{\mathbf{Y}}^{\top }\mathbf{B}\Vert$$

$$\underset{{\mathbf{Q}}^{\top }\mathbf{Q}={\mathbf{I}}_k}{\min }\Vert {\mathbf{B}}^{1/2}\mathbf{X}-{\mathbf{B}}^{1/2}\mathbf{YQ}\Vert$$

$$p_s\left(\mathbf{A}\right)$$

$$\underset{x\in \left[-1,1\right]}{\sup }\left|p_{s,d}\left(x\right)-x^s\right|\le 2e^{-d^2/2s}.$$

$$\underset{x\in \left[-1,1\right]}{\sup }\left|p_{s,d}\left(x\right)-x^s\right|\le \delta .$$

$$f\left(x\right):=\frac{1}{n}\sum _{i=1}^n{f}_i\left(x\right),$$

$$f\left(\mathbf{w}\right):=\frac{1}{2}{\mathbf{w}}^{\top }\mathbf{B}\mathbf{w}-{\mathbf{w}}^{\top }\mathbf{b},$$

$$f_i(\mathbf{w})$$

$$\begin{array}{rlrl}{\mathbf{w}}_i\in & \underset{\mathbf{w}}{\mathrm{argmax}}& & |{\mathbf{w}}^{\top }\mathbf{A}\mathbf{w}|\\ & \text{mbox}subjectto& & {\mathbf{w}}^{\top }\mathbf{B}\mathbf{w}=1\\ & & & {\mathbf{w}}^{\top }\mathbf{B}{\mathbf{w}}_j=0\forall j\in \{1,2,\dots ,i-1\}.\end{array}$$

$$T_0\left(z\right)$$

$$T_t\left(z\right)$$

$${\mathbf{y}}_{t+1}=2{\mathbf{B}}^{-1}\mathbf{A}{\mathbf{y}}_t-{\mathbf{y}}_{t-1}.$$

$${\mathbf{x}}_{t+1}={\mathbf{B}}^{-1}\mathbf{A}{\mathbf{x}}_t-\beta {\mathbf{x}}_{t-1},$$

$${\tilde{\mathbf{w}}}_{t+1}:={\mathbf{B}}^{-1}\mathbf{A}{\mathbf{w}}_t-\beta {\mathbf{w}}_{t-1},{\gamma }_{t+1}={\Vert {\tilde{\mathbf{w}}}_{t+1}\Vert }_{\mathbf{B}}$$

$${\mathbf{w}}_t:={\mathbf{w}}_t{\gamma }_{t+1}^{-1},{\mathbf{w}}_{t+1}:={\tilde{\mathbf{w}}}_{t+1}{\gamma }_{t+1}^{-1},$$

$${\mathbf{w}}_T=\frac{C_T\left({\mathbf{B}}^{-1}\mathbf{A}\right){\mathbf{w}}_0}{{\Vert C_T\left({\mathbf{B}}^{-1}\mathbf{A}\right){\mathbf{w}}_0\Vert }_{\mathbf{B}}},$$

$${\mathbf{w}}_{t+1}=\left({\mathbf{B}}^{-1}\mathbf{A}{\mathbf{w}}_t-\beta {\mathbf{w}}_{t-1}{\gamma }_t^{-1}\right){\gamma }_{t+1}^{-1}.$$

$${\mathbf{w}}_{t+1}$$

$$=\left({\mathbf{B}}^{-1}\mathbf{A}{\mathbf{x}}_t{\alpha }_t^{-1}-\beta {\mathbf{x}}_{t-1}{\alpha }_t^{-1}\right){\gamma }_{t+1}^{-1}$$

$$=\left({\mathbf{B}}^{-1}\mathbf{A}{\mathbf{x}}_t-\beta {\mathbf{x}}_{t-1}\right){\alpha }_t^{-1}{\gamma }_{t+1}^{-1}$$

$$={\mathbf{x}}_{t+1}{\alpha }_{t+1}^{-1},$$

$${\sin }^2\theta \left({\mathbf{w}}_T,{\mathbf{u}}_1\right)\le 4{\left(\frac{2\sqrt{\beta }}{\left|{\lambda }_1\right|+\sqrt{{\lambda }_1^2-4\beta }}\right)}^{2T}{\tan }^2\theta \left({\mathbf{w}}_0,{\mathbf{u}}_1\right).$$

$$\underset{\mathbf{w}\in \text{mathbbmss}{R}^d}{\min }f\left(\mathbf{w}\right):=\frac{1}{2}{\mathbf{w}}^{\top }\mathbf{B}\mathbf{w}-{\mathbf{w}}^{\top }\mathbf{A}{\mathbf{w}}_t,$$

$$f\left(\alpha {\mathbf{w}}_t\right)=\frac{1}{2}{\alpha }^2{\mathbf{w}}_t^{\top }\mathbf{B}{\mathbf{w}}_t-\alpha {\mathbf{w}}_t^{\top }\mathbf{A}{\mathbf{w}}_t.$$

$${\mathbf{\vartheta }}_{t+1}:=\left[\begin{array}{c}{\mathbf{w}}_{t+1}\\ {\mathbf{w}}_t\end{array}\right]$$

$$=\frac{1}{{\gamma }_{t+1}}\left(\mathbf{C}{\mathbf{\vartheta }}_t+{\mathbf{\xi }}_t\right),$$

$$\mathbf{C}=\left[\begin{array}{cc}{\mathbf{B}}^{-1}\mathbf{A}& -\beta \mathbf{I}\\ \mathbf{I}& 0\end{array}\right]\text{and}{\mathbf{\xi }}_t=\left[\begin{array}{c}{\mathbf{e}}_t\\ 0\end{array}\right].$$

$${\mu }_i^2-{\lambda }_i{\mu }_i+\beta =0.$$

Full text

1 Introduction

Computing the leading eigenvectors of a matrix is one of the most fundamental problems in scientific computing. The problem has been studied extensively in a classical setting and is used in countless applications. For generalized eigenvector problems, standard iterative methods such as the power method and the Lanczos method require multiple matrix-vector products of the form $\mathbf{B}^{-1}\bm{x}$ at each iteration [22, Section 9.2.6], which becomes prohibitively expensive in high dimensions. The poor scalability of the standard techniques make them unfit for many applications in machine learning and data science. In these applictions, the ability to deal efficiently with large-scale data is the main concern, while machine precision accuracy is no longer essential since the problem data is often uncertain and subject to noise. Novel algorithms for scalable and approximate eigenvalue computations are therefore urgently needed.

In this paper, we consider the problem of finding the solution $(\bm{w},\lambda)$ to the following equation

[TABLE]

where $\mathbf{A}$ and $\mathbf{B}$ are given real, symmetric matrices and $\mathbf{B}$ is positive definite. Such eigenproblems appear in a wide range of machine learning tasks, such as finding a set of directions in the data-embedding space which contains the maximum amount of variance (principal components analysis), finding a hyperplane that separates two classes of data (Fisher discriminant), and finding correlations between two different representations of the same data (canonical correlation analysis), among many others.

Let $\mathbf{X},\mathbf{Y}\in\mathbbmss{R}^{n\times d}$ be two (potentially large) sets of data with empirical covariance matrices $\mathbf{\Sigma}_{11}=\frac{1}{n}\mathbf{X}^{\top}\mathbf{X}$, $\mathbf{\Sigma}_{22}=\frac{1}{n}\mathbf{Y}^{\top}\mathbf{Y}$, and $\mathbf{\Sigma}_{12}=\frac{1}{n}\mathbf{X}^{\top}\mathbf{Y}$. The canonical correlation analysis (CCA) problem is concerned with finding the features $\mathbf{x}$ and $\mathbf{y}$ that best encapsulate the similarity of $\mathbf{X}$ and $\mathbf{Y}$. In particular, CCA maximizes the empirical cross-correlation between $\mathbf{X}$ and $\mathbf{Y}$ by solving the following problem

[TABLE]

Generalized eigenvalue (GEV) problems, on the other hand, compute features that maximize discrepancies between the data sets as measured by the the quantities

[TABLE]

Both problems can be reduced to performing principal component analysis (PCA) on complicated matrices, e.g., $\mathbf{\Sigma}_{22}^{-1/2}\mathbf{\Sigma}_{12}^{\top}\mathbf{\Sigma}_{11}^{-1}\mathbf{\Sigma}_{12}\mathbf{\Sigma}_{22}^{-1/2}$ for CCA and $\mathbf{\Sigma}_{22}^{-1/2}\mathbf{\Sigma}_{11}\mathbf{\Sigma}_{22}^{-1/2}$ for GEV. However, this approach is not scalable, as the formation of $\mathbf{\Sigma}_{11}^{-1/2}$ and $\mathbf{\Sigma}_{22}^{-1/2}$ becomes prohibitively expensive for large data sets.

1.1 Related Work

While there have been significant recent progress in understanding PCA, efficient algorithms and theoretical guarantees remain limited for generalized eigenvalue problems and CCA. Below, we summarize recent work in optimization and approximation theory that have influenced the developments in this paper.

1.1.1 Stochastic optimization for PCA

Over the past few years, the challenges of dealing with huge data sets have inspired the development of novel optimization algorithms for empirical risk minimization (see, e.g., [3] and references therein). Leveraging ideas from stochastic optimization, several scalable algorithms have also been proposed for different eigenvalue problems [8, 26, 27, 14, 21]. These algorithms try to combine the low iteration cost of the stochastic methods and the high accuracy of the standard deterministic techniques. Notably, borrowing the idea of variance reduction [15], the author in [26] studied a variant of Oja’s algorithm [19], giving the first linearly convergent algorithm for stochastic PCA. The authors of [8] improved the result further by using the well-known shift-and-invert technique from numerical linear algebra [11, 22]. Inspired by the momentum method as known as the Heavy ball method in the optimization literature [20], an accelerated stochastic PCA was proposed in [21].

1.1.2 Approximation theory and acceleration

Approximation theory is ubiquitous in machine learning, optimization, and numerical linear algebra. This is because many problems in these areas can be reduced to finding a faster way to compute primitives such as $x^{s}$, $x^{-1}$, or $e^{x}$ [23]. For example, for the solution to linear systems of equations of the form $Ax=b$ with $A\in\mathbbmss{R}^{n\times n}$, Chebyshev iteration improves the iteration complexity from $\mathcal{O}\left(n\kappa\log\frac{1}{\epsilon}\right)$ of the Richardson iteration to $\mathcal{O}\left(n\sqrt{\kappa}\log\frac{1}{\epsilon}\right)$, where $\kappa$ is the condition number of $A$ [11, Chapter 11]. This improved rate is similar to the speed-ups by Nesterov’s fast gradient method over the classical gradient descent [18], which was recently discussed in terms of polynomial approximations in [24]. For eigenproblems, the Lanczos method can be considered as an accelerated version of the power method and improves the iteration complexity from $\mathcal{O}\left(\frac{1}{\Delta}\log\frac{1}{\epsilon}\right)$ to $\mathcal{O}\left(\frac{1}{\sqrt{\Delta}}\log\frac{1}{\epsilon}\right)$. Here, $\Delta$ is the relative difference between the two eigenvalues of largest magnitude of the relevant matrix, which tends to be very small for many practical data sets. Although polynomial approximation is not explicit in the Lancoz method, its convergence proof is heavily based on Chebyshev polynomials [22]. We remark that the notion of noisy Chebyshev iterations was studied for linear system of equations in [10]; however, the proof technique is not directly applicable to eigenproblems.

1.1.3 Canonical correlation analysis

Recently, a series of papers have appeared which establish fundamental theoretical guarantees and propose efficient algorithms for the CCA problem [5, 17, 9, 31]. One of the first scalable algorithm for solving CCA is AppGrad [17], which is an alternating least-squares method followed by projection steps. However, the algorithm is only guaranteed to converge locally. Inspired by the noisy power method [12], the authors in [9] proposed CCALin, a globally convergent algorithm designed to efficiently approximate the power method by breaking the problem into several subproblems, each of them is solved inexactly in an controllable manner using some of the most advanced convex optimization solvers. The main drawback of CCALin is that it has the same iteration complexity $\mathcal{O}\left(\frac{1}{\Delta}\log\frac{1}{\epsilon}\right)$ as the power method. Recently, an accelerated rate for CCA has been achieved by the use of carefully tuned iterative methods [31, 2]. In particular, for solving top-1 CCA problems, the shift-and-invert (SI) [31] and the LazyCCA [2] methods can achieve a so-called doubly-accelerated rate, i.e., having running time depends on $\sqrt{\tilde{\kappa}}$ and $1/\sqrt{\Delta}$ (see, Table 1). For top-$k$ CCA problems, only LazyCCA, which relies on the shift-and-invert framework coupled with an extensive analysis, can achieve doubly-accelerated rate. In this work, we aim to achieve such doubly-accelerated rate based on a simple modification of the power method. Note that to establish such improved convergence bounds, LazyCCA, SI, and our method require more information about the CCA problem than CCALin.

1.2 Contributions

We propose a linearly convergent algorithm for the GEV and CCA problems which combines the simplicity of the power method with the asymptotic performance of the Lanczos method. Our method—called NAPI—builds on an acceleration technique inspired by polynomial approximation and an efficient inexact matrix-vector product computation via advanced convex optimization methods. Our work is motivated by the work in [10, 9, 21], and the main contributions are summarized as follows:

• •

We explicitly characterize the convergence properties of NAPI. Our result matches the asymptotic iteration complexity of the Lanczos method but with much lower computational cost per iteration. In particular, NAPI attains the worst-case complexity bound $\mathcal{O}\left(\frac{1}{\sqrt{\Delta}}\log\frac{1}{\epsilon}\right)$ which is not improvable in terms of the eigenvalue gap $\Delta$. While the Lanczos method can not operate in a stochastic setting, each iteration of NAPI involves solving a least squares problem inexactly, hence can be handled very efficiently by some of the most advanced convex optimization solvers. By appropriately controlling the errors in solving the least squares problem, and by using a well-designed initialization, NAPI obtains a linear running time in the input-size which allows it to operate on large-scale data sets.

• •

Our result can be used to obtain improved running time for several downstream applications. We demonstrate how NAPI can be applied to the CCA problem, resulting in practical and theoretical improvements over the state-of-the-art algorithms (see, Table 1). Experiments on real-world data sets confirm the effectiveness of our method.

1.3 Notation

The symbols $\mathbbmss{R}$ and $\mathbbmss{N}$ denote the set of real and nonnegative natural numbers, respectively. We use $\left\|{\cdot}\right\|$ to denote the Euclidean norm for vectors and the spectral norm for matrices. For a symmetric positive definite matrix $\mathbf{B}$, we denote $\bm{x}^{\top}\mathbf{B}\bm{y}$ as the $\mathbf{B}$-inner product of $\bm{x}$ and $\bm{y}$ and $\left\|{\bm{x}}\right\|_{\mathbf{B}}=\sqrt{\bm{x}^{\top}\mathbf{B}\bm{x}}$ as the $\mathbf{B}$-norm of $\bm{x}$. For a matrix $\mathbf{X}$, we denote by $\left\|{\mathbf{X}}\right\|_{\mathbf{B}}=\left\|{\mathbf{B}^{1/2}\mathbf{X}}\right\|$ and $\left\|{\mathbf{X}}\right\|_{\mathbf{B},\rm{F}}=\left\|{\mathbf{B}^{1/2}\mathbf{X}}\right\|_{\rm{F}}$ the spectral and Frobenius norms, respectively. The number of nonzero elements of a matrix $\mathbf{A}$ is denoted by $\mathrm{nnz}\left(\mathbf{A}\right)$.

1.4 Preliminaries

Principal angles are a useful tool for studying the convergence of subspaces in iterative eigenvalue methods. We rely on the following definition:

Definition 1** (Principal Angles [33]).**

Let $\mathcal{X}$ and $\mathcal{Y}$ be subspaces of $\mathbbmss{R}^{d}$ of dimension at least $k$. The principal angles $0\leq\theta^{(1)}\leq\theta^{(2)}\cdots\leq\theta^{(k)}\leq\pi/2$ between $\mathcal{X}$ and $\mathcal{Y}$ are defined recursively via:

[TABLE]

subject to $\left\|{\bm{x}}\right\|_{\mathbf{B}}=\left\|{\bm{y}}\right\|_{\mathbf{B}}=1$, $\bm{x}^{\top}\mathbf{B}\bm{x}_{j}=0$, $\bm{y}^{\top}\mathbf{B}\bm{y}_{j}=0$, $j=1,\ldots,i-1$. The vectors $\{\bm{x}_{1},\bm{x}_{2},\ldots,\bm{x}_{k}\}$ and $\{\bm{y}_{1},\bm{y}_{2},\ldots,\bm{y}_{k}\}$ are called the principal vectors. For matrices $\mathbf{X}$ and $\mathbf{Y}$, we use $\theta^{(i)}\left(\mathbf{X},\mathbf{Y}\right)$ to denote the $i$th principal angle between their ranges.

Since we are only interested in the largest principal angle, we make slight abuse of notation and write $\theta\left(\mathbf{X},\mathbf{Y}\right)$ instead of $\theta^{(k)}\left(\mathbf{X},\mathbf{Y}\right)$.

Proposition 1**.**

Let $\mathbf{X}$ and $\mathbf{Y}$ be orthogonal bases for subspaces $\mathcal{X}$ and $\mathcal{Y}$ w.r.t $\mathbf{B}$, respectively. Let $\mathbf{X}_{\bot}$ be an orthogonal basis w.r.t $\mathbf{B}$ for the orthogonal complement of $\mathcal{X}$. Then,

[TABLE]

If, in addition, $\mathbf{X}^{\top}\mathbf{B}\mathbf{Y}$ is invertible, then

[TABLE]

The following lemma presents a nice connection between the principal angles and the distances between subspaces. The proof of this lemma when $\mathbf{B}=\mathbf{I}$ can be found in [29, Section II.4]; its generalization to general positive definite matrices is straightforward, hence omitted here.

Lemma 1** (Distance Between Subspaces).**

Let $\mathbf{X}$ and $\mathbf{Y}$ be orthogonal bases for subspaces $\mathcal{X}$ and $\mathcal{Y}$ w.r.t $\mathbf{B}$, respectively. Let $\mathbf{Q}$ be an orthogonal matrix, then

[TABLE]

2 Polynomial Approximation, Least Squares Solvers and Eigenproblems

2.1 What does polynomial approximation have to do with eigenvalues?

The power method [11] is one of the oldest techniques for computing the largest eigenvalue of a matrix, and the conceptual starting point for many contemporary numerical algorithms. It works as follows: choose an initial vector $\bm{x}_{0}\in\mathbbmss{R}^{d}$ uniformly at random from the unit sphere and multiply it repeatedly by the symmetric matrix $\mathbf{A}\in\mathbbmss{R}^{d\times d}$. This generates a sequence of iterates $\{\bm{x}_{s}\}$, related via $\bm{x}_{s}=\mathbf{A}\bm{x}_{s-1}=\mathbf{A}^{s}\bm{x}_{0}$. Under certain conditions, this sequence, normalized appropriately, will converge to the leading eigenvector of $A$, i.e, the one associated with the eigenvalue of largest magnitude.

The power method can be viewed as applying of a degree-$s$ polynomial, namely $p_{s}\left(z\right)=z^{s}$, to $\mathbf{A}$, and then multiplying it by the initial vector $\bm{x}_{0}$. The polynomial $p_{s}\left(z\right)$ evaluates $1$ to $1$ and moves every $z$ with $\left|{z}\right|<1$ closer to [math] at the rate determined by $\left|{z}\right|$. Let $\mathbf{A}=\mathbf{U}\mathbf{\Lambda}\mathbf{U}^{\top}$ be the eigenvalue value decomposition of $\mathbf{A}$, where $\mathbf{\Lambda}=\mathrm{diag}\left({\lambda_{1}},{\ldots},{\lambda_{d}}\right)$ is a matrix of eigenvalues satisfying $\left|{\lambda_{1}}\right|>\left|{\lambda_{2}}\right|\geq\ldots\geq\left|{\lambda_{d}}\right|$. Then, the polynomial $p_{s}\left(z\right)$ applied to $\mathbf{A}$ is simply

[TABLE]

Since $\left|{\lambda_{1}}\right|>\left|{\lambda_{2}}\right|$, $\mathbf{A}^{s}/\lambda_{1}^{s}$ converges to $\bm{u}_{1}\bm{u}_{1}^{\top}$, where $\bm{u}_{1}$ is the first column of $\mathbf{U}$. This means that an estimate of $\bm{u}_{1}$ can be found by simply picking the first column of $\mathbf{A}^{s}$ and normalizing it to have a unit norm. If we suppose that $\left|{\lambda_{2}}\right|=\left(1-\Delta\right)\left|{\lambda_{1}}\right|$, then it will take roughly $s=\mathcal{O}\left(\frac{1}{\Delta}\log\frac{1}{\epsilon}\right)$ iterations until $|\lambda_{2}|^{s}\leq\epsilon|\lambda_{1}|^{s}$. Since the eigengap $\Delta$ tends to be very small for practical large-scale matrices, the power method converges very slowly. Fortunately, the next proposition indicates that it is possible to to improve the running time from $\mathcal{O}\left(\frac{1}{\Delta}\log\frac{1}{\epsilon}\right)$ to $\mathcal{O}\left(\frac{1}{\sqrt{\Delta}}\log\frac{1}{\epsilon}\right)$ by a simple modification of the power method.

Proposition 2**.**

[23, Theorem 3.3]** For any positive integers $s$ and $d$, there is a degree-$d$ polynomial $p_{s,d}\left(x\right)$ satisfying

[TABLE]

The result implies that it is possible to approximate $x^{s}$ to any accuracy $\delta>0$ i.e.

[TABLE]

using a polynomial $p_{s,d}(x)$ of degree $d=\lceil\sqrt{2\ln\left(2/\delta\right)s}\rceil$. This suggests that it should be possible to develop an approximate power iteration which reduces the number of matrix multiplications by $A$ from $s$ to $\sqrt{s}$. Of course, the existence result in Proposition 2 is only useful if there is an efficient way to construct the approximating polynomial. The proof of the proposition relies on Chebyshev polynomials, which are ubiquitous in numerical optimization. We also note that the result above is essentially optimal in the sense that polynomial approximations of $x^{s}$ over $\left[-1,1\right]$ require degree $\Omega(\sqrt{s})$ [23, Chapter 5].

2.2 Stochastic optimization of least squares

Over the past few years, there has been a great attention in developing stochastic algorithms for solving empirical risk minimization problems in machine learning. Specifically, one is interested in the unconstrained minimization of

[TABLE]

where $f_{1},\ldots,f_{n}$ is $L$-smooth and convex functions, and $f$ is $\mu$-strongly convex function. Solving this problem using standard first order methods such as gradient descent requires $\mathcal{O}\left(n\kappa\log\frac{1}{\epsilon}\right)$ passes over the data set to achieve an $\epsilon$-optimal solution, where the condition number $\kappa=L/\mu$ can be as large as $\Omega\left(n\right)$ in machine learning applications [4]. This running time can be further improved to $\mathcal{O}\left(n\sqrt{\kappa}\log\frac{1}{\epsilon}\right)$ if one uses Nesterov’s accelerated scheme [18]. Many linearly convergent stochastic methods, such as SDCA [25], SAGA [6], or SVRG [15], have been introduced and shown to outperform the standard first order methods under mild assumptions. These methods only access to the gradient of one individual component function $f_{i}$ at each step of the algorithm, instead of the full gradient as in standard first order methods. This results in an improved computational complexity of $\mathcal{O}\left(\left(n+\kappa\right)\log\frac{1}{\epsilon}\right)$ passes over the data set to achieve an $\epsilon$-optimal solution in expectation. When these methods are combined with Nesterov acceleration, the expected complexity becomes $O\left((n+\sqrt{n\kappa})\log(1/\epsilon)\right)$ (see, e.g., [7] and [1]).

In this work, we approximately compute the matrix-vector product $\mathbf{B}^{-1}\bm{b}$ by minimizing the least-squares loss

[TABLE]

over $\bm{w}$. This problem has the unique solution is $\bm{w}^{\star}=\mathbf{B}^{-1}\bm{b}$. In our applications of interest, $\mathbf{B}$ is usually an empirical covariance matrix of the form $\frac{1}{n}\sum_{i=1}^{n}\bm{x}_{i}\bm{x}_{i}^{\top}$, hence $f\left(\bm{w}\right)$ can be expressed as a finite-sum as in (3) with

[TABLE]

Thus, we can find an approximate solution to $\mathbf{B}^{-1}\bm{b}$ using either standard first-order methods or more advanced stochastic solvers. As we will see, using the advanced solvers together with good initializations will lead to a significant improvement in iteration complexity of our method.

2.3 The top-k generalized eigenvector problem

Definition 2** (Top-$k$ generalized eigenvectors).**

*Given symmetric matrices $\mathbf{A}$ and $\mathbf{B}$ where $\mathbf{B}$ is positive definite, the top-$k$ generalized eigenvectors $\bm{w}_{1},\ldots,\bm{w}_{k}$ are defined by *

[TABLE]

Problem (4) can be reduced to the problem of computing the top-$k$ eigenvectors associated with $k$ eigenvalues of largest magnitude of $\mathbf{M}=\mathbf{B}^{-1/2}\mathbf{A}\mathbf{B}^{-1/2}$ and then multiplying by $\mathbf{B}^{-1/2}$ (see, Lemma 2). Let us use the simple power method for this task. Let $\bm{y}=\mathbf{B}^{-1/2}\mathbf{M}^{T}\mathbf{x}$ be the result of $T$ iterations of the power method followed by multiplying with $\mathbf{B}^{-1/2}$, we can write $\bm{y}=\left(\mathbf{B}^{-1}\mathbf{A}\right)^{T}\mathbf{B}^{-1/2}\bm{x}$. Since the power method converges from a random initial vector $\bm{x}$ with high probability, one can ignore $\mathbf{B}^{-1/2}$ and compute instead $\bm{y}=\left(\mathbf{B}^{-1}\mathbf{A}\right)^{T}\bm{x}$. Thus, we can compute our desired eigenvectors by simply alternating between applying $\mathbf{A}$ and $\mathbf{B^{-1}}$ to a random initial vector. The authors of [9] suggested a procedure (GenELin) for this task which applies $\mathbf{B}^{-1}$ approximately by least squares solvers. It is easy to see that GenLin simulates the noisy power method in [12], and hence it will converge as long as the noise is under control.

3 Computing the Leading Generalized Eigenvector

In this section, we introduce NAPI, an noisy accelerated power method for solving (1). We then characterize the convergence rate of the proposed algorithm for the special case of computing the leading generalized eigenvector.

3.1 Description of the algorithm

Our algorithm is designed to run in an inner-outer fashion and it hinges upon: (i) an acceleration technique inspired by Chebyshev polynomials for approximating matrix powers; and (ii) an efficient inexact matrix-vector product computation. In particular, the acceleration technique in the outer loop allows to reduce the number of outer iterations, while the subproblems in the inner loop are efficiently solved by some of the most advanced convex optimization solvers. These solvers are capable of exploiting sparsity and problem structure and scale to truly large-scale data sets. The last critical component of our algorithm is a well-designed warm start procedure which reduces the practical running times even further.

3.2 Deriving acceleration through Chebyshev magic

As discussed above, the power method can be accelerated by finding a good low-degree polynomial approximation to $x^{T}$. However, to evaluate a polynomial of degree $T$, we may need to store up to $T$ terms before adding them up. Fortunately, the next result demonstrates that we can construct a Chebyshev polynomial recursively using only the two previous terms.

Definition 3** (Chebyshev Polynomials).**

For a nonnegative integer $t$, the degree $t$ Chebyshev polynomial of the first kind, $T_{t}\left(z\right)$, is defined recursively as follows:

[TABLE]

To motivate our algorithm, let us focus on the unconstrained version of Problem (4), where we relax the condition that the solution should have unit $\mathbf{B}$-norm. Our goal is to find the leading eigenvector of $\mathbf{B}^{-1}\mathbf{A}$ by approximating $\left(\mathbf{B}^{-1}\mathbf{A}\right)^{T}$ using a Chebyshev polynomial. Using Definition 3 with $z=\mathbf{B}^{-1}\mathbf{A}$ and letting $\bm{y}_{t}=T_{t}\left(\mathbf{B}^{-1}\mathbf{A}\right)\bm{w}_{0}$, we can approximate $\left(\mathbf{B}^{-1}\mathbf{A}\right)^{T}\bm{w}_{0}$ by recursively building a new iterate using only the previous two iterates:

[TABLE]

Note that the result in Proposition 2 is stated for approximating over the interval $[-1,1]$. Therefore, we will consider a scaled Chebyshev polynomial of the form $C_{t}(\mathbf{B}^{-1}\mathbf{A})=T_{t}\left(\frac{\mathbf{B}^{-1}\mathbf{A}}{c}\right)$ for some constant $c$. Letting $\bm{x}_{t}=C_{t}(\mathbf{B}^{-1}\mathbf{A})\bm{w}_{0}$, we will thus study the recursion

[TABLE]

where $\beta$ is a constant depending on $c$. As we shall see, determining the best possible $\beta$ leads to an improved convergence rate for our method.

Returning to the constrained case, we wish to use a (scaled) Chebyshev polynomial to get an improved convergence rate while ensuring that the iterates have unit $\mathbf{B}$-norm. To this end, in our algorithm, we normalize the intermediate iterate $\tilde{\bm{w}}_{t+1}$ to have unit $\mathbf{B}$-norm and rescale the current iterate $\bm{w}_{t}$ by the same factor as that for $\tilde{\bm{w}}_{t+1}$. The iterates generated by Algorithm 1 can be written as

[TABLE]

and the latest iterate $\bm{w}_{T}$ can be expressed as

[TABLE]

where $C_{T}$ is a scaled Chebyshev polynomial of degree $T$. That is, the output vector obtained by applying (7)-(8) recursively is equivalent to performing only recursion (6) without normalizing the iterates at each step, but only at the last step. To see that, suppose that $\bm{w}_{-1}=\bm{x}_{-1}=0$ and $\bm{w}_{0}=\bm{x}_{0}$, then we will show by induction that $\bm{w}_{t}=\bm{x}_{t}\alpha_{t}^{-1}$ for all $t\geq 1$, where $\alpha_{t}=\gamma_{t}\gamma_{t-1}\ldots\gamma_{0}>0$. Assume that $\bm{w}_{i}=\bm{x}_{i}\alpha_{i}^{-1}$ for all $i\leq t$, the sequence generated by (7)-(8) can be written as

[TABLE]

It follows that

[TABLE]

as desired. The following result, whose proof can be found in [21, Theorem 8], characterizes the convergence rate of Algorithm 1 in the absence of noise.

Proposition 3** (Convergence of the exact method).**

Let $\lambda_{i}$ be the eigenvalues of $\mathbf{B}^{-1}\mathbf{A}$ satisfying $\left|{\lambda_{1}}\right|>\left|{\lambda_{2}}\right|\geq\ldots\geq\left|{\lambda_{d}}\right|$, and let $\Delta=1-\left|{\lambda_{2}}\right|/\left|{\lambda_{1}}\right|$ be its relative eigenvalue gap. If the least squares subproblems in Step 4 of Algorithm 1 are solved exactly and the parameter $\beta$ satisfies $\left|{\lambda_{2}}\right|\leq 2\sqrt{\beta}<\left|{\lambda_{1}}\right|$, then

[TABLE]

Moreover, if $\beta=\lambda_{2}^{2}/4$, then one can achieve an $\epsilon$-optimal solution, i.e, $\sin\theta\left(\bm{w}_{T},\bm{u}_{1}\right)\leq\epsilon$, after at most $T=\mathcal{O}\left(\frac{1}{\sqrt{\Delta}}\log\frac{\tan\theta_{0}}{\epsilon}\right)$ iterations.

The above result matches the iteration complexity of the Lanczos method, a significantly more complex algorithm [11]. However, note that the result requires exact computation of $\mathbf{B}^{-1}\mathbf{A}\bm{w}$ for some vector $\bm{w}$, which (as already discussed) is prohibitively expensive for large scale data sets.

3.3 Warm-start

Initialization plays a critical role in designing an efficient method, cf. the work on the standard power iteration in [9]. Recall that the matrix-vector product $\mathbf{B}^{-1}\mathbf{A}\bm{w}_{t}$ is computed approximately by solving the least squares problem

[TABLE]

whose unique solution is $\bm{w}^{\star}=\mathbf{B}^{-1}\mathbf{A}\bm{w}_{t}$. As the iterates converge, the solutions to the least squares problems in consecutive iterations become increasingly similar. It is natural to use a scaled version of the previous iterate as initialization. This gives the following initial objective

[TABLE]

Minimizing $f\left(\alpha\bm{w}_{t}\right)$ over $\alpha$ gives us the optimal scaling $\alpha_{t}^{\star}=\bm{w}_{t}^{\top}\mathbf{A}\bm{w}_{t}/\bm{w}_{t}^{\top}\mathbf{B}\bm{w}_{t}$. As will be shown later, with this initialization, we only need to run the least squares solver until the initial error has been reduced by a constant factor, independent of the final suboptimality required.

3.4 Convergence argument

In this subsection, we will show that as long as the errors in the matrix-vector multiplications are controlled in an appropriate way, NAPI enjoys the same convergence rate as the idealized method where subproblems are solved exactly.

Due to the inexact solution of subproblems, the output of Algorithm 1 is no longer equivalent to applying the Chebyshev polynomial to an initial vector, hence the proof of Proposition 3 is not applicable anymore. However, the iterates defined by Algorithm 1 can be seen as a second order iterative system and be cast into the following form

[TABLE]

where $\bm{e}_{t}=\tilde{\bm{w}}_{t+1}-\mathbf{B}^{-1}\mathbf{A}\bm{w}_{t}$ is the approximation error in Step 4 of Algorithm 1,

[TABLE]

To establish the convergence rate of NAPI, we will thus study the behaviour of this extended system. . Note that the extended vectors $\bm{\vartheta}_{t}$ do not have unit norm w.r.t any fixed matrix, hence care must be taken when performing the convergence analysis for this sequence. We begin with the existence of eigenvectors and eigenvalues of $\mathbf{B}^{-1}\mathbf{A}$ and $\mathbf{C}$.

Lemma 2** (Existence of Eigenbasis).**

Let $\left(\bm{p}_{i},\lambda_{i}\right)$ be the eigenpairs of the symmetric matrix $\mathbf{B}^{-1/2}\mathbf{A}\mathbf{B}^{-1/2}$. Then, $\left(\bm{u}_{i}=\mathbf{B}^{-1/2}\bm{p}_{i},\lambda_{i}\right)$ is an eigenpair of $\mathbf{B}^{-1}\mathbf{A}$. In addition, $\left([\mu_{i}\bm{u}_{i}^{\top},\bm{u}_{i}^{\top}]^{\top},\mu_{i}\right)$ is an eigenpair of $\mathbf{C}$, where the eigenvalues $\mu_{i}$ satisfy

[TABLE]

Proof.

See Appendix A.1. ∎

The next lemma shows that under a certain condition on the suboptimality in solving the least-squares subproblems, the angle between $\bm{\vartheta}_{t}$ and the optimal eigenvector $\bm{v}_{1}$ shrinks geometrically in each iteration.

Lemma 3**.**

Let $\bm{v}_{1}$ be a leading eigenvector of the augmented matrix $\mathbf{C}$ and let $\mu_{i}$ be the eigenvalues of $\mathbf{C}$ such that $\left|{\mu_{1}}\right|>\left|{\mu_{2}}\right|\geq\ldots\geq\left|{\mu_{2d}}\right|$. If the errors in solving the least squares problems satisfy

[TABLE]

for all $t\in\mathbbmss{N}$, then

[TABLE]

Proof.

See Appendix A.2. ∎

We are now ready to state our main theorem, quantifying the iteration complexity of NAPI:

Theorem 1**.**

Let $\lambda_{i}$ be the eigenvalues of $\mathbf{B}^{-1}\mathbf{A}$ satisfying $\left|{\lambda_{1}}\right|>\left|{\lambda_{2}}\right|\geq\ldots\geq\left|{\lambda_{d}}\right|$, and let $\Delta=1-\left|{\lambda_{2}}\right|/\left|{\lambda_{1}}\right|$ be its relative eigenvalue gap. For $\left|{\lambda_{2}}\right|\leq 2\sqrt{\beta}<\left|{\lambda_{1}}\right|$, the output of Algorithm 1 satisfies

[TABLE]

where $\theta_{0}=\theta\left(\bm{\vartheta}_{0},\bm{v}_{1}\right)$. Furthermore, if the parameter $\beta$ is set to $\left|{\lambda_{2}}\right|^{2}/4$, we can achieve an $\epsilon$-optimal solution after at most

[TABLE]

iterations. The total running time of the algorithm is

[TABLE]

where $\mathcal{T}\left(\tau\right)$ is the time that a least squares solver takes to reduce the residual error by a factor $\tau$.

Proof.

See Appendix B ∎

Remark 1**.**

The result in Theorem 1 for the extended system readily implies the convergence of $\bm{w}_{t}$ to an optimal solution to Problem 4. More specifically, invoking Lemma 1 with

[TABLE]

yields

[TABLE]

It follows that

[TABLE]

Since $\left\|{\bm{w}_{T}}\right\|_{\mathbf{B}}=\left\|{\bm{u}_{1}}\right\|_{\mathbf{B}}=1$, the above inequality indicates that $\bm{w}_{T}$ converges to a leading eigenvector ($\bm{u}_{1}$ or $-\bm{u}_{1}$) of $\mathbf{B}^{-1}\mathbf{A}$, an optimal solution of Problem 4.

Remark 2**.**

Theorem 1 shows that as long as the errors arising from solving the least squares problems approximately are controlled appropriately, the proposed algorithm enjoys the same convergence rate as the idealized method(c.f. Proposition 3). Oddly, our analysis shows that the accelerated method is less sensitive to errors than the standard power method in [9]. In particular, the condition (10) in Lemma 3 requires the noises to scale with $\sqrt{\Delta}$ while the condition in [9] requires the noises scale with $\Delta$. It is worth mentioning that compared to [9], our method has one additional parameter (the momentum gain), whose optimal value is not known in general. This situation is quite common in momentum-like methods such as the heavy-ball method [20] in convex optimization. However, it has been widely observed that even with a suboptimal estimate of the momentum parameter, one often can achieve significantly speed-ups over the original algorithm.

Corollary 1**.**

If we use Nesterov’s accelerated method as a least squares solver, the running time of Algorithm 1 to achieve an $\epsilon$-optimal solution is bounded by

[TABLE]

*where *

[TABLE]

This is a linear running time in the input size of the problem, i.e., the number of nonzero elements of $\mathbf{A}$ and $\mathbf{B}$.

3.5 Extension to the Top-$k$ Generalized Eigenvectors

We consider an extension of our accelerated method to the problem of finding the top-$k$ generalized eigenvectors. The algorithm is formally given as Algorithm 2. It is a natural generalization of Algorithm 1, where the iterates are now matrices of size $d\times k$, and the normalization steps are replaced by orthogonalization steps using the Gram-Schmidt procedure with the inner product $\left\langle\cdot,\cdot\right\rangle_{\mathbf{B}}$. Similarly as in the top-$1$ case, we will study the convergence behaviour of the sequence generated by Algorithm 2 via the extended system given by

[TABLE]

Let $\mathbf{U}=[\bm{u}_{1},\ldots,\bm{u}_{k}]$ with $\mathbf{U}^{\top}\mathbf{B}\mathbf{U}=\mathbf{I}$ be the matrix of top-$k$ eigenvectors of $\mathbf{B}^{-1}\mathbf{A}$ and let $\lambda_{1},\ldots,\lambda_{k}$ be the corresponding top-$k$ eigenvalues. Then, by Lemma 2, the matrices of top-$k$ eigenvectors $\bar{\mathbf{V}}$ and eigenvalues $\mathbf{\Lambda}$ of the extended matrix $\mathbf{C}$ are $\bar{\mathbf{V}}=[\mathbf{\Lambda}\mathbf{U}^{\top},\mathbf{U}^{\top}]^{\top}$ and $\mathbf{\Lambda}=\mathrm{diag}\left(\mu_{1},\ldots,\mu_{k}\right)$ with $\mu_{i}$ being solutions to the equations $\mu_{i}^{2}-\lambda_{i}\mu_{i}+\beta=0$, $i=1\ldots,k$. Let $\mathbf{V}=\bar{\mathbf{V}}(\mathbf{I}+\mathbf{\Lambda}^{2})^{-1/2}$, and let

[TABLE]

then $\mathbf{V}^{\top}\mathbf{\mathcal{B}}\,\mathbf{V}=\mathbf{I}_{k}$. Furthermore, using QR decomposition w.r.t $\mathbf{\mathcal{B}}$, one can write $\mathbf{Y}_{t}$ as $\mathbf{Y}_{t}=\bar{\mathbf{Y}}_{t}\mathbf{G}_{t}$ for some matrix $\bar{\mathbf{Y}}_{t}$ satisfying $\bar{\mathbf{Y}}_{t}^{\top}\mathbf{\mathcal{B}}\,\bar{\mathbf{Y}}_{t}=\mathbf{I}$ and for some nonsingular matrix $\mathbf{G}_{t}\in\mathbbmss{R}^{k\times k}$.

The following theorem characterizes the convergence properties of the above extended system. Note that the angles in the theorem are measured w.r.t $\mathbf{\mathcal{B}}$.

Theorem 2**.**

*Let $\lambda_{i}$ be the eigenvalues of $\mathbf{B}^{-1}\mathbf{A}$ satisfying $\left|{\lambda_{1}}\right|\geq\ldots\geq\left|{\lambda_{k}}\right|>\left|{\lambda_{t+1}}\right|\geq\ldots\geq\left|{\lambda_{d}}\right|$, and let $\Delta_{k}=1-\left|{\lambda_{k+1}}\right|/\left|{\lambda_{k}}\right|$ be its relative eigenvalue gap. If the errors in solving the least squares problems satisfy *

[TABLE]

for all $t\in\mathbbmss{N}$, and if $\left|{\lambda_{k+1}}\right|\leq 2\sqrt{\beta}<\left|{\lambda_{k}}\right|$, the output of Algorithm 2 satisfies

[TABLE]

where $\theta_{0}=\theta(\bar{\mathbf{Y}}_{0},\mathbf{V})$. Furthermore, if the parameter $\beta$ is set to $\left|{\lambda_{k+1}}\right|^{2}/4$, we can achieve an $\epsilon$-optimal solution after at most

[TABLE]

iterations. The total running time of the algorithm is

[TABLE]

where $\gamma=\lambda_{1}/\lambda_{k}$ and $\mathcal{T}\left(\tau\right)$ is the time a least squares solver takes to reduce the residual error by a factor $\tau$.

Proof.

See Appendix C. ∎

For $k=1$, we have shown in Remark 1 that the convergence of the extended sequence translates to the convergence of the original one. However, for a general case, we need a more involved analysis to derive a similar conclusion. We start by invoking Lemma 1 with $\mathbf{X}=\bar{\mathbf{Y}}_{T}$, $\mathbf{Y}=\mathbf{V}$, $\mathbf{B}=\mathbf{\mathcal{B}}$ to obtain

[TABLE]

It follows that there exists an orthogonal matrix $\mathbf{Q}$ such that

[TABLE]

Therefore, if we let $\mathbf{M}=\mathbf{\Lambda}(\mathbf{I}+\mathbf{\Lambda}^{2})^{-1/2}\mathbf{Q}$, then $\mathbf{W}_{T}$ can be written as

[TABLE]

where $\mathbf{E}$ is a perturbation matrix satisfying $\left\|{\mathbf{E}}\right\|_{\mathbf{B}}\leq 2\sin\frac{\theta\left(\bar{\mathbf{Y}}_{T},\mathbf{V}\right)}{2}$. Note that since $\mathbf{M}\mathbf{G}_{T}^{-1}$ is nonsingular, the range space of $\mathbf{W}_{T}$ is identical to that of $\mathbf{U}+\mathbf{E}\mathbf{M}^{-1}$. Consider the following matrix

[TABLE]

then it can be verified that $\mathbf{P}^{\top}\mathbf{B}\mathbf{P}=\mathbf{I}$, thereby being an orthogonal basis for the range space of $\mathbf{W}_{T}$. Using the binomial expansion identity [13]

[TABLE]

for $\left\|{\mathbf{A}}\right\|<1$ and the fact that $\mathbf{U}^{\top}\mathbf{B}\mathbf{U}=\mathbf{I}$, one can write the matrix $\mathbf{P}$ as \mathbf{P}=\big{(}\mathbf{U}+\mathbf{E}\mathbf{M}^{-1}\big{)}\big{(}\mathbf{I}+\mathbf{E}^{\prime}\big{)} with $\left\|{\mathbf{E}^{\prime}}\right\|_{\mathbf{B}}=O(\left\|{\mathbf{E}}\right\|_{\mathbf{B}})$. Therefore, if we denote by $\mathbf{U}_{\perp}$ an orthogonal basis (w.r.t $\mathbf{B}$) of the subspace spanned by $\bm{u}_{k+1},\ldots,\bm{u}_{d}$, then it holds that

[TABLE]

which implies that

[TABLE]

for some constant $c$. The last inequality indicates that the original sequence converges at the same rate as the extended one, as desired.

4 Application to CCA

Our main result can readily lead to improvement of several downstream applications. In this section, we will discuss the benefits of our proposed method in solving the CCA problem.

In CCA, we are provided with pairs of data points from two views $(\bm{x}_{1},\bm{y}_{1}),\ldots,(\bm{x}_{n},\bm{y}_{n})$ where $\bm{x}_{i}\in\mathbbmss{R}^{1\times d_{1}}$, $\bm{y}_{i}\in\mathbbmss{R}^{1\times d_{2}}$, and $n$ is the size of the training set. For example, the views can be the visual image data and audio data in a video of people speaking. Let $\mathbf{X}=\left[\bm{x}_{1},\ldots,\bm{x}_{n}\right]^{\top}\in\mathbbmss{R}^{n\times d_{1}}$ and $\mathbf{Y}=\left[\bm{y}_{1},\ldots,\bm{y}_{n}\right]^{\top}\in\mathbbmss{R}^{n\times d_{2}}$ be the data matrices for each view, respectively. We define the relevant empirical covariance matrices as $\mathbf{\Sigma}_{11}=\frac{1}{n}\mathbf{X}^{\top}\mathbf{X}+\gamma_{1}\mathbf{I}$, $\mathbf{\Sigma}_{22}=\frac{1}{n}\mathbf{Y}^{\top}\mathbf{Y}+\gamma_{2}\mathbf{I}$, and $\mathbf{\Sigma}_{12}=\frac{1}{n}\mathbf{X}^{\top}\mathbf{Y}$, where $\gamma_{1}$ and $\gamma_{2}$ are regularization parameters. The goal is to find a linear transformation of the points in each view that retains as much as the redundant information between the views. More specifically, CCA finds the solution of the following problem

[TABLE]

It is easy to check that any stationary point $\left(\bm{\phi}_{i},\bm{\varphi}_{i}\right)$ of Problem (11) satisfies

[TABLE]

where $\lambda_{i}$ and $\lambda^{\prime}_{i}$ are two constants. Using the constraints, it can be verified that $\lambda_{i}=\lambda^{\prime}_{i}$. The system (12) above can therefore be cast into the following generalized eigenproblem

[TABLE]

The solution to this problem has $d_{1}+d_{2}$ eigenvalues $\lambda_{1}>\lambda_{2}>\ldots>-\lambda_{2}>-\lambda_{1}$; for each eigenvalue $\lambda_{i}$ with the corresponding eigenvector $[\bm{\phi}_{i}^{\top},\bm{\varphi}_{i}^{\top}]^{\top}$, $-\lambda_{i}$ is also an eigenvalue with the corresponding eigenvector $[\bm{\phi}_{i}^{\top},-\bm{\varphi}_{i}^{\top}]^{\top}$. Let $\mathbf{A}$ and $\mathbf{B}$ be the matrices on the left and right in (13), respectively, then the eigenvalue gap of $\mathbf{B}^{-1}\mathbf{A}$ is zero. Thus, iterative methods for finding the leading eigenvector of $\mathbf{B}^{-1}\mathbf{A}$ may not converge. Since the top $2k$-dimensional eigen-space of $\mathbf{B}^{-1}\mathbf{A}$ is the subspace spanned by $[\bm{\phi}_{i}^{\top},\bm{\varphi}_{i}^{\top}]^{\top}$ and $[\bm{\phi}_{i}^{\top},-\bm{\varphi}_{i}^{\top}]^{\top}$, $i=1,\ldots,k$, one can solve the above problem for the top-$2k$ eigenvector problem and then choose any orthogonal basis that spans the output subspace and a random $k$-dimensional projection of those vectors. Finally, to fulfil the constraints in CCA problem, we perform the last orthogonalization steps w.r.t $\mathbf{\Sigma}_{11}$ and $\mathbf{\Sigma}_{22}$, respectively. The formal description of the discussions above are given in Algorithm 3.

We remark that when running Algorithm 3, we do not need to form $\mathbf{A}$ and $\mathbf{B}$ explicitly thanks to their diagonal forms. In fact, we only need to apply $\mathbf{\Sigma}_{11}$, $\mathbf{\Sigma}_{12}$, and $\mathbf{\Sigma}_{22}$ as operators, which can be done efficiently by applying $\mathbf{X}$ and $\mathbf{Y}$ appropriately. For example, at some points of the algorithm, one needs to compute $\bm{x}^{\top}\mathbf{\Sigma}_{11}\bm{y}$ for some vectors $\bm{x},\bm{y}$, one can compute $\mathbf{X}\bm{y}$ and $\mathbf{X}\bm{x}$ and then compute the inner product of those vectors without forming $\mathbf{\Sigma}_{11}$. This also allows to exploit the sparsity of the data matrices $\mathbf{X}$ and $\mathbf{Y}$, which is critical for large-scale problems.

The least squares problem at the $t$th iteration of procedure NAPI can be written explicitly as follows

[TABLE]

The condition number of this problem is

[TABLE]

If one wishes to use stochastic optimization solvers, one can be express the objective function above as

[TABLE]

where

[TABLE]

In this case, the relevant condition number depends on the individual Lipschitz constant of the component functions, and can be computed as

[TABLE]

The following lemma shows that Steps 5 and 6 in Algorithm 3 do not cause much loss in the alignment with the true canonical space.

Lemma 4**.**

[9, Lemma 14]** If the output of procedure NAPI has the angle at most $\theta$ with the true top-$2k$ generalized eigenspace of $\mathbf{B}$, $\mathbf{A}$, then with probability at least $1-\delta$, both $\bm{\Phi}_{T}$ and $\bm{\Psi}_{T}$ has angle at most $\mathcal{O}\left(k^{2}\theta/\delta^{2}\right)$ with the true top-$k$ canonical space of $\mathbf{X}$ and $\mathbf{Y}$.

Finally, we state the following theorem charactering the iteration complexity of Algorithm 3.

Theorem 3**.**

*Given data matrices $\mathbf{X}$ and $\mathbf{Y}$, let $\mathbf{A}$ and $\mathbf{B}$ be the matrices defined in Algorithm 3. Let $\lambda_{i}$ be the eigenvalues of $\mathbf{B}^{-1}\mathbf{A}$ satisfying $\lambda_{1}>\lambda_{2}>\ldots>-\lambda_{2}>-\lambda_{1}$, and let $\Delta_{k}=1-{\left|{\lambda_{k+1}}\right|}/{\left|{\lambda_{k}}\right|}$ be its relative eigenvalue gap. Let $\bm{\Phi}^{\star}$ and $\bm{\Psi}^{\star}$ be the top-$k$ true canonical space of $\mathbf{X}$ and $\mathbf{Y}$, respectively. If $\beta=\left|{\lambda_{k+1}}\right|^{2}/4$, then with probability at least $1-\delta$, the output of Algorithm 3 satisfying *

[TABLE]

in time

[TABLE]

where $\mathrm{nnz}\left(\mathbf{X},\mathbf{Y}\right)=\mathrm{nnz}\left(\mathbf{X}\right)+\mathrm{nnz}\left(\mathbf{Y}\right)$, $\gamma=\lambda_{1}/\lambda_{k}$ and $c$ is a universal constant.

Proof.

The proof follows readily from Lemma 4 and Theorem 2. ∎

Corollary 2**.**

If we use Nesterov’s accelerated method as a least squares solver, the running time of Algorithm 3 to achieve an $\epsilon$-optimal solution is bounded by

[TABLE]

On the other hand, if we choose to use SVRG, then the running time is bounded by

[TABLE]

where

[TABLE]

5 Experimental Results

In this section, we perform experiments on real-word data sets to validate the effectiveness of the proposed algorithms. In the plots, we illustrate the performance of NAPI-based CCA on the following real-world data sets.

Mediamill is an annotated video data set from the Mediamill Challenge for automated detection of semantic concepts, which contains 85 hours of international broadcast news data [28]. These news videos are first automatically segmented into 30,993 subshots. Each image is a representative key-frame of a video subshot annotated with 101 labels and consists of 120 features. CCA is used to learn the correlation structure between the images and its labels.

MNIST is a data set of 60,000 handwritten digits from 0 to 9 [16]. Each digit is represented by an image of $392\times 392$ values in $\left[0,1\right]$. CCA is performed to explore the correlated representations between the left and right halves of the handwritten digit images.

XRMB is Wisconsin Xray Microbeam database consits of simultaneous acoustic and articulatory recordings [32, 30]. The inputs to CCA are the acoustic and articulatory features concatenated over a 7-frame window around each frame, giving acoustic vectors $\bm{x}\in\mathbbmss{R}^{273}$ and articulatory vectors $\bm{y}\in\mathbbmss{R}^{112}$.

Figure 1 shows the suboptimality in objective versus the number of epochs (passes over the full data set) for different algorithms solving CCA problem with $k=1$ starting from a random initial vector. We compare the following algorithms: AppGrad [17], ALS-VR [31], CCALin [9], and SI-VR [31] where SVRG is used as the least squares solver. For CCALin, and NAPI, we set the SVRG stepsizes to be $1/\max_{i}\left(\left\|{\bm{x}_{i}}\right\|^{2},\left\|{\bm{y}_{i}}\right\|^{2}\right)$. We use the implementation of SI-VR provided by the authors in [31]. This implementation uses slightly different parameters than analyzed in [31], but outperforms the theoretically suggested values. For AppGrad, we tune the stepsizes and report the result from using the best ones. For each dataset, we vary the regularization parameters to influence the condition numbers of the least squares problems, where larger regularization implies better conditioning.

It can be seen that our proposed method converges linearly to an optimal solution and significantly outperforms the other algorithms, especially for ill-conditioned problems. As predicted by the theoretical results, NAPI and SI-VR have superior performance compared to other non-accelerated methods. It should be emphasized that for well-conditioned problems such as the case of the Mediamill data set with $\gamma_{x}=\gamma_{y}=10^{-3}$, CCALin and AlS-VR can have comparable performance to NAPI and SI-VR. Finally, we observe that for ill-conditioned problems, the convergence of SI-VR can be quite slow with a similar slope as ALS-VR, while NAPI still maintains an accelerated rate and quickly finds a high accuracy solution.

6 Conclusion

We have proposed and analyzed a simple algorithm for finding the dominant generalized eigenvectors of a pair of symmetric matrices. The algorithm admits a linear convergence rate and was shown to outperform the current state-of-the-art algorithms. To achieve that result, our algorithm exploits the problem structure and decomposes the overall problem into a sequence of subproblems. The approximation theory is used to accelerate the overall process which helps to achieve the asymptotic iteration complexity of the Lanczos method, while each of the subproblem can be solved efficiently by advanced iterative solvers and with a well-design initialization. This translates into a linear time algorithm in the input sizes, suitable for many large-scale applications in machine learning and data science.

Acknowledgment

This research was sponsored in part by the Knut and Alice Wallenberg Foundation and the Swedish Research Council. We thank Weiran Wang for sharing with us the implementation of the shift-and-invert method in [31] and two anonymous reviewers for their useful comments and suggestions.

Appendix A Proof of Auxiliary Lemmas

A.1 Proof of Lemma 2

For the first part, we have

[TABLE]

For the second part, for any eigenpair $\left(\bm{u},\lambda\right)$ of $\mathbf{B}^{-1}\mathbf{A}$, let $\mu$ be a solution of $\mu^{2}-\lambda\mu+\beta=0$. Then,

[TABLE]

Therefore, the vector $\bm{v}=[\mu\bm{u}^{\top},\bm{u}^{\top}]^{\top}$ is an eigenvector of $\mathbf{C}$ with eigenvalue $\mu$. The proof is complete.

A.2 Proof of Lemma 3

The proof follows closely those in [9, 12] with a few modifications to take the additional normalization steps in Algorithm 1 into account. First, we recall that

[TABLE]

where $\bm{\vartheta}_{t}=[\bm{w}_{t}^{\top},\bm{w}_{t-1}^{\top}]^{\top}$, $\bm{\xi}_{t}=[\bm{e}_{t}^{\top},\mathbf{0}]^{\top}$, and $\gamma_{t+1}$ is a normalization term. Let

[TABLE]

then $\left\|{\bm{v}_{1}}\right\|_{\mathbf{\mathcal{B}}}=1$ since $\left\|{\bm{u}_{1}}\right\|_{\mathbf{B}}=1$. Note that $\mathbf{\mathcal{B}}$ is positive definite, hence the $\mathbf{\mathcal{B}}$-norm is well defined. Let $\mathrm{P}_{\bm{v}_{1}}$ be a projection onto the space spanned by the top eigenvector of $\mathbf{C}$ and let $\mathrm{P}^{\perp}_{\bm{v}_{1}}$ be a projection onto the space spanned by the remaining ones, with respect to the $\mathbf{\mathcal{B}}$-norm. It can be verified that $\mathrm{P}_{\bm{v}_{1}}=\bm{v}_{1}\bm{v}_{1}^{\top}\mathbf{\mathcal{B}}$ and $\mathrm{P}^{\perp}_{\bm{v}_{1}}=\mathbf{I}-\bm{v}_{1}\bm{v}_{1}^{\top}\mathbf{\mathcal{B}}$. By the definition of principal angles under the $\mathcal{B}$-norm for $k=1$, and the facts that $\mathrm{P}_{\bm{v}_{1}}=\mathrm{P}_{\bm{v}_{1}}^{2}$ and $\mathrm{P}^{\perp}_{\bm{v}_{1}}=(\mathrm{P}^{\perp}_{\bm{v}_{1}})^{2}$, it is readily verified that $\cos\theta\left(\bm{\vartheta}_{t},\bm{v}_{1}\right)=\frac{\left\|{\mathrm{P}_{\bm{v}_{1}}\bm{\vartheta}_{t}}\right\|_{\mathbf{\mathcal{B}}}}{\left\|{\bm{\vartheta}_{t}}\right\|_{\mathbf{\mathcal{B}}}}$ and $\sin\theta\left(\bm{\vartheta}_{t},\bm{v}_{1}\right)=\frac{\left\|{\mathrm{P}^{\perp}_{\bm{v}_{1}}\bm{\vartheta}_{t}}\right\|_{\mathbf{\mathcal{B}}}}{\left\|{\bm{\vartheta}_{t}}\right\|_{\mathbf{\mathcal{B}}}}$.

We now start by expanding $\tan\theta\left(\bm{\vartheta}_{t+1},\bm{v}_{1}\right)$ as follows

[TABLE]

Since

[TABLE]

it follows that

[TABLE]

To get a geometric rate of convergence, it suffices to choose the noise level such that the rightmost term in (A.2) is less than 1. Since $\left\|{\bm{\vartheta}_{t}}\right\|_{\mathbf{\mathcal{B}}}=\sqrt{1+\gamma_{t}^{-2}}>1$, choosing

[TABLE]

yields

[TABLE]

completing the proof of Lemma 3.

Appendix B Proof of Theorem 1

First, by applying inequality (16) recursively, we obtain

[TABLE]

where the last step follows since $\left|{\mu_{1}}\right|>\left|{\mu_{2}}\right|$. We next examine the relative eigenvalue gap of the extended matrix $\mathbf{C}$. For a fixed $i$, from Lemma 2, the eigenvalues $\mu_{i}$ of $\mathbf{C}$ are the solutions of the characteristic equation

[TABLE]

If $2\sqrt{\beta}>\left|{\lambda_{1}}\right|$, the roots of the characteristic equation are imaginary, and both have magnitude $2\sqrt{\beta}$, which results in a zero eigenvalue gap. On the other hand, if $2\sqrt{\beta}<\left|{\lambda_{2}}\right|$, then

[TABLE]

while, for $\left|{\lambda_{2}}\right|\leq 2\sqrt{\beta}<\left|{\lambda_{1}}\right|$, we get

[TABLE]

We focus on the last case, since this is the range for $\beta$ which can offer accleration. Note that the right-hand side of (18) is a decreasing function of $\beta$, hence achieving its optimal solution at $\beta=\left|{\lambda_{2}}\right|^{2}/4$. Combining (B) and (18) yields

[TABLE]

We now pay attention to the running time of Algorithm 1 required to achieve an $\epsilon$-optimal solution, when the parameter $\beta$ is chosen optimally. For $\beta=\left|{\lambda_{2}}\right|^{2}/4$, we have

[TABLE]

From (B), it suffices to choose $k$ such that

[TABLE]

We now turn to control the errors arising from solving the least squares problems inexactly. Define the error in solving the least squares problem as

[TABLE]

Let $\textsf{r}_{\mathrm{init}}=\textsf{r}\left(\alpha_{t}\bm{w}_{t}\right)$ and $\textsf{r}_{\mathrm{des}}=\left\|{\bm{e}_{t}}\right\|_{\mathbf{B}}^{2}=\textsf{r}\left(\tilde{\bm{w}}_{t+1}\right)$ be the initial and required residual errors in Step 4 of Algorithm 1, respectively. Then, with our choice of initialization, we have by [9] that $\textsf{r}_{\mathrm{init}}\leq\lambda_{1}^{2}\sin^{2}\theta\left(\bm{w}_{t},\bm{u}_{1}\right)$. Therefore, to guarantee the noise condition in (10), it suffices to enforce that

[TABLE]

Thus, we only need to solve the least squares problem until the ratio of the final to initial error satisfies

[TABLE]

For $\left|{\lambda_{2}}\right|\leq 2\sqrt{\beta}<\left|{\lambda_{1}}\right|$, we have $\left|{\mu_{1}}\right|\geq\left|{\lambda_{1}}\right|$, hence if the errors above satisfy

[TABLE]

the condition (22) is automatically satisfied.

We can decompose the running time into two phases: i) initial phase due to large initial angle; and ii) convergence phase due to high accuracy $\epsilon$ required. In the initial phase, the angle are large, the min in (B) is determined by the second term. Thus, we can set the ratio to be $\frac{\Delta\cos^{2}\theta\left(\bm{\vartheta}_{0},\bm{v}_{1}\right)}{32}$, until $\tan\theta\left(\bm{\vartheta}_{t},\bm{v}_{1}\right)$ reduces to 1. This phase takes $\mathcal{O}\left(\frac{\log\tan\theta\left(\bm{\vartheta}_{0},\bm{v}_{1}\right)}{\sqrt{\Delta}}\right)$ iterations corresponding to $\epsilon=1$ in (21). In the convergence phase, the angles are small, the first term determines the min in (B), we can set the ratio to be $\frac{\Delta}{32}$ until we reach the target accuracy. This phase takes $\mathcal{O}\left(\frac{1}{\sqrt{\Delta}}\log\frac{1}{\epsilon}\right)$ iterations corresponding to $\tan\theta\left(\bm{\vartheta}_{0},\bm{v}_{1}\right)=1$ in (B).

Therefore, if we let $\mathcal{T}\left(\tau\right)$ be the time that a least squares solver takes to reduce the residual error by a factor $\tau$, then the total running time of Step 4 in Algorithm 1 is given by

[TABLE]

Finally, the running time for Steps 3 and 6 is

[TABLE]

which completes the proof of Theorem 1.

Appendix C Proof of Theorem 2

To begin with, recall that $\mathbf{U}$ is the matrix of top-$k$ eigenvectors of $\mathbf{B}^{-1}\mathbf{A}$ and $\bar{\mathbf{V}}=[\mathbf{\Lambda}\mathbf{U}^{\top},\mathbf{U}^{\top}]^{\top}$ is the top-$k$ eigenvectors of the extended matrix $\mathbf{C}$ with $\mathbf{\Lambda}=\mathrm{diag}\left(\mu_{1},\ldots,\mu_{k}\right)$ being the corresponding matrix of eigenvalues. Let

[TABLE]

then $\mathbf{V}^{\top}\mathbf{\mathcal{B}}\,\mathbf{V}=\mathbf{I}_{k}$. If we further let $\mathbf{V}_{\perp}$ be an orthogonal basis w.r.t $\mathbf{\mathcal{B}}$ of the orthogonal complement of $\rm{span}(\mathbf{V})$, then one can decompose $\mathbf{C}$ as

[TABLE]

where $\mathbf{\Lambda}_{\perp}=\mathrm{diag}\left(\mu_{k+1},\ldots,\mu_{2d-k}\right)$. Let $\bm{\Xi}_{t}=[\mathbf{E}_{t}^{\top},\mathbf{0}]^{\top}$, then the iterates generated by Algorithm 2 can be expressed as

[TABLE]

where $\mathbf{Y}_{t}=[\mathbf{W}_{t}^{\top},\mathbf{W}_{t-1}^{\top}]^{\top}$ with $\mathbf{W}_{t}^{\top}\mathbf{B}\mathbf{W}_{t}=\mathbf{I}$ and $\mathbf{W}_{t-1}^{\top}\mathbf{B}\mathbf{W}_{t-1}=\mathbf{R}_{t}^{-T}\mathbf{R}_{t}^{-1}$. Therefore, if we let $\bar{\mathbf{Y}}_{t}=\mathbf{Y}_{t}(\mathbf{I}+\mathbf{R}_{t}^{-T}\mathbf{R}_{t}^{-1})^{-1/2}$, then $\bar{\mathbf{Y}}_{t}^{\top}\mathbf{\mathcal{B}}\,\bar{\mathbf{Y}}_{t}=\mathbf{I}.$

By Proposition 1, \tan\theta\big{(}\rm{span}(\mathbf{Y}_{t+1}),\rm{span}(\bar{\mathbf{V}})\big{)} w.r.t $\mathbf{\mathcal{B}}$ is given by

[TABLE]

which implies that

[TABLE]

where the third equality follows from (24) and the facts that $\mathbf{V}^{\top}\mathbf{\mathcal{B}}\,\mathbf{V}=\mathbf{I}$ and $\mathbf{V}_{\perp}^{\top}\mathbf{\mathcal{B}}\,\mathbf{V}_{\perp}=\mathbf{I}$. Since

[TABLE]

where $\sigma_{k}(\cdot)$ denotes the $k$th largest singular value, it follows that

[TABLE]

We have

[TABLE]

where the last step follows from the definition of $\cos\theta(\bar{\mathbf{Y}}_{t},\mathbf{V})$. Substituting (28) into (27) yields

[TABLE]

where we have used the fact that $\left\|{\mathbf{V}_{\perp}^{\top}\mathbf{\mathcal{B}}\,\bm{\Xi}_{t}}\right\|\leq\left\|{\bm{\Xi}_{t}}\right\|_{\mathbf{\mathcal{B}}}=\left\|{\mathbf{E}_{t}}\right\|_{\mathbf{B}}$. As a consequence, if $\mathbf{E}_{t}$ satisfies

[TABLE]

then we obtain

[TABLE]

Having established the inequality (31), we can now use similar steps as in the proof of Theorem 1 to obtain the following iteration complexity:

[TABLE]

where $\Delta_{k}=1-\lambda_{k+1}/\lambda_{k}$.

Similar to the proof of Theorem 1, we define the error in solving the least squares problem as

[TABLE]

Let $\textsf{r}_{\mathrm{init}}=\textsf{r}\left(\mathbf{W}_{t}\mathbf{Z}_{t}\right)$ and $\textsf{r}_{\mathrm{des}}=\left\|{\mathbf{E}_{t}}\right\|_{\mathbf{B}}^{2}=\textsf{r}\left(\tilde{\mathbf{W}}_{t+1}\right)$, then it has been shown in [9] that $\textsf{r}_{\mathrm{init}}\leq 4k\left|{\lambda_{1}}\right|^{2}\tan\theta(\mathbf{W}_{t},\mathbf{U})$. Combining with (30), it is thus sufficient to decrease the error until

[TABLE]

For $\left|{\lambda_{k+1}}\right|\leq 2\sqrt{\beta}<\left|{\lambda_{k}}\right|$, we have $\left|{\mu_{k}}\right|\geq\left|{\lambda_{k}}\right|$, hence if the errors above satisfy

[TABLE]

where $\gamma=\left|{\lambda_{1}}\right|/\left|{\lambda_{k}}\right|$, then the condition (30) will be satisfied. Finally, using the same two-phase analysis as in the proof of Theorem 1, we arrive at the desired running time.