Global Convergence of a Grassmannian Gradient Descent Algorithm for Subspace Estimation

Introduction

Low-rank matrix factorization is one of the foundational tools of signal processing, numerical methods, and data analysis. Suppose we wish to factorize a matrix $M=UW^{T}$, imposing orthogonality constraints on $U$ or $W$. Solving for such matrix factorizations can be computationally burdensome, and many algorithms that attempt to speed up computation are actually solving a non-convex optimization problem, therefore coming with few guarantees.

Incremental gradient descent is our focus, motivated by streaming data applications. There are many applications of subspace estimation and tracking in medical imaging, communications, and environmental science; see more in [edelman1998geometry, balzano2014local, balzano2012handling]. Matrix factors with orthogonality constraints, such as those given by the SVD, are also used in several data applications: they provide a unique collection of low-dimensional projections for data visualization, capture directions of maximal variance so as to give useful insights into data structure, and allow compressed storage of massive datasets with a precise notion of loss in compression.

Formulation and Related Work

i.e., the span of the data vectors or the range of $\bar{U}$, denoted $R(\bar{U})$. When $\xi_{t}\neq 0$ we still discuss results in terms of the distance from $\bar{U}$. If we consider only $t=1,\dots,N$, Problem (2) is identical to Problem \eqrefprob:batch. The GROUSE algorithm (Grassmannian Rank-One Update Subspace Estimation) we analyze is shown as Algorithm 1, where we generate a sequence $\{U_{t}\}_{t=0,1,\dots}$ of $n\times d$ matrices with orthonormal columns with the goal that $R(U_{t})\rightarrow R(\bar{U})$ as $t\rightarrow\infty$. Each observed vector is used to update $U_{t}$ to $U_{t+1}$, and we constrain the gradient descent method to the Grassmannian using a geodesic update [edelman1998geometry].

Because of the importance of the problem, it has been studied for decades, and there is a great deal of related work. We direct the reader to [edelman1998geometry, balzano2012handling] for in-depth descriptions of algorithms and guarantees. We focus here on recent results that have global convergence guarantees to the global minimizer and study either gradient-type algorithms, algorithms that handle streaming data, or algorithms that maintain orthogonality constraints with manifold optimization.

First we discuss incremental methods. [de2014global] established the global convergence of a stochastic gradient descent method for the recovery of a positive definite matrix $M$ in the undersampled case, where the matrix $M$ is not measured directly but instead via linear measurements. They propose a step size scheme under which they prove global convergence results from a randomly generated initialization. Similarly, [balsubramani2013fast] invokes a martingale-based argument to show the global convergence rate of the proposed incremental PCA method to the single top eigenvector in the fully sampled case. In contrast, [arora2013stochastic] estimates the best $d$-dimensional subspace in the fully sampled case and provides a global convergence result by relaxing the non-convex problem to a convex one. We seek to identify the $d$ dimensional subspace by solving the non-convex problem directly. Finally, our work is most related to [balzano2014local], which provides local convergence guarantees for GROUSE in both the fully sampled and undersampled case. Our work focuses on global convergence but only in the fully sampled case; we will extend the global convergence results to the undersampled case in future work.

Turning to batch methods, [RH2012, jain2013low] provided the first theoretical guarantee for an alternating minimization algorithm for low-rank matrix recovery in the undersampled case. Under typical assumptions required for the matrix recovery problems [recht2010guaranteed], they established geometric convergence to the global optimal solution. Earlier work [keshavan2010matrix, ngo2012scaled] considered the same undersampled problem formulation and established convergence guarantees for a steepest descent method (and a preconditioned version) on the full gradient, performed on the Grassmannian. [chen2015fast, bhojanapalli2015dropping, zheng2015convergent] considered low rank semidefinite matrix estimation problems, where they reparamterized the underlying matrix as $M=UU^{T}$, and update $U$ via a first order gradient descent method. However, all these results require batch processing and a decent initialization that is close enough to the optimal point, resulting in a heavy computational burden and precluding problems with streaming data. We study random initialization, and our algorithm has fast, computationally efficient updates that can be performed in an online context.

Lastly, several convergence results for optimization on general Riemannian manifolds, including several special cases for the Grassmannian, can be found in [absil2009optimization]. Most of the results are very general; they include global convergence rates to local optima for steepest descent, conjugate gradient, and trust region methods, to name a few. We instead focus on solving the problem in \eqrefeq:cost and provide global convergence rates to the global minimum.

Convergence analysis

Before we present our main results on the convergence of the GROUSE algorithm, we first call out the following definitions and condition that will be used throughout our analysis.

We use $\phi_{i}\left(\bar{U},U_{t}\right),i=1,\dots,d$ to denote the principal angles between subspaces $\text{R}(U_{t})$ and $\text{R}(\bar{U})$, which are defined [[stewart1990matrix], Chapter 5] by $\cos\phi_{i}(\bar{U},U_{t})=\sigma_{i}(\bar{U}^{T}U_{t})$.

Our first metric is $\zeta_{t}\in$, which measures the similarity between two subspaces and is defined as

Our second metric is $\epsilon_{t}\in[0,d]$, which measures the discrepancy between $R(U_{t})$ and $R(\bar{U})$, and is defined as

In this section, we first derive a greedy step size scheme for each iteration $t$ that maximizes the improvement on the defined metrics ($\epsilon_{t},\zeta_{t}$) of convergence for the noiseless case, i.e., $x_{t}=v_{t}$. Let $v_{t,\parallel}$ and $v_{t,\perp}$ denote the projection and residual of $v_{t}$ onto $R(U_{t})$. Then after each update we have the following (Appendix LABEL:sec:proof_of_supporting_theory):