# On Approximation Guarantees for Greedy Low Rank Optimization

We provide new approximation guarantees for greedy low rank matrix estimation under standard assumptions of restricted strong convexity and smoothness. Our novel analysis also uncovers previously unknown connections between the low rank estimation and combinatorial optimization, so much so that our bounds are reminiscent of corresponding approximation bounds in submodular maximization. Additionally, we also provide statistical recovery guarantees. Finally, we present empirical comparison of greedy estimation with established baselines on two important real-world problems.

## Authors

• 13 publications
• 3 publications
• 43 publications
• 10 publications
• ### Improved Algorithms for Matrix Recovery from Rank-One Projections

We consider the problem of estimation of a low-rank matrix from a limite...
05/21/2017 ∙ by Mohammadreza Soltani, et al. ∙ 0

• ### Regularized Weighted Low Rank Approximation

The classical low rank approximation problem is to find a rank k matrix ...
11/16/2019 ∙ by Frank Ban, et al. ∙ 0

• ### Learning Non-Parametric Basis Independent Models from Point Queries via Low-Rank Methods

We consider the problem of learning multi-ridge functions of the form f(...
10/07/2013 ∙ by Hemant Tyagi, et al. ∙ 0

• ### Global Optimality in Distributed Low-rank Matrix Factorization

We study the convergence of a variant of distributed gradient descent (D...
11/07/2018 ∙ by Zhihui Zhu, et al. ∙ 0

• ### Sketchy Decisions: Convex Low-Rank Matrix Optimization with Optimal Storage

This paper concerns a fundamental class of convex matrix optimization pr...
02/22/2017 ∙ by Alp Yurtsever, et al. ∙ 0

• ### An improved analysis and unified perspective on deterministic and randomized low rank matrix approximations

We introduce a Generalized LU-Factorization (GLU) for low-rank matrix ap...
10/01/2019 ∙ by James Demmel, et al. ∙ 0

• ### Robust High-Dimensional Linear Regression

The effectiveness of supervised learning techniques has made them ubiqui...
08/07/2016 ∙ by Chang Liu, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Low rank matrix estimation stands as a major tool in modern dimensionality reduction and unsupervised learning. The singular value decomposition can be used when the optimization objective is rotationally invariant to the parameters. However, if we wish to optimize over more complex objectives we must choose to either optimize over the non-convex space (which have seen recent theoretical success)

[1, 2, 3, 4, 5] or rely on convex relaxations to the non-convex optimization [6, 7, 8].

More concretely, in the low-rank matrix optimization problem we wish to solve

 argmaxΘℓ(Θ)s.t. rank(Θ)≤r, (1)

and rather than perform the computationally intractable optimization above researchers have studied convex relaxations of the form

 argmaxΘℓ(Θ)−λ|||Θ|||nuc.

Unfortunately, the above optimization can be computationally taxing. General purpose solvers for the above optimization problem that rely on semidefinite programming require computation, which is prohibitive. Gradient descent techniques require computational cost for an epsilon accurate solution. This improvement is sizeable in comparison to SDP solvers. Unfortunately, it is still prohibitive for large scale matrix estimation.

To alleviate some of the computational issues an alternate vein of research has focused on directly optimizing the non-convex problem in (1). To that end, authors have studied the convergence properties of

 argmaxU∈Rn×r,V∈Rd×rℓ(UVT).

Solving the problem above automatically forces the solution to be low rank and recent results have shown promising behavior. An alternative approach is to optimize via rank one updates to the current estimate [9, 10]. This approach has also been studied in more general contexts such as boosting [11], coordinate descent [12, 13], and incremental atomic norm optimization [14, 15, 16, 17].

### 1.1 Set Function Optimization and Coordinate Descent

The perspective that we take is treating low rank matrix estimation as a set optimization over an infinite set of atoms. Specifically, we wish to optimize

 argmax{X1,…Xk}∈Aℓ(k∑i=1αiXi),

where the set of atoms is the set of all rank one matrices with unit operator norm. This settings is analogous to that taken in the results studying atomic norm optimization, coordinate descent via the norm in total variation, and Frank-Wolfe style algorithms for atomic optimization. This formulation allows us to connect the problem of low rank matrix estimation to that of submodular set function optimization, which we discuss in the sequel. Before proceeding we discuss related work and an informal statement of our result.

### 1.2 Informal Result and Related Work

Our result demonstrates an exponential decrease in the amount of error suffered by greedily adding rank one matrices to the low rank matrix approximation.

###### Theorem 1 (Approximation Guarantee, Informal).

If we let be our estimate of the rank matrix at iteration , then for some universal constant related to the restricted condition number of the problem we have

 ℓ(Θk)−ℓ(0)≥(1−exp(−ck/r))(ℓ(Θ∗)−ℓ(0)).

Note that after iterations the matrix will be at most rank . Now, we can contrast this result to related work.

##### Related work:

There has been a wide array of studies looking at the computational and statistical benefits of rank one updates to estimating a low rank matrix. At its most basic, the singular value decomposition will keep adding rank one approximations through deflation steps. Below we discuss a few of the results.

The work can be generally segmented into to sets of results. Those results that present sublinear rates of convergence and those that obtain linear rates. Interestingly, parallel lines of work have also demonstrated similar convergence bounds for more general atomic or dictionary element approximations [11, 14, 15, 16]. For space constraints, we will summarize these results into two categories rather than explicitly state the results for each individual paper.

If we define the atomic norm of a matrix to be to be the sum of the singular values of that matrix, then the bounds establish in the sublinear convergence cases behave as

 ℓ(Θ∗)−ℓ(Θk)≤|||Θ∗|||2nuck,

where we take to be the best rank solution. What we then see is convergence towards the optimal bound. However, we expect our statistical error to behave as where is the number of samples that we have received from our statistical model and is rank  [7, 8]. We can take , which would then imply that we would need to behave as . However, that would then imply that the rank of our matrix should grow linearly in the number of observations in order to achieve the same statistical error bounds. The above error bounds are “fast.” If we consider a model that yields slow error bounds, then we expect the error to behave like . In that case, we can take , which looks better, but still requires significant growth in as a function of .

To overcome the above points, some authors have aimed to study similar greedy algorithms that then enjoy exponential rates of convergence as we show in our paper. These results share the most similarities with our own and behave as

 ℓ(Θk)≥(1−γk)ℓ(Θ∗)

where is the best over all set of parameters. This result decays exponentially. However, when one looks at the behavior of it will typically act as , for an matrix. As a result, we would need to take on the order of the number of the dimensionality of the problem in order to begin to see gains. In contrast, for our result listed above, if we seek to only compare to the best rank solution, then the gamma we find is . Of course, if we wish to find a solution with full-rank, then the bounds we stated above match the existing bounds.

In order to establish our results we rely on a notion introduced in the statistical community called restricted strong convexity. This assumption has connections to ideas such as the Restricted Isometry Property, Restricted Eigenvalue Condition, and Incoherence. In the work by Shalev-Shwartz, Gonen, and Shamir

[9] they present results under a form of strong convexity condition imposed over matrices. Under that setting, the authors demonstrate that

 ℓ(Θk)≥ℓ(Θ∗)−ℓ(0)rk

where is the rank of . In contrast, our bound behaves as

 ℓ(Θk)≥ℓ(Θ∗)+(ℓ(Θ∗)−ℓ(0))exp(\nicefrac−kr)
##### Our contributions:

We improve upon the linear rates of convergence for low-rank approximation using rank one updates by connecting the coordinate descent problem to that of submodular optimization. We present this result in the sequel along with the algorithmic consequences. We demonstrate the good performance of these rank one updates in the experimental section.

## 2 Background

We begin by fixing some notation. We represent sets using sans script fonts e.g.

. Vectors are represented using lower case bold letters

e.g. , and matrices are represented using upper case bold letters e.g. . Non-bold face letters are used for scalars e.g. and function names e.g. . The transpose of a vector or a matrix is represented by e.g. . Define . For singleton sets, we write . Size of a set is denoted by . is used for matrix inner product.

Our goal is to analyze greedy algorithms for low rank estimation. Consider the classic greedy algorithm that picks up the next element myopically i.e. given the solution set built so far, the algorithm picks the next element as the one which maximizes the gain obtained by adding the said element into the solution set. Approximation guarantees for the greedy algorithm readily imply for the subclass of functions called submodular functions which we define next.

###### Definition 1.

A set function is submodular if for all ,

 f(A)+f(B)≥f(A∪B)+f(A∩B).

Submodular set functions are well studied and have many desirable properties that allow for efficient minimization, and maximization with approximation guarantees. Our low rank estimation problem also falls under the purview of another class of functions called monotone functions. A function is called monotone if and only if for all . For the specific case of maximizing monotone submodular set functions, it is known that the greedy algorithm run for (say) iterations is guaranteed to return a solution that is within of the optimum set of size  [18]. Moreover, without further assumptions or knowledge of the function, no other polynomial time algorithm can provide a better approximation guarantee unless P=NP [19].

More recently, a line of works have shown that the greedy approximation guarantee that is typically applicable to monotone submodular functions can be extended to a larger class of functions called weakly submodular functions [20, 21]. Central to the notion of weak submodularity is the quantity submodularity ratio which we define next.

###### Definition 2 (Submodularity Ratio [22]).

Let be two disjoint sets, and . The submodularity ratio of with respect to is given by

 γL,S:=∑j∈S[f(L∪{j})−f(L)]f(L∪S)−f(L). (2)

The submodularity ratio of a set with respect to an integer is given by

 (3)

It is easy to show that is submodular if and only if for all sets and . However, as noted by Das and Kempe [22], Elenberg et al. [20], an approximation guarantee is guaranteed when , thereby extending the applicability of the greedy algorithm to a much larger class of functions. The subset of monotone functions which have are called weakly submodular functions in the sense that even though the function is not submodular, it still provides provable bound for greedy selections.

Also vital to our analysis is the notion of restricted strong concavity and smoothness [23, 24].

###### Definition 3 (Low Rank Restricted Strong Concavity, Restricted Smoothness).

A function is said to be restricted strong concave with parameter and restricted smooth with parameter if for all ,

 −mΩ2∥Y−X∥2F ≥ℓ(Y)−ℓ(X)−⟨∇ℓ(X),Y−X⟩ ≥−MΩ2∥Y−X∥2F.
###### Remark 1.

If a function has restricted strong concavity parameter , then its negative has restricted strong convexity parameter . We choose to use the nomenclature of concavity for ease of exposition in terms of relationship to submodular maximization. Further, note that we define RSC/RSM conditions on the space of matrices rather than vectors.

It is straightforward to see that if , then and .

## 3 Setup

In this section, we delineate our setup of low rank estimation. For the sake of convenience of relating to the framework of weak submodular maximization, we operate in the setting of maximization of a concave matrix variate function under a low rank constraint. This is equivalent to minimizing a convex matrix variate function under the low rank constraint as considered by Shalev-Shwartz et al. [9] or under nuclear norm constraint or regularization as considered by Jaggi and Sulovský [13]. The goal is to maximize a function under a low rank constraint:

 maxrank(X)≤rℓ(X). (4)

Instead of using a convex relaxation of the constrained problem (4), our approach is to enforce the rank constraint directly by adding rank matrices greedily until is of rank . The rank matrices to be added are obtained as outer product of vectors from the given vector sets and . e.g. and .

The problem (4) can be interpreted in a context of sparsity as long as and are enumerable. For example, by using the SVD theorem, it is known that we can rewrite as , where , and . By using enumeration of the sets and under a finite precision representation of real values, one can rethink of the optimization (4) as finding a sparse solution for the infinite dimensional vector  [9]. As a first step, we can define an optimization over specified support sets, similar to choosing support for classical sparsity in vectors. For a support set , let and be the matrices formed by stacking the chosen elements of and respectively. For the support , we can define a set function that maximizes over :

 f(L)=maxH∈R|L|×|L|ℓ(U⊤LHVL)−ℓ(0). (5)

We will denote the optimizing matrix for a support set as . In other words, let be the argmax obtained in (5), then .

Thus, the low rank matrix estimation problem (4) can be reinterpreted as the following equivalent combinatorial optimization problem:

 max|S|≤kf(S). (6)

### 3.1 Algorithms

We briefly state the algorithms. Our greedy algorithm is illustrated in Algorithm 1. The greedy algorithm builds the support set incrementally – adding a rank matrix one at a time, so that at iteration for the size of the chosen support set and hence rank of the current iterate is . We assume access to a subroutine GreedySel for the greedy selection (Step 4). This subroutine solves an inner optimization problem by calling a subroutine GreedySel which returns an atom from the candidate support set that ensures

 f(SGi−1∪{s})−f(SGi−1)≥τ(f(SGi−1∪{s⋆})−f(SGi−1)),

where

 s⋆←argmaxa∈(U×V)⊥SGi−1f(SGi−1∪{a})−f(SGi−1).

In words, the subroutine GreedySel ensures that the gain in obtained by using the selected atom is within multiplicative approximation to the atom with the best possible gain in

. The hyperparameter

governs a tradeoff allowing a compromise in myopic gain for a possibly quicker selection.

The greedy selection requires to fit and score every candidate support, which is prohibitively expensive. An alternative is to choose the next atom by using the linear maximization oracle used by Frank-Wolfe [12] or Matching Pursuit algorithms [14]. This step replaces Step 4 of Algorithm 1 as illustrated in Algorithm 2. Let be the set constructed by the algorithm at iteration . The linear oracle OMPSel returns an atom for iteration ensuring

 ⟨∇ℓ(B(L)),usv⊤s⟩≥τmax(u,v)∈(U×V)⊥SOi−1⟨∇ℓ(B(L)),uv⊤⟩.

The linear problem OMPSel can be considerably faster that GreedySel. OMPSel  reduces to finding the left and right singular vectors of corresponding to its largest singular value, which is , where is the number of non-zero entries of .

Algorithm 2 is the same as considered by Shalev-Shwartz et al. [9] as GECO (Greedy Efficient Component Optimization). However, as we shall see, our analysis provides stronger bounds than their Theorem 2.

###### Remark 2.

We note that Step 5 of Algorithms 1,2 requires solving the RHS of (5) which is a matrix variate problem of size at iteration . This refitting is equivalent to the “fully-corrective” versions of Frank-Wolfe/Matching Pursuit algorithms which, intuitively speaking, extract out all the information w.r.t from the chosen set of atoms, thereby ensuring that the next rank atom chosen has row and column space orthogonal to the previously chosen atoms. Thus the constrained maximization on the orthogonal complement of in subroutines OMPSel and GreedySel need not be explicitly enforced, but is still shown for clarity.

## 4 Analysis

In this section, we prove that low rank matrix optimization over the rank one atoms satisfies weak submodularity. This helps us bound the function value obtained till greedy iterations vis-a-vis the function value at the optimal sized selection.

We explicitly delineate some notation and assumptions. With slight abuse of notation, we assume is -strongly concave and -smooth over matrices of rank . For , note that and . Additionally, let , and assume is -smooth over . It is easy to see .

Since we obtain approximation bounds for the greedy algorithm similar to ones obtained using classical methods, we must also mention the corresponding assumptions in the submodular literature that further draws parallels to our analysis. Submodularity guarantees that greedy maximization of monotone normalized functions yields a approximation. Since we are doing support selection, increasing the support size does not decrease the function value. Hence the set function we consider is monotone. Further, we also subtract to make sure (see (5)) so that our set function is also normalized. We shall see that our bounds are of similar flavor as the classical submodularity bound.

As the first step, we prove that if the low rank RSC holds, then the submodularity ratio (Definition 2) is lower-bounded by the inverse condition number.

###### Theorem 2.

Let be a set of rank atoms and be a set of rank atoms where we sequentially orthogonalize the atoms against . If is -strongly concave over matrices of rank , and -smooth over the set , then

 γL,r:=∑a∈S[f(L∪{a})−f(L)]f(L∪S)−f(L)≥mr+k~M1.

The proof of Theorem 2 is structured around individually obtaining a lower bound for the numerator and an upper bound for the denominator of the submodularity ratio by exploiting the concavity and convexity conditions. Bounding the submodularity ratio is crucial to obtaining the approximation bounds for the greedy algorithm as we shall see in the sequel.

### 4.1 Greedy Improvement

In this section, we obtain approximation guarantees for Algorithm 1.

###### Theorem 3.

Let be the greedy solution set obtained by running Algorithm 1 for iterations, and let be an optimal support set of size . Let be strongly concave on the set of matrices with rank less than or equal to , and smooth on the set of matrices in the set . Then,

 f(S) ≥(1−1ec1)f(S⋆) ≥(1−1ec2)f(S⋆),

where and .

The proof technique for the first inequality of Theorem 3 relies on lower bounding the progress made in each iteration of Algorithm 1. Intuitively, it exploits the weak submodularity to make sure that each iteration makes enough progress, and then applying an induction argument for iterations. We also emphasize that the bounds in Theorem 3 are for normalized set function (which means ). A more detailed proof is presented in the appendix.

###### Remark 3.

Theorem 3 provides the approximation guarantees for running the greedy selection algorithm up to iterations to obtain a rank matrix iterate vis-a-vis the best rank approximation. For , and , we get an approximation bound which is reminiscent of the greedy bound of under the framework of submodularity. Note that our analysis can not be used to establish classical submodularity. However, establishing weak submodularity that lower bounds is sufficient to provide slightly weaker than classical submodularity guarantees.

###### Remark 4.

Theorem 3 implies that to obtain approximation guarantee in the worst case, running Algorithm 1 for iterations suffices. This is useful when the application allows a tradeoff: compromising on the low rank constraint a little to achieve tighter approximation guarantees.

###### Remark 5.

Das and Kempe [22] considered the special case of greedily maximizing

statistic for linear regression, which corresponds to classical sparsity in vectors. They also obtain a bound of

, where is the submodularity ratio for their respective setup. This was generalized by Elenberg et al. [20] to general concave functions under sparsity constraints. Our analysis is for the low rank constraint, as opposed to sparsity in vectors that was considered by them.

### 4.2 GECO Improvement

In this section, we obtain the approximation guarantees for Algorithm 2. The greedy search over the infinitely many candidate atoms is infeasible, especially when . Thus while Algorithm 1 establishes interesting theoretical connections with submodularity, it is, in general, not practical. To obtain a tractable and practically useful algorithm, the greedy search is replaced by a Frank Wolfe or Matching Pursuit style linear optimization which can be easily implemented as finding the top singular vectors of the gradient at iteration . In this section, we show that despite the speedup, we lose very little in terms of approximation guarantees. In fact, if the approximation factor in OMPSel() is , we get the same bounds as those obtained for the greedy algorithm.

We now present our main result for Algorithm 2.

###### Theorem 4.

Let be the greedy solution set obtained using Algorithm 2 for iterations, and let be the optimum size support set. Let be strongly concave on the set of matrices with rank less than or equal to , and smooth on the set of matrices with rank in the set . Then,

 f(S) ≥(1−1ec3)f(S⋆),

where .

The proof of Theorem 4 follows along the lines of Theorem 3. The central idea is similar - to exploit the RSC conditions to make sure that each iteration makes sufficient progress, and then provide an induction argument for iterations. Unlike the greedy algorithm, however, using the weak submodularity is no longer required. Note that the bound obtained in Theorem 4 is similar to Theorem 3, except the exponent on the approximation factor .

###### Remark 6.

Our proof technique for Theorem 4 can be applied for classical sparsity to improve the bounds obtained by Elenberg et al. [20] for OMP for support selection under RSC, and by Das and Kempe [22] for statistic. If , their bounds involve terms of the form in the exponent, as opposed to our bounds which only has in the exponent.

## 5 Recovery Guarantees

While understanding approximation guarantees are useful, providing parameter recovery bounds can further help us understand the practical utility of the greedy algorithm. In this section, we present a general theorem that provides us with recovery bounds of the true underlying low rank structure.

###### Theorem 5.

Suppose that an algorithm achieves the approximation guarantee:

 f(Sk)≥Cr,kf(S⋆r),

where is the set of size at iteration of the algorithm, be the optimal solution for -cardinality constrained maximization of , and be the corresponding approximation ratio guaranteed by the algorithm. Recall that we represent by the matrices formed by stacking the vectors represented by the support set chosen from respectively, s.t. . Then under RSC, with for any , we have

 ∥B(Sk)−Br∥2F ≤4(k+r)∥∇ℓ(Br)∥22m2k+r +4(1−Cr,k)mk+r[ℓ(Br)−ℓ(0)]

Theorem 5 can be applied for , which is the argmax for maximizing under the low rank constraint. It is general - in the sense that it can be applied for getting recovery bounds from approximation guarantees for any algorithm, and hence is applicable for both Algorithms 1 and 2.

For specific function and statistical model, statistical recovery guarantees guarantees can be obtained from Theorem 5 for specific and statistical model, Consider the case of low rank matrix estimation from linear measurements. Say for are generated so that each entry of is . We observe , where is low rank, and say . Let , and let be the linear operator so that . Our corresponding function is now . For this function, using arguments by  Negahban et al. [23], we know and

with high probability. It is also straightforward to apply their results to bound

, and , which gives explicit bounds as per Theorem 5 for Algorithms 12 for the considered function and the design matrix.

## 6 Experiments

In this section, we empirically evaluate the proposed algorithms.

### 6.1 Clustering under Stochastic Block Model

In this section, we test empirically the performance of GECO (Algorithm 2) for a clustering task. We are provided with a graph with nodes and the respective edges between the nodes. The observed graph is assumed to have been noisily generated from a true underlying clustering. The goal is to recover the underlying clustering structure from the noisy graph provided to us. The adjacency matrix of the true underlying graph is low rank. As such, our greedy framework is applicable . We compare performance of Algorithm 2

on simulated data against standard baselines of spectral clustering which are commonly used for this task. We begin by describing a generative model for creating edges between nodes given the ground truth.

The Stochastic Block Model is a model to generate random graphs. It takes its input the set of nodes, and a partition of which form a set of disjoint clusters, and returns the graph with nodes and the generated edges. The model is also provided with generative probabilities – so that a pair of nodes within the same cluster have an edge between them with probability , while a pair of nodes belonging to different clusters have an edge between them with probability . For simplicity we assume . The model then iterates over each pair of nodes. For each such pair that belongs to same cluster, it samples an edge as Bernoulli(), otherwise as Bernoulli(). This provides us with a adjacency matrix.

For baselines, we compare against two versions of spectral clustering, which is a standard technique applied to find communities in a graph. The method takes as input the adjacency matrix , which is a matrix with an entry if there is an edge between node and , and is otherwise. From the adjacency matrix, the graph Laplacian is constructed. The Laplacian may be unnormalized, in which case it is simply , where is the diagonal matrix of degrees of nodes. A normalized Laplacian is computed as . After calculating the Laplacian, the algorithm solves for bottom eigenvectors of the Laplacian, and then apply -means clustering on the rows of the thus obtained eigenvector matrix. We refer to the works of Shi and Malik [25], Ng et al. [26] for the specific details of clustering algorithms using unnormalized and normalized graph Laplacian respectively.

We compare the spectral clustering algorithms with logistic PCA, which is a special case of the exponential family PCA [27]. The exponential family extension of classical PCA is analogous to the extension of the linear regression to generalized linear models (GLMs). For a given matrix , the GLM generative model assumes that each cell is independently drawn with likelihood proportional to , where is the true underlying parameter, and is the partition function. It is easy to see we can apply our framework of greedy selection by defining as the log-likelihood:

 ℓ(Θ)=⟨Θ,X⟩−∑i,jG(Θij),

where is the true parameter matrix of and that generates a realization of . Since the true is low rank, we get the low rank constrained optimization problem:

 maxrank(Θ)≤kℓ(Θ),

where is the hyperparameter that is suggestive of true number of clusters. Note that lack of knowledge of true value of is not more restrictive than spectral clustering algorithms which typically also require the true value of , albeit some subsequent works have tried to address tuning for .

Having cast the clustering problem in the same form as (4), we can apply our greedy selection algorithm as opposed to the more costly alternating minimizing algorithms suggested by Collins et al. [27]. Since the given matrix is with each entry sampled from a Bernoulli, we use which gives us logistic PCA.

We generate the data as follows. For nodes, and fixed number of cluster , we vary the within cluster edge generation probability from to in increments of , and use the Stochastic Block model to generate a noisy graph with each . Note that smaller means that the sampled graph will be more noisy and likely to be more different than the underlying clustering.

We compare against the spectral clustering algorithm using unnormalized Laplacian of Shi and Malik [25] which we label “Spectral_unnorm{k}” for , and the spectral clustering algorithm using normalized Laplacian of Ng et al. [26] which we label “Spectral_norm{k}” for . We use Algorithm 2 which we label “Greedy{k}” for . For each of these models, the referred is the supplied hyperparameter. We report the least squares error of the output from each model to the true underlying (which we call the generalization error), and to the instantiation used for training (which we call the reconstruction error). The results are presented in Figure 1.

Figure 1 shows that the greedy logistic PCA performs well in not only re-creating the given noisy matrix (reconstruction) but also captures the true low rank structure better (generalization). Further, note that providing the true hyper parameter is vital for spectral clustering algorithms, while on the other hand greedy is less sensitive to which is very useful in practice as is typically not known. So the spectral clustering algorithms typically would involve taking an SVD and re-running the for different values of to choose the best performing hyperparameter. The greedy factorization on the other hand is more robust, and moreover is incremental - it does not require to be re-run from scratch for different values of .

### 6.2 Word Embeddings

The task of embedding text into a vector space yields a representation that can have many advantages, such as using them as features for subsequent tasks as sentiment analysis

Mikolov et al. [28] proposed a context-based embedding called skip-gram or word2vec). The context of a word can be defined as a set of words before, around, or after the respective word. Their model strives to find an embedding of each word so that the representation predicts the embedding of each context word around it. In a recent paper, Levy and Goldberg [29] showed that the word embedding model proposed by Mikolov et al. [28] can be re-interpreted as matrix factorization of the PMI matrix constructed as follows. A word is in context of if it lies within the respective window of . The PMI matrix is then calculated as

 PMIw,c=log(p(w,c)p(w)p(c)).

In practice the probabilities are replaced by their empirical counterparts. Further, note that is if words and do not co-exist in the same, context which yields for PMI.  Levy and Goldberg [29] suggest using an alternative: . They also suggest variations of PMI hyper parameterized by which corresponds to the number of negative samples in the training of skip gram model of Mikolov et al. [28].

We employ the binomial model on the normalized count matrix (instead of the PMI), in a manner similar to the clustering approach in Section 6.1. The normalized counts matrix is calculated simply as , without taking explicit logs unlike the PMI  matrix. This gives us a probability matrix which has each entry between and , which can be factorized under the binomial model greedily as per Algorithm 2, similar to the way we do it in Section 6.1.

We note that embeddings using the SVD is more scalable than our greedy approach because of advancements in linear algebraic techniques for SVD on sparse matrices that PPMI yields. Our experiments show that binomial PCA can be competitive to other existing embedding methods. Since our current implementation is not as scalable, we are further investigating this as on-going work.

We empirically study the embeddings obtained by binomial factorization on two tasks - word similarity and analogies. For word similarity, we use the W353 dataset [30] which has 353 queries and the MEN data [31]

which has 3000 queries. Both these datasets contain words with human assigned similarity scores. We evaluate the embeddings by their cosine similarity, and measuring the correlation with the available human ratings. For the analogy task, we use the Microsoft Research (MSR) syntactic analogies

[32] which has 8000 queries, and the Google mixed analogies dataset [33] with 19544 queries. To compute accuracy, we use the multiplication similarity metric as used by Levy and Goldberg [29]. To train the word embeddings, we use the 2013 news crawl datasethttp://www.statmt.org/wmt14/training-monolingual-news-crawl. We filter out stop words and non-ascii characters, and keep only the words which occurs atleast 2000 times which yields vocabulary of 6713. Note that since we filter only the most common words, several queries from the datasets are invalid because we do not have embeddings for words appearing in them. However, we do include them and report the overall average over the entire dataset, with metric being by default for each query we are not able to process.

Table 1 shows the empirical evaluation. SVD and PPMI are the models proposed by Levy and Goldberg [29], while SGNS is skipgram with negative sampling model of Mikolov et al. [28]. We run each of these for and report the best numbers. The results show that alternative factorizations such as our application of binomial PCA to those of taking SVD of the PPMI  matrix can be more consistent and competitive with other embedding methods.

## 7 Conclusion

We have connected the problem of greedy low-rank matrix estimation to that of submodular optimization. Through that connection we have provided improved exponential rates of convergence for the algorithm. An interesting area of future study will be to connect these ideas to general atoms or dictionary elements.

## References

• Park et al. [2016] D. Park, A. Kyrillidis, S. Bhojanapalli, C. Caramanis, and S. Sanghavi, “Provable non-convex projected gradient descent for a class of constrained matrix optimization problems,” arXiv, vol. abs/1606.01316, 2016.
• Jain et al. [2013] P. Jain, P. Netrapalli, and S. Sanghavi, “Low-rank matrix completion using alternating minimization,” in

Symposium on Theory of Computing Conference, STOC’13, Palo Alto, CA, USA, June 1-4, 2013

, 2013, pp. 665–674.
• Chen and Wainwright [2015] Y. Chen and M. J. Wainwright, “Fast low-rank estimation by projected gradient descent: General statistical and algorithmic guarantees,” arXiv, vol. abs/1509.03025, 2015.
• Lee and Bresler [2013] K. Lee and Y. Bresler, “Corrections to ”admira: Atomic decomposition for minimum rank approximation”,” IEEE Trans. Information Theory, vol. 59, no. 7, pp. 4730–4732, 2013.
• Jain et al. [2014] P. Jain, A. Tewari, and P. Kar, “On iterative hard thresholding methods for high-dimensional m-estimation,” in Advances in Neural Information Processing Systems 27: Annual Conference on Neural Information Processing Systems 2014, December 8-13 2014, Montreal, Quebec, Canada, 2014, pp. 685–693.
• Recht et al. [2010] B. Recht, M. Fazel, and P. A. Parrilo, “Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization,” SIAM Review, vol. 52, no. 3, p. 471–501, Jan 2010.
• Negahban and Wainwright [2011] S. Negahban and M. J. Wainwright, “Estimation of (near) low-rank matrices with noise and high-dimensional scaling,” The Annals of Statistics, vol. 39, no. 2, pp. 1069–1097, 2011.
• Rohde and Tsybakov [2011] A. Rohde and A. B. Tsybakov, “Estimation of high-dimensional low-rank matrices,” The Annals of Statistics, vol. 39, no. 2, p. 887–930, Apr 2011.
• Shalev-Shwartz et al. [2011] S. Shalev-Shwartz, A. Gonen, and O. Shamir, “Large-scale convex minimization with a low-rank constraint,” in ICML, 2011.
• Wang et al. [2015] Z. Wang, M.-J. Lai, Z. Lu, W. Fan, H. Davulcu, and J. Ye, “Orthogonal rank-one matrix pursuit for low rank matrix completion,” SIAM Journal on Scientific Computing, vol. 37, no. 1, p. A488–A514, Jan 2015.
• Buhlmann and Yu [2009] P. Buhlmann and B. Yu, “Boosting,” Wiley Interdisciplinary Reviews: Computational Statistics, vol. 2, no. 1, p. 69–74, Dec 2009.
• Jaggi [2013] M. Jaggi, “Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization,” in ICML, 2013, pp. 427–435.
• Jaggi and Sulovský [2010] M. Jaggi and M. Sulovský, “A simple algorithm for nuclear norm regularized problems,” in

Proceedings of the 27th International Conference on Machine Learning (ICML-10)

, J. Fürnkranz and T. Joachims, Eds.   Omnipress, 2010, pp. 471–478.
• Gribonval and Vandergheynst [2006] R. Gribonval and P. Vandergheynst, “On the exponential convergence of matching pursuits in quasi-incoherent dictionaries,” IEEE Trans. Inform. Theory, vol. 52, no. 1, pp. 255–261, 2006.
• Barron et al. [2008] A. R. Barron, A. Cohen, W. Dahmen, and R. A. DeVore, “Approximation and learning by greedy algorithms,” The Annals of Statistics, vol. 36, no. 1, p. 64–94, Feb 2008.
• Khanna et al. [2016] R. Khanna, M. Tschannen, and M. Jaggi, “Pursuits in Structured Non-Convex Matrix Factorizations,” arXiv, 2016, 1602.04208v1.
• Rao et al. [2015] N. Rao, P. Shah, and S. Wright, “Forward backward greedy algorithms for atomic norm regularization,” IEEE Transactions on Signal Processing, vol. 63, no. 21, p. 5798–5811, Nov 2015.
• Nemhauser et al. [1978] G. L. Nemhauser, L. A. Wolsey, and M. L. Fisher, “An analysis of approximations for maximizing submodular set functions—i,” Mathematical Programming, vol. 14, no. 1, pp. 265–294, 1978.
• Feige [1998] U. Feige, “A threshold of ln n for approximating set cover,” Journal of the ACM (JACM), vol. 45, no. 4, pp. 634–652, 1998.
• Elenberg et al. [2016] E. Elenberg, R. Khanna, A. Dimakis, and S. Negahban, “Restricted Strong Convexity Implies Weak Submodularity,” Proc. NIPS Workshop on Learning in High Dimensions with Structure, Dec. 2016.
• Khanna et al. [2017] R. Khanna, E. Elenberg, A. Dimakis, S. Neghaban, and J. Ghosh, “Scalable Greedy Support Selection via Weak Submodularity,” AISTATS, 2017.
• Das and Kempe [2011] A. Das and D. Kempe, “Submodular meets Spectral: Greedy Algorithms for Subset Selection, Sparse Approximation and Dictionary Selection,” in ICML, Feb. 2011.
• Negahban et al. [2012] S. Negahban, P. Ravikumar, B. Yu, and M. J. Wainwright, “A Unified Framework for High-Dimensional Analysis of M-Estimators with Decomposable Regularizers,” Statistica Sinica, vol. 27, no. 4, pp. 538–557, 2012.
• Loh and Wainwright [2015] P.-L. Loh and M. J. Wainwright, “Regularized m-estimators with nonconvexity: Statistical and algorithmic theory for local optima,” J. Mach. Learn. Res., vol. 16, no. 1, pp. 559–616, Jan. 2015.
• Shi and Malik [2000] J. Shi and J. Malik, “Normalized cuts and image segmentation,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 22, no. 8, pp. 888–905, Aug. 2000.
• Ng et al. [2001]

A. Y. Ng, M. I. Jordan, and Y. Weiss, “On spectral clustering: Analysis and an algorithm,” in

ADVANCES IN NEURAL INFORMATION PROCESSING SYSTEMS.   MIT Press, 2001, pp. 849–856.
• Collins et al. [2001]

M. Collins, S. Dasgupta, and R. E. Schapire, “A generalization of principal component analysis to the exponential family,” in

Advances in Neural Information Processing Systems.   MIT Press, 2001.
• Mikolov et al. [2013a]

T. Mikolov, I. Sutskever, K. Chen, G. S. Corrado, and J. Dean, “Distributed representations of words and phrases and their compositionality,” in

Advances in Neural Information Processing Systems 26, C. J. C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K. Q. Weinberger, Eds.   Curran Associates, Inc., 2013, pp. 3111–3119.
• Levy and Goldberg [2014] O. Levy and Y. Goldberg, “Neural word embedding as implicit matrix factorization,” in Advances in Neural Information Processing Systems 27, Z. Ghahramani, M. Welling, C. Cortes, N. D. Lawrence, and K. Q. Weinberger, Eds.   Curran Associates, Inc., 2014, pp. 2177–2185.
• Finkelstein et al. [2001] L. Finkelstein, E. Gabrilovich, Y. Matias, E. Rivlin, Z. Solan, G. Wolfman, and E. Ruppin, “Placing search in context: the concept revisited,” 2001, pp. 406–414.
• Bruni et al. [2012] E. Bruni, G. Boleda, M. Baroni, and N.-K. Tran, “Distributional semantics in technicolor,” in Proceedings of the 50th Annual Meeting of the Association for Computational Linguistics: Long Papers - Volume 1, ser. ACL ’12.   Stroudsburg, PA, USA: Association for Computational Linguistics, 2012, pp. 136–145.
• Mikolov et al. [2013b] T. Mikolov, S. W.-t. Yih, and G. Zweig, “Linguistic regularities in continuous space word representations,” in Proceedings of the 2013 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies (NAACL-HLT-2013).   Association for Computational Linguistics, May 2013.
• Mikolov et al. [2013c] T. Mikolov, K. Chen, G. Corrado, and J. Dean, “Efficient estimation of word representations in vector space,” CoRR, vol. abs/1301.3781, 2013.

## Appendix A Supplement

In this section, we provide the missing proofs.

### a.1 Proof of Theorem 2

###### Proof.

An important aspect of the assumptions is that the space of atoms spanned by is orthogonal to the span of . Furthermore, . Let . We will first upper bound the denominator in the submodularity ratio. From strong concavity,

 m¯k2∥B(L∪S)−B(L)∥2F

Rearranging

 0≤ℓ(B(L∪S))−ℓ(B(L)) ≤⟨∇ℓ(B(L)),B(L∪S)−B(L)⟩−m¯k2∥B(L∪S)−B(L)∥2F

where the last equality holds because . Solving the argmax problem, we get . Plugging in, we get,

 ℓ(B(L∪S))−ℓ(B(L))≤12m¯k∥PUS(∇ℓ(B(L)))PVS∥2F

We next bound the numerator. Recall that the atoms in are orthogonal to each other i.e. and are both orthonormal.

For clarity, we define the shorthand, , for

With an arbitrary , and arbitrary scalars for ,

 ℓ(B(L∪{i}))−ℓ(B(L)) ≥ℓ(B(L)+αiiB(L∪S)ii+∑j∈LαijB(L∪S)ij+∑j∈LαjiB(L∪S)ji)−ℓ(B(L)) ≥⟨∇ℓ(B(L)),αiiB(L∪S)ii+∑j∈LαijB(L∪S)ij+∑j∈LαjiB(L∪S)ji⟩ −~M12⎡⎣α2ii∥B(L∪S)ii∥2F+∑j∈Lα2ij∥B(L∪S)ij∥2F+∑j∈Lα2ji∥B(L∪S)ji∥2F⎤⎦.

where the last inequality follows by setting for , and for .

Summing up for all , we get

 ∑i∈Sℓ(B(L∪{i}))−ℓ(B(L)) ≥∑i∈S⎡⎢⎣⟨∇ℓ(B(L)),B(L∪S)ii⟩22~M1∥B(L∪S)ii∥2F+∑j∈L⎛⎜⎝⟨∇ℓ(B(L)),B(L∪S)ij⟩22~M1∥B(L∪S)ij∥2F+⟨∇ℓ(B(L)),B(L∪S)ji⟩22~M1∥B(L∪S)ji∥2F⎞⎟⎠⎤⎥⎦ =12~M1∥PUS∇ℓ(B(L))PVS∥2F

### a.2 Proofs for greedy improvement

Let be the support set formed by Algorithm 1 at iteration . Define with as the greedy improvement. We also define to be the remaining amount to improve, where is the optimum -sized solution. We provide an auxiliary Lemma that uses the submodularity ratio to lower bound the greedy improvement in terms of best possible improvement from step .

###### Lemma 1.

At iteration , the incremental gain of the greedy method (Algorithm 1) is

 A(i+1)≥τγSGi,rrB(i).
###### Proof.

Let . Let be the sequential orthogonalization of the atoms in relative to . Thus,

 rA(i+1)≥|SR|A(i+1) ≥τ|SR|maxa∈SRf(S∪{a})−f(S) ≥τ∑a∈SR[f(S∪{a})−f(S)] ≥τγS,|SR|[f(S∪SR)−f(S)] ≥τγS,|SR|B(i)

Note that the last inequality follows because . The penultimate inequality follows by the definition of weak submodularity, which applies in this case because the atoms in are orthogonal to eachother and are also orthogonal to . ∎

Using Lemma 1, one can prove an approximation guarantee for Algorithm 1.

#### a.2.1 Proof of Theorem 3

###### Proof.

From the notation used for Lemma 1, . Let . From Lemma 1, we have,

 B(i+1)≤(1−C)B(i)≤(1−C)i+1B(0).

From its definition, . So we get,

 [f(S⋆)−f(∅)]−[f(SGi)−f(∅)]≤(1−C)i[f(S⋆)−f(∅)] ⟹ [f(SGi)−f(∅)]≥(1−(1−C)i)[f(S⋆)−f(∅)]≥⎛⎜⎝1−1eτγSGi,rkr⎞⎟⎠[f(S⋆)−f(∅)]

from which the result follows. ∎

### a.3 Proof for GECO bounds

Let be the support set selected by the GECO procedure (Algorithm 2) at iteration . Similar to the section on greedy improvement, we define some notation. Let be the improvement made at step , and as before we have be the remaining amount to improve.

We prove the following auxiliary lemma which lower bounds the gain after adding the atom selected by the subroutine OMPSelin terms of operator norm of the gradient of the current iterate and smoothness of the function.

###### Lemma 2.

Assume that is -strongly concave and -smooth over matrices of in the set . Then,

 D(i+1)≥τmr+kr~M1B(i).
###### Proof.

For simplicity, say . Recall that for a given support set , i.e. we denote by the argmax for for a given support set . Hence, by the optimality of ,

 D(i+1) =ℓ(B(L∪{i}))−ℓ(B(L)) ≥ℓ(B(L)+αuv⊤)−ℓ(B(L))

for an arbitrary , and the vectors selected by OMPSel. Using the smoothness of the , we get,

 D(i+1)≥α⟨∇ℓ(B(L)),uv⊤⟩−α2~M12

Putting in , and by -optimality of OMPSel, we get,

 D(i+1)≥τ22~M1∥∇ℓ(B(L))∥22

Let be obtained from after sequentially orthogonalizing w.r.t. . By definition of the operator norm, we further get,

 D(i+1) ≥τ22~M1∥∇ℓ(B(L))∥22 ≥τ22r~M1∑i∈SR⟨uiv⊤i,∇ℓ(B(L))⟩2 =∥PUSR∇ℓ(B(L))PVSR∥2F ≥τ2mr+kr~M1(ℓ(BL∪SR)−ℓ(B(L))) ≥τ2mr+kr~M1(ℓ(BS⋆)−ℓ(B(L))) =τ2mr+kr~M1B(i)

The proof for Theorem 4 from Lemma 2 now follows using the same steps as for Theorem 3 from Lemma 2.

### a.4 Proof for recovery bounds

#### a.4.1 Proof of Theorem 5

For clarity of representation, let , and for an arbitrary , let , and . Note that has rank atmost . Recall that by the RSC (Definition 3),

 ℓ(B(Sk))−ℓ(Br)−⟨∇ℓ(Br),Δ⟩≤−mk+r2∥Δ∥2F.

From the approximation guarantee, we have,

 ℓ(B(Sk))−ℓ(Br)≥(1−C)[ℓ(0)−ℓ(Br)] ⟹ ℓ(B(Sk))−ℓ(Br)−⟨∇ℓ(Br),Δ⟩≥(1