LOCAL SEARCH ALGORITHMS FOR RANK- CONSTRAINED CONVEX OPTIMIZATION
Year: 2021

Abstract:
We propose greedy and local search algorithms for rank-constrained convex optimization, namely solving min rank(A)≤r∗ R(A) given a convex function R : Rm×n → R and a parameter r∗. These algorithms consist of repeating two steps: (a) adding a new rank-1 matrix to A and (b) enforcing the rank constraint on A. We refine and improve the theoretical analysis of Shalev-Shwartz et al. (2011), and show that if the rank-restricted condition number of R is κ, a solution A with rank O(r∗ · min{κ log R(0)−R(A ∗) , κ }) and R(A) ≤ R(A∗) + can be recovered, where A∗ is the optimal solution. This significantly generalizes associated results on sparse convex optimization, as well as rank-constrained convex optimization for smooth functions. We then introduce new practical variants of these algorithms that have superior runtime and recover better solutions in practice. We demonstrate the versatility of these methods on a wide range of applications involving matrix completion and robust principal component analysis.

text: rank(A)≤r∗ R(A) given a convex function R : Rm×n →
R and a parameter r∗. These algorithms consist of repeating two steps: (a) adding a new rank-1 matrix to A and (b) enforcing the rank constraint on A. We refine and improve the theoretical analysis of Shalev-Shwartz et al. (2011), and show that if the rank-restricted condition number of R is κ, a solution A with rank O(r∗ · min{κ log R(0)−R(A ∗) , κ
2}) and R(A) ≤ R(A∗) + can be recovered, where A∗ is the optimal solution. This significantly generalizes associated results on sparse convex optimization, as well as rank-constrained convex optimization for smooth functions. We then introduce new practical variants of these algorithms that have superior runtime and recover better solutions in practice. We demonstrate the versatility of these methods on a wide range of applications involving matrix completion and robust principal component analysis.

heading: 1 INTRODUCTION
text: Given a real-valued convex function R : Rm×n → R on real matrices and a parameter r∗ ∈ N, the rank-constrained convex optimization problem consists of finding a matrix A ∈ Rm×n that minimizes R(A) among all matrices of rank at most r∗:
min rank(A)≤r∗ R(A) (1)
Even though R is convex, the rank constraint makes this problem non-convex. Furthermore, it is known that this problem is NP-hard and even hard to approximate (Natarajan (1995); Foster et al. (2015)).
In this work, we propose efficient greedy and local search algorithms for this problem. Our contribution is twofold:
1. We provide theoretical analyses that bound the rank and objective value of the solutions returned by the two algorithms in terms of the rank-restricted condition number, which is the natural generalization of the condition number for low-rank subspaces. The results are significantly stronger than previous known bounds for this problem.
2. We experimentally demonstrate that, after careful performance adjustments, the proposed general-purpose greedy and local search algorithms have superior performance to other methods, even for some of those that are tailored to a particular problem. Thus, these algorithms can be considered as a general tool for rank-constrained convex optimization and a viable alternative to methods that use convex relaxations or alternating minimization.
The rank-restricted condition number Similarly to the work in sparse convex optimization, a restricted condition number quantity has been introduced as a reasonable assumption on R. If we let ρ+r be the maximum smoothness bound and ρ − r be the minimum strong convexity bound only along rank-r directions of R (these are called rank-restricted smoothness and strong convexity respectively), the rank-restricted condition number is defined as κr =
ρ+r ρ−r . If this quantity is bounded,
one can efficiently find a solutionA withR(A) ≤ R(A∗)+ and rank r = O(r∗ ·κr+r∗ R(0) ) using a greedy algorithm (Shalev-Shwartz et al. (2011)). However, this is not an ideal bound since the rank scales linearly with R(0) , which can be particularly high in practice. Inspired by the analogous literature on sparse convex optimization by Natarajan (1995); Shalev-Shwartz et al. (2010); Zhang (2011); Jain et al. (2014) and more recently Axiotis & Sviridenko (2020), one would hope to achieve a logarithmic dependence or no dependence at all on R(0) . In this paper we achieve this goal by providing an improved analysis showing that the greedy algorithm of Shalev-Shwartz et al. (2011) in fact returns a matrix of rank of r = O(r∗ · κr+r∗ log R(0) ). We also provide a new local search algorithm together with an analysis guaranteeing a rank of r = O(r∗ · κ2r+r∗). Apart from significantly improving upon previous work on rank-restricted convex optimization, these results directly generalize a lot of work in sparse convex optimization, e.g. Natarajan (1995); Shalev-Shwartz et al. (2010); Jain et al. (2014). Our algorithms and theorem statements can be found in Section 2.
Runtime improvements Even though the rank bound guaranteed by our theoretical analyses is adequate, the algorithm runtimes leave much to be desired. In particular, both the greedy algorithm of Shalev-Shwartz et al. (2011) and our local search algorithm have to solve an optimization problem in each iteration in order to find the best possible linear combination of features added so far. Even for the case that R(A) = 12 ∑ (i,j)∈Ω (M − A)2ij , this requires solving a least squares problem on |Ω| examples and r2 variables. For practical implementations of these algorithms, we circumvent this issue by solving a related optimization problem that is usually much smaller. This instead requires solving n least squares problems with total number of examples |Ω|, each on r variables. This not only reduces the size of the problem by a factor of r, but also allows for a straightforward distributed implementation. Interestingly, our theoretical analyses still hold for these variants. We propose an additional heuristic that reduces the runtime even more drastically, which is to only run a few (less than 10) iterations of the algorithm used for solving the inner optimization problem. Experimental results show that this modification not only does not significantly worsen results, but for machine learning applications also acts as a regularization method that can dramatically improve generalization. These matters, as well as additional improvements for making the local search algorithm more practical, are addressed in Section 2.3.
Roadmap In Section 2, we provide the descriptions and theoretical results for the algorithms used, along with several modifications to boost performance. In Section 3, we evaluate the proposed greedy and local search algorithms on optimization problems like robust PCA. Then, in Section 4 we evaluate their generalization performance in machine learning problems like matrix completion.

heading: 2 ALGORITHMS & THEORETICAL GUARANTEES
text: In Sections 2.1 and 2.2 we state and provide theoretical performance guarantees for the basic greedy and local search algorithms respectively. Then in Section 2.3 we state the algorithmic adjustments that we propose in order to make the algorithms efficient in terms of runtime and generalization performance. A discussion regarding the tightness of the theoretical analysis is deferred to Appendix A.4.
When the dimension is clear from context, we will denote the all-ones vector by 1, and the vector that is 0 everywhere and 1 at position i by 1i. Given a matrix A, we denote by im(A) its column span. One notion that we will find useful is that of singular value thresholding. More specifically,
given a rank-k matrix A ∈ Rm×n with SVD k∑ i=1 σiu ivi> such that σ1 ≥ · · · ≥ σk, as well as an
integer parameter r ≥ 1, we define Hr(A) = r∑ i=1 σiu ivi> to be the operator that truncates to the r highest singular values of A.

heading: 2.1 GREEDY
text: Algorithm 1 (Greedy) was first introduced in Shalev-Shwartz et al. (2011) as the GECO algorithm. It works by iteratively adding a rank-1 matrix to the current solution. This matrix is chosen as the
rank-1 matrix that best approximates the gradient, i.e. the pair of singular vectors corresponding to the maximum singular value of the gradient. In each iteration, an additional procedure is run to optimize the combination of previously chosen singular vectors.
In Shalev-Shwartz et al. (2011) guarantee on the rank of the solution returned by the algorithm is r∗κr+r∗ R(0) . The main bottleneck in order to improve on the R(0) factor is the fact that the analysis is done in terms of the squared nuclear norm of the optimal solution. As the worst-case discrepancy between the squared nuclear norm and the rank is R(0)/ , their bounds inherit this factor. Our analysis works directly with the rank, in the spirit of sparse optimization results (e.g. Shalev-Shwartz et al. (2011); Jain et al. (2014); Axiotis & Sviridenko (2020)). A challenge compared to these works is the need for a suitable notion of “intersection” between two sets of vectors. The main technical contribution of this work is to show that the orthogonal projection of one set of vectors into the span of the other is such a notion, and, based on this, to define a decomposition of the optimal solution that is used in the analysis.
Algorithm 1 Greedy 1: procedure GREEDY(r ∈ N : target rank) 2: function to be minimized R : Rm×n → R 3: U ∈ Rm×0 . Initially rank is zero 4: V ∈ Rn×0 5: for t = 0 . . . r − 1 do 6: σuv> ← H1(∇R(UV >)) . Max singular value σ and corresp. singular vectors u, v 7: U ← (U u) . Append new vectors as columns 8: V ← (V v) 9: U, V ← OPTIMIZE(U, V ) 10: return UV > 11: procedure OPTIMIZE(U ∈ Rm×r, V ∈ Rn×r) 12: X ← arg min
X∈Rr×r R(UXV >)
13: return UX, V
Theorem 2.1 (Algorithm 1 (greedy) analysis). Let A∗ be any fixed optimal solution of (1) for some function R and rank bound r∗, and let > 0 be an error parameter. For any integer r ≥ 2r∗ · κr+r∗ log R(0)−R(A∗) , if we let A = GREEDY(r) be the solution returned by Algorithm 1, then R(A) ≤ R(A∗) + . The number of iterations is r.
The proof of Theorem 2.1 can be found in Appendix A.2.

heading: 2.2 LOCAL SEARCH
text: One drawback of Algorithm 1 is that it increases the rank in each iteration. Algorithm 2 is a modification of Algorithm 1, in which the rank is truncated in each iteration. The advantage of Algorithm 2 compared to Algorithm 1 is that it is able to make progress without increasing the rank of A, while Algorithm 1 necessarily increases the rank in each iteration. More specifically, because of the greedy nature of Algorithm 1, some rank-1 components that have been added to A might become obsolete or have reduced benefit after a number of iterations. Algorithm 2 is able to identify such candidates and remove them, thus allowing it to continue making progress.
Theorem 2.2 (Algorithm 2 (local search) analysis). Let A∗ be any fixed optimal solution of (1) for some function R and rank bound r∗, and let > 0 be an error parameter. For any integer r ≥ r∗ · (1 + 8κ2r+r∗), if we let A = LOCAL SEARCH(r) be the solution returned by Algorithm 2, then R(A) ≤ R(A∗) + . The number of iterations is O ( r∗κr+r∗ log R(0)−R(A∗) ) .
The proof of Theorem 2.2 can be found in Appendix A.3.
Algorithm 2 Local Search 1: procedure LOCAL SEARCH(r ∈ N : target rank) 2: function to be minimized R : Rm×n → R 3: U ← 0m×r . Initialize with all-zero solution 4: V ← 0n×r 5: for t = 0 . . . L− 1 do . Run for L iterations 6: σuv> ← H1(∇R(UV >)) . Max singular value σ and corresp. singular vectors u, v 7: U, V ← TRUNCATE(U, V ) . Reduce rank of UV > by one 8: U ← (U u) . Append new vectors as columns 9: V ← (V v) 10: U, V ← OPTIMIZE(U, V ) 11: return UV > 12: procedure TRUNCATE(U ∈ Rm×r, V ∈ Rn×r) 13: UΣV > ← SVD(Hr−1(UV >)) . Keep all but minimum singular value 14: return UΣ, V

heading: 2.3 ALGORITHMIC ADJUSTMENTS
text: Inner optimization problem The inner optimization problem that is used in both greedy and local search is:
min X∈Rr×r
R(UXV >) . (2)
It essentially finds the choice of matrices U ′ and V ′, with columns in the column span of U and V respectively, that minimizes R(U ′V ′>). We, however, consider the following problem instead:
min V ∈Rn×r
R(UV >) . (3)
Note that the solution recovered from (3) will never have worse objective value than the one recovered from (2), and that nothing in the analysis of the algorithms breaks. Importantly, (3) can usually be solved much more efficiently than (2). As an example, consider the following objective that appears in matrix completion: R(A) = 12 ∑ (i,j)∈Ω (M − A)2ij for some Ω ⊆ [m] × [n]. If we let ΠΩ(·) be an operator that zeroes out all positions in the matrix that are not in Ω, we have ∇R(A) = −ΠΩ(M −A). The optimality condition of (2) now is U>ΠΩ(M −UXV >)V = 0 and that of (3) is U>ΠΩ(M − UV >) = 0. The former corresponds to a least squares linear regression problem with |Ω| examples and r2 variables, while the latter can be decomposed into n independent
systems U> ( ∑ i:(i,j)∈Ω 1i1 > i ) UV j = U>ΠΩ (M1j), where the variable is V j which is the j-th column of V . The j-th of these systems corresponds to a least squares linear regression problem with |{i : (i, j) ∈ Ω}| examples and r variables. Note that the total number of examples in all systems is
∑ j∈[n] |{i : (i, j) ∈ Ω}| = |Ω|. The choice of V here as the variable to be optimized is
arbitrary. In particular, as can be seen in Algorithm 3, in practice we alternate between optimizing U and V in each iteration. It is worthy of mention that the OPTIMIZE FAST procedure is basically the same as one step of the popular alternating minimization procedure for solving low-rank problems. As a matter of fact, when our proposed algorithms are viewed from this lens, they can be seen as alternating minimization interleaved with rank-1 insertion and/or removal steps.
Singular value decomposition As modern methods for computing the top entries of a singular value decomposition scale very well even for large sparse matrices (Martinsson et al. (2011); Szlam et al. (2014); Tulloch (2014)), the “insertion” step of greedy and local search, in which the top entry of the SVD of the gradient is determined, is quite fast in practice. However, these methods are not suited for computing the smallest singular values and corresponding singular vectors, a step required for the local search algorithm that we propose. Therefore, in our practical implementations we opt to perfom the alternative step of directly removing one pair of vectors from the representation UV >. A simple approach is to go over all r possible removals and pick the one that increases the
Algorithm 3 Fast inner Optimization
1: procedure OPTIMIZE FAST(U ∈ Rm×r, V ∈ Rn×r, t ∈ N : iteration index of algorithm) 2: if t mod 2 = 0 then 3: X ← arg min
X∈Rm×r R(XV >)
4: return X,V 5: else 6: X ← arg min
X∈Rn×r R(UX>)
7: return U,X
objective by the least amount. A variation of this approach has been used by Shalev-Shwartz et al. (2011). However, a much faster approach is to just pick the pair of vectors U1i, V 1i that minimizes ‖U1i‖2‖V 1i‖2. This is the approach that we use, as can be seen in Algorithm 4.
Algorithm 4 Fast rank reduction
1: procedure TRUNCATE FAST(U ∈ Rm×r, V ∈ Rn×r) 2: i← arg min
i∈[r] ‖U1i‖2‖V 1i‖2 3: return ( U[m],[1,i−1] U[m],[i+1,r] ) , ( V[n],[1,i−1] V[n],[i+1,r] ) . Remove column i
After the previous discussion, we are ready to state the fast versions of Algorithm 1 and Algorithm 2 that we use for our experiments. These are Algorithm 2.3 and Algorithm 5. Notice that we initialize Algorithm 5 with the solution of Algorithm 2.3 and we run it until the value ofR(·) stops decreasing rather than for a fixed number of iterations.
Algorithm 2.3 (Fast Greedy). The Fast Greedy algorithm is defined identically as Algorithm 1, with the only difference that it uses the OPTIMIZE FAST routine as opposed to the OPTIMIZE routine.
Algorithm 5 Fast Local Search 1: procedure FAST LOCAL SEARCH(r ∈ N : target rank) 2: function to be minimized R : Rm×n → R 3: U, V ← solution returned by FAST GREEDY(r) 4: do 5: Uprev, Vprev ← U, V 6: σuv> ← H1(∇R(UV >)) . Max singular value σ and corresp. singular vectors u, v 7: U, V ← TRUNCATE FAST(U, V ) . Reduce rank of UV > by one 8: U ← (U u) . Append new vectors as columns 9: V ← (V v) 10: U, V ← OPTIMIZE FAST(U, V, t) 11: while R(UV >) < R(UprevV >prev) 12: return UprevV >prev

heading: 3 OPTIMIZATION APPLICATIONS
text: An immediate application of the above algorithms is in the problem of low rank matrix recovery. Given any convex distance measure between matrices d : Rm×n × Rm×n → R≥0, the goal is to find a low-rank matrix A that matches a target matrix M as well as possible in terms of d:
min rank(A)≤r∗ d(M,A) This problem captures a lot of different applications, some of which we go over in the following sections.

heading: 3.1 LOW-RANK APPROXIMATION ON OBSERVED SET
text: A particular case of interest is when d(M,A) is the Frobenius norm of M − A, but only applied to entries belonging to some set Ω. In other words, d(M,A) = 12‖ΠΩ(M −A)‖ 2 F . We have compared our Fast Greedy and Fast Local Search algorithms with the SoftImpute algorithm of Mazumder et al. (2010) as implemented by Rubinsteyn & Feldman (2016), on the same experiments as in Mazumder et al. (2010). We have solved the inner optimization problem required by our algorithms by the LSQR algorithm Paige & Saunders (1982). More specifically, M = UV > + η ∈ R100×100, where η is some noise vector. We let every entry of U, V, η be i.i.d. normal with mean 0 and the entries of Ω are chosen i.i.d. uniformly at random over the set [100] × [100]. The experiments have three parameters: The true rank r∗ (of UV >), the percentage of observed entries p = |Ω|/104, and the signal-to-noise ratio SNR. We measure the normalized MSE, i.e. ‖ΠΩ(M − A)‖2F /‖ΠΩ(M)‖2F . The results can be seen in Figure 1, where it is illustrated that Fast Local Search sometimes returns significantly more accurate and lower-rank solutions than Fast Greedy, and Fast Greedy generally returns significantly more accurate and lower-rank solutions than SoftImpute.

heading: 3.2 ROBUST PRINCIPAL COMPONENT ANALYSIS (RPCA)
text: The robust PCA paradigm asks one to decompose a given matrix M as L + S, where L is a lowrank matrix and S is a sparse matrix. This is useful for applications with outliers where directly computing the principal components of M is significantly affected by them. For a comprehensive survey on Robust PCA survey one can look at Bouwmans et al. (2018). The following optimization problem encodes the above-stated requirements:
min rank(L)≤r∗ ‖M − L‖0 (4)
where ‖X‖0 is the sparsity (i.e. number of non-zeros) of X . As neither the rank constraint or the `0 function are convex, Candès et al. (2011) replaced them by their usual convex relaxations, i.e. the nuclear norm ‖ · ‖∗ and `1 norm respectively. However, we opt to only relax the `0 function but not the rank constraint, leaving us with the problem:
min rank(L)≤r∗ ‖M − L‖1 (5)
In order to make the objective differentiable and thus more well-behaved, we further replace the `1 norm by the Huber loss Hδ(x) = { x2/2 if |x| ≤ δ δ|x| − δ2/2 otherwise , thus getting: minrank(L)≤r∗ ∑ ij Hδ(M − L)ij . This is a problem on which we can directly apply our algorithms. We solve the inner optimization problem by applying 10 L-BFGS iterations.
In Figure 2 one can see an example of foreground-background separation from video using robust PCA. The video is from the BMC 2012 dataset Vacavant et al. (2012). In this problem, the lowrank part corresponds to the background and the sparse part to the foreground. We compare three
algorithms: Our Fast Greedy algorithm, standard PCA with 1 component (the choice of 1 was picked to get the best outcome), and the standard Principal Component Pursuit (PCP) algorithm (Candès et al. (2011)), as implemented in Lin et al. (2010), where we tuned the regularization parameter λ to achieve the best result. We find that Fast Greedy has the best performance out of the three algorithms in this sample task.

heading: 4 MACHINE LEARNING APPLICATIONS

heading: 4.1 REGULARIZATION TECHNIQUES
text: In the previous section we showed that our proposed algorithms bring down different optimization objectives aggressively. However, in applications where the goal is to obtain a low generalization error, regularization is needed. We considered two different kinds of regularization. The first method is to run the inner optimization algorithm for less iterations, usually 2-3. Usually this is straightforward since an iterative method is used. For example, in the caseR(A) = 12‖ΠΩ(M−A)‖ 2 F the inner optimization is a least squares linear regression problem that we solve using the LSQR algorithm. The second one is to add an `2 regularizer to the objective function. However, this option did not provide a substantial performance boost in our experiments, and so we have not implemented it.

heading: 4.2 MATRIX COMPLETION WITH RANDOM NOISE
text: In this section we evaluate our algorithms on the task of recovering a low rank matrix UV > after observing ΠΩ(UV > + η), i.e. a fraction of its entries with added noise. As in Section 3.1, we use the setting of Mazumder et al. (2010) and compare with the SoftImpute method. The evaluation metric is the normalized MSE, defined as (
∑ (i,j)/∈Ω (UV > −A)2ij)/( ∑ (i,j)/∈Ω (UV >)2ij), where A is the
predicted matrix andUV > the true low rank matrix. A few example plots can be seen in Figure 3 and a table of results in Table 1. We have implemented the Fast Greedy and Fast Local Search algorithms with 3 inner optimization iterations. In the first few iterations there is a spike in the relative MSE of the algorithms that use the OPTIMIZE FAST routine. We attribute this to the aggressive alternating minimization steps of this procedure and conjecture that adding a regularization term to the objective might smoothen the spike. However, the Fast Local Search algorithm still gives the best overall performance in terms of how well it approximates the true low rank matrix UV >, and in particular with a very small rank—practically the same as the true underlying rank.

heading: 4.3 RECOMMENDER SYSTEMS
text: In this section we compare our algorithms on the task of movie recommendation on the Movielens datasets Harper & Konstan (2015). In order to evaluate the algorithms, we perform random 80%- 20% train-test splits that are the same for all algorithms and measure the mean RMSE in the test set. If we let Ω ⊆ [m] × [n] be the set of user-movie pairs in the training set, we assume that the true user-movie matrix is low rank, and thus pose (1) with R(A) = 12‖ΠΩ(M − A)‖ 2 F . We make the following slight modification in order to take into account the range of the ratings [1, 5]: We clip the entries of A between 1 and 5 when computing ∇R(A) in Algorithm 2.3 and Algorithm 5. In other words, instead of ΠΩ(A −M) we compute the gradient as ΠΩ(clip(A, 1, 5) −M). This is similar to replacing our objective by a Huber loss, with the difference that we only do so in the steps that we mentioned and not the inner optimization step, mainly for runtime efficiency reasons.
The results can be seen in Table 2. We do not compare with Fast Local Search, as we found that it only provides an advantage for small ranks (< 30), and otherwise matches Fast Greedy. For the inner optimization steps we have used the LSQR algorithm with 2 iterations in the 100K and 1M datasets, and with 3 iterations in the 10M dataset. Note that even though the SVD algorithm by Koren et al. (2009) as implemented by Hug (2020) (with no user/movie bias terms) is a highly tuned algorithm for recommender systems that was one of the top solutions in the famous Netflix prize, it has comparable performance to our general-purpose Algorithm 2.3.
Finally, Table 3 demonstrates the speedup achieved by our algorithms over the basic greedy implementation. It should be noted that the speedup compared to the basic greedy of Shalev-Shwartz et al. (2011) (Algorithm 1) is larger as rank increases, since the fast algorithms scale linearly with rank, but the basic greedy scales quadratically.
It is important to note that our goal here is not to be competitive with the best known algorithms for matrix completion, but rather to propose a general yet practically applicable method for rankconstrained convex optimization. For a recent survey on the best performing algorithms in the Movielens datasets see Rendle et al. (2019). It should be noted that a lot of these algorithms have significant performance boost compared to our methods because they use additional features (meta information about each user, movie, timestamp of a rating, etc.) or stronger models (user/movie biases, ”implicit” ratings). A runtime comparison with these recent approches is an interesting avenue for future work. As a rule of thumb, however, Fast Greedy has roughly the same runtime as SVD (Koren et al. (2009)) in each iteration, i.e. O(|Ω|r), where Ω is the set of observable elements and r is the rank. As some better performing approaches have been reported to be much slower than SVD (e.g. SVD++ is reported to be 50-100x slower than SVD in the Movielens 100K and 1M datasets (Hug (2020)), this might also suggest a runtime advantage of our approach compared to some better performing methods.

heading: 5 CONCLUSIONS
text: We presented simple algorithms with strong theoretical guarantees for optimizing a convex function under a rank constraint. Although the basic versions of these algorithms have appeared before, through a series of careful runtime, optimization, and generalization performance improvements that we introduced, we have managed to reshape the performance of these algorithms in all fronts. Via our experimental validation on a host of practical problems such as low-rank matrix recovery with missing data, robust principal component analysis, and recommender systems, we have shown that the performance in terms of the solution quality matches or exceeds other widely used and even specialized solutions, thus making the argument that our Fast Greedy and Fast Local Search routines can be regarded as strong and practical general purpose tools for rank-constrained convex optimization. Interesting directions for further research include the exploration of different kinds of regularization and tuning for machine learning applications, as well as a competitive implementation and extensive runtime comparison of our algorithms.

heading: A APPENDIX
text: A.1 PRELIMINARIES AND NOTATION
Given an positive integer k, we denote [k] = {1, 2, . . . , k}. Given a matrix A, we denote by ‖A‖F its Frobenius norm, i.e. the `2 norm of the entries of A (or equivalently of the singular values of A). The following lemma is a simple corollary of the definition of the Frobenius norm: Lemma A.1. Given two matrices A ∈ Rm×n,B ∈ Rm×n, we have ‖A + B‖2F ≤ 2 ( ‖A‖2F + ‖B‖2F ) .
Proof. ‖A+B‖2F = ∑ ij (A+B)2ij ≤ 2 ∑ ij (A2ij +B 2 ij) = 2(‖A‖2F + ‖B‖2F )
Definition A.2 (Rank-restricted smoothness, strong convexity, condition number). Given a convex function R ∈ Rm×n → R and an integer parameter r, the rank-restricted smoothness of R at rank r is the minimum constant ρ+r ≥ 0 such that for any two matrices A ∈ Rm×n, B ∈ Rm×n such that rank(A−B) ≤ r, we have
R(B) ≤ R(A) + 〈∇R(A), B −A〉+ ρ + r
2 ‖B −A‖2F .
Similarly, the rank-restricted strong convexity of R at rank r is the maximum constant ρ−r ≥ 0 such that for any two matrices A ∈ Rm×n, B ∈ Rm×n such that rank(A−B) ≤ r, we have
R(B) ≥ R(A) + 〈∇R(A), B −A〉+ ρ − r
2 ‖B −A‖2F .
Given that ρ+r , ρ − r exist and are nonzero, the rank-restricted condition number of R at rank r is then defined as
κr = ρ+r ρ−r
Note that ρ+r is increasing and ρ − r is decreasing in r. Therefore, even though our bounds are proven in terms of the constants ρ + 1
ρ−r and ρ
+ 2
ρ−r , these quantities are always at most ρ
+ r
ρ−r = κr as long as r ≥ 2,
and so they directly imply the same bounds in terms of the constant κr.
Definition A.3 (Spectral norm). Given a matrix A ∈ Rm×n, we denote its spectral norm by ‖A‖2. The spectral norm is defined as
‖A‖2 = max x∈Rn ‖Ax‖2 ‖x‖2 ,
Definition A.4 (Singular value thresholding operator). Given a matrix A ∈ Rm×n of rank k, a singular value decomposition A = UΣV > such that Σ11 ≥ Σ22 ≥ · · · ≥ Σkk, and an integer 1 ≤ r ≤ k, we define Hr(A) = UΣ′V >, here Σ′ is a diagonal matrix with
Σ′ii = { Σii if i ≤ r 0 otherwise

heading: In other words, Hr(·) is an operator that eliminates all but the top r highest singular values of a matrix.
text: Lemma A.5 (Weyl’s inequality). For any matrix A and integer i ≥ 1, let σi(A) be the i-th largest singular value ofA or 0 if i > rank(A). Then, for any two matricesA,B and integers i ≥ 1, j ≥ 1:
σi+j−1(A+B) ≤ σi(A) + σj(B)
A proof of the previous fact can be found e.g. in Fisk (1996). Lemma A.6 (Hr(·) optimization problem). Let A ∈ Rm×n be a rank-k matrix and r ∈ [k] be an integer parameter. Then M = 1λHr(A) is an optimal solution to the following optimization problem:
max rank(M)≤r {〈A,M〉 − λ 2 ‖M‖2F } (6)
Proof. Let UΣV > = ∑ i ΣiiUiV > i be a singular value decomposition of A. We note that (6) is equivalent to
min rank(M)≤r
‖A− λM‖2F := f(M) (7)
Now, note that f( 1λHr(A)) = ‖A−Hr(A)‖ 2 F = k∑ i=r+1 Σ2ii. On the other hand, by applying Weyl’s inequality (Lemma A.5) for j = r + 1,
f(M) = ‖A− λM‖2F = k+r∑ i=1 σ2i (A− λM) ≥ k+r∑ i=1 (σi+r(A)− σr+1(λM))2 = k∑ i=r+1 Σ2ii ,
where the last equality follows from the fact that rank(A) = k and rank(M) ≤ r. Therefore, M = 1λHr(A) minimizes (7) and thus maximizes (6).
A.2 PROOF OF THEOREM 2.1 (GREEDY)
We will start with the following simple lemma about the Frobenius norm of a sum of matrices with orthogonal columns or rows: Lemma A.7. Let U ∈ Rm×r, V ∈ Rn×r, X ∈ Rm×r, Y ∈ Rn×r be such that the columns of U are orthogonal to the columns of X or the columns of V are orthogonal to the columns of Y . Then ‖UV > +XY >‖2F = ‖UV >‖2F + ‖XY >‖2F .
Proof. If the columns of U are orthogonal to those of X , then U>X = 0 and if the columns of V are orthogonal to those of Y , then Y >V = 0. Therefore in any case 〈UV >, XY >〉 = Tr(V U>XY >) = Tr(U>XY >V ) = 0, implying
‖UV > +XY >‖2F = ‖UV >‖2F + ‖XY >‖2F + 2〈UV >, XY >〉 = ‖UV >‖2F + ‖XY >‖2F
Additionally, we have the following lemma regarding the optimality conditions of (2):
Lemma A.8. Let A = UXV > where U ∈ Rm×r, X ∈ Rr×r, and V ∈ Rn×r, such that X is the optimal solution to (2). Then for any u ∈ im(U) and v ∈ im(V ) we have that 〈∇R(A), uv>〉 = 0.
Proof. By the optimality condition of 2, we have that
U>∇R(A)V = 0
Now, for any u = Ux and v = V y we have
〈∇R(A), uv>〉 = u>∇R(A)v = x>U>∇R(A)V y = 0
We are now ready for the proof of Theorem 2.1.
Proof. Let At−1 be the current solution UV > before iteration t− 1 ≥ 0. Let u ∈ Rm and v ∈ Rm be left and right singular vectors of matrix ∇R(A), i.e. unit vectors maximizing |〈∇R(A), uv>〉|. Let Bt = {B|B = At−1 + ηuvT , η ∈ R}. By smoothness we have
R(At−1)−R(At) ≥ max B∈Bt {R(At−1)−R(B)}
≥ max B∈Bt
{ −〈∇R(At−1), B −At−1〉 −
ρ+1 2 ‖B −At−1‖2F } ≥ max
η
{ η〈∇R(At−1), uv>〉 − η2
ρ+1 2 } = max
η
{ η‖∇R(At−1)‖2 − η2
ρ+1 2 } = ‖∇R(At−1)‖22
2ρ+1
where ‖ · ‖2 is the spectral norm (i.e. maximum magnitude of a singular value). On the other hand, by strong convexity and noting that
rank(A∗ −At−1) ≤ rank(A∗) + rank(At−1) ≤ r∗ + r ,
R(A∗)−R(At−1) ≥ 〈∇R(At−1), A∗ −At−1〉+ ρ−r+r∗
2 ‖A∗ −At−1‖2F . (8)
Let At−1 = UV > and A∗ = U∗V ∗>. We let Πim(U) = U(U>U)+U> and Πim(V ) = V (V >V )+V > denote the orthogonal projections onto the images of U and V respectively. We now write
A∗ = U∗V ∗> = (U1 + U2)(V 1 + V 2)> = U1V 1> + U1V 2> + U2V ∗>
where U1 = Πim(U)U∗ is a matrix where every column of U∗ is replaced by its projection on im(U) and U2 = U∗−U1 and similarly V 1 = Πim(V )V ∗ is a matrix where every column of V ∗ is replaced by its projection on im(V ) and V 2 = V ∗− V 1. By setting U ′ = (−U | U1) and V ′ = (V | V 1) we can write
A∗ −At−1 = U ′V ′> + U1V 2> + U2V ∗>
where im(U ′) = im(U) and im(V ′) = im(V ). Also, note that
rank(U1V 2>) ≤ rank(V 2) ≤ rank(V ∗) = rank(A∗) ≤ r∗
and similarly rank(U2V ∗>) ≤ r∗. So now the right hand side of (8) can be reshaped as
〈∇R(At−1), A∗ −At−1〉+ ρ−r+r∗
2 ‖A∗ −At−1‖2F
= 〈∇R(At−1), U ′V ′> + U1V 2> + U2V ∗>〉+ ρ−r+r∗
2 ‖U ′V ′> + U1V 2> + U2V ∗>‖2F
Now, note that since by definition the columns of U ′ are in im(U) and the columns of V ′ are in im(V ), Lemma A.8 implies that 〈∇R(At−1), U ′V ′>〉 = 0. Therefore the above is equal to
〈∇R(At−1), U1V 2> + U2V ∗>〉+ ρ−r+r∗
2 ‖U ′V ′> + U1V 2> + U2V ∗>‖2F
≥ 〈∇R(At−1), U1V 2>〉+ 〈∇R(At−1), U2V ∗>〉+ ρ−r+r∗
2
( ‖U1V 2>‖2F + ‖U2V ∗>‖2F ) ≥ 2 min
rank(M)≤r∗
{ 〈∇R(At−1),M〉+ ρ−r+r∗
2 ‖M‖2F
}
= −2‖Hr ∗(∇R(At−1))‖2F
2ρ−r+r∗
≥ −r∗ ‖∇R(At−1)‖ 2 2
ρ−r+r∗
where the first equality follows by noticing that the columns of V ′ and V 1 are orthogonal to those of V 2 and the columns of U ′ and U1 are orthogonal to those of U2, and applying Lemma A.7. The last equality is a direct application of Lemma A.6 and the last inequality states that the largest squared singular value is not smaller than the average of the top r∗ squared singular values. Therefore we have concluded that
‖∇R(At−1)‖22 ≥ ρ−r+r∗
r∗ (R(At−1)−R(A∗))
Plugging this back into the smoothness inequality, we get
R(At−1)−R(At) ≥ 1
2r∗κ (R(At−1)−R(A∗))
or equivalently
R(At)−R(A∗) ≤ (
1− 1 2r∗κ
) (R(At−1)−R(A∗)) .
Therefore after L = 2r∗κ log R(A0)−R(A ∗) iterations we have
R(AT )−R(A∗) ≤ (
1− 1 2r∗κ
)L (R(A0)−R(A∗))
≤ e− L2r∗κ (R(A0)−R(A∗)) ≤
Since A0 = 0, the result follows.
A.3 PROOF OF THEOREM 2.2 (LOCAL SEARCH)
Proof. Similarly to Section A.3, we let At−1 be the current solution before iteration t − 1 ≥ 0. Let u ∈ Rm and v ∈ Rm be left and right singular vectors of matrix ∇R(A), i.e. unit vectors maximizing |〈∇R(A), uv>〉| and let
Bt = {B|B = At−1 + ηuvT − σminxy>, η ∈ R},
where σminxy> = At−1 −Hr−1(At−1) is the rank-1 term corresponding to the minimum singular value of At−1. By smoothness we have
R(At−1)−R(At) ≥ max B∈Bt {R(At−1)−R(B)}
≥ max B∈Bt
{ −〈∇R(At−1), B −At−1〉 −
ρ+2 2 ‖B −At−1‖2F } = max
η∈R
{ −〈∇R(At−1), ηuv> − σminxy>〉 −
ρ+2 2 ‖ηuv> − σminxy>‖2F } ≥ max
η∈R
{ −〈∇R(At−1), ηuv>〉 − η2ρ+2 − σ2minρ + 2 } = max
η∈R
{ η‖∇R(At−1)‖2 − η2ρ+2 − σ2minρ + 2 } = ‖∇R(At−1)‖22
4ρ+2 − σ2minρ+2 ,
where in the last inequality we used the fact that 〈∇R(At−1), xy>〉 = 0 following from Lemma A.8, as well as Lemma A.1.
On the other hand, by strong convexity,
R(A∗)−R(At−1) ≥ 〈∇R(At−1), A∗ −At−1〉+ ρ−r+r∗
2 ‖A∗ −At−1‖2F .
Let At−1 = UV > and A∗ = U∗V ∗>. We write
A∗ = U∗V ∗> = (U1 + U2)(V 1 + V 2)> = U1V 1> + U1V 2> + U2V ∗>
where U1 is a matrix where every column of U∗ is replaced by its projection on im(U) and U2 = U∗ − U1 and similarly V 1 is a matrix where every column of V ∗ is replaced by its projection on im(V ) and V 2 = V ∗ − V 1. By setting U ′ = (−U | U1) and V ′ = (V | V 1) we can write
A∗ −At−1 = U ′V ′> + U1V 2> + U2V ∗>
where im(U ′) = im(U) and im(V ′) = im(V ). Also, note that
rank(U1V 2>) ≤ rank(V 2) ≤ rank(V ∗) = rank(A∗) ≤ r∗
and similarly rank(U2V ∗>) ≤ r∗. So we now have
〈∇R(At−1), A∗ −At−1〉+ ρ−r+r∗
2 ‖A∗ −At−1‖2F
= 〈∇R(At−1), U ′V ′> + U1V 2> + U2V ∗>〉+ ρ−r+r∗
2 ‖U ′V ′> + U1V 2> + U2V ∗>‖2F
= 〈∇R(At−1), U1V 2> + U2V ∗>〉+ ρ−r+r∗
2 ‖U ′V ′> + U1V 2> + U2V ∗>‖2F
= 〈∇R(At−1), U1V 2> + U2V ∗>〉+ ρ−r+r∗
2
( ‖U ′V ′>‖2F + ‖U1V 2>‖2F + ‖U2V ∗>‖2F )
≥ 〈∇R(At−1), U1V 2>〉+ 〈∇R(At−1), U2V ∗>〉+ ρ−r+r∗
2
( ‖U1V 2>‖2F + ‖U2V ∗>‖2F ) + ρ−r+r∗
2 ‖U ′V ′>‖2F
≥ 2 min rank(M)≤r∗
{ 〈∇R(At−1),M〉+ ρ−r+r∗
2 ‖M‖2F
} + ρ−r+r∗
2 ‖U ′V ′>‖2F
= −2‖Hr ∗(∇R(At−1))‖2F 2ρ−r+r∗ + ρ−r+r∗ 2 ‖U ′V ′>‖2F
≥ −r∗ ‖∇R(At−1)‖ 2 2 ρ−r+r∗ + ρ−r+r∗ 2 ‖U ′V ′>‖2F
where the second equality follows from the fact that 〈∇R(At−1), uv>〉 = 0 for any u ∈ im(U), v ∈ im(V ), the third equality from the fact that im(U2) ⊥ im(U ′) ∪ im(U1) and im(V 2) ⊥ im(V ′) and by applying Lemma A.7, and the last inequality from the fact that the largest squared singular value is not smaller than the average of the top r∗ squared singular values. Now, note that since rank(U1V 1>) ≤ r∗ < r = rank(UV >),
‖U ′V ′>‖2F = ‖U1V 1> − UV >‖2F
= r∑ i=1 σ2i (U 1V 1> − UV >)
≥ r∑ i=1 (σi+r∗(UV >)− σr∗+1(U1V 1>))2
= r∑ i=r∗+1 σ2i (UV >)
≥ (r − r∗)σ2min(UV >) = (r − r∗)σ2min(At−1) ,
where we used the fact that rank(U1V 1>) ≤ r∗ together with Lemma A.5. Therefore we have concluded that
‖∇R(At−1)‖22 ≥ ρ−r+r∗
r∗ (R(At−1)−R(A∗)) +
(ρ−r+r∗) 2(r − r∗)
2r∗ σ2min
Plugging this back into the smoothness inequality and setting κ̃ = ρ + 2
ρ− r+r∗
, we get
R(At−1)−R(At) ≥ 1
4r∗κ̃ (R(At−1)−R(A∗)) +
( ρ−r+r∗(r − r∗)
8r∗κ̃ − ρ+2
) σ2min(At−1)
≥ 1 4r∗κ̃ (R(At−1)−R(A∗))
as long as r ≥ r∗(1 + 8κ̃2), or equivalently, R(At)−R(A∗) ≤ (
1− 1 4r∗κ̃
) (R(At−1)−R(A∗)) .
Therefore after L = 4r∗κ̃ log R(A0)−R(A ∗) iterations we have
R(AT )−R(A∗) ≤ (
1− 1 4r∗κ̃
)L (R(A0)−R(A∗))
≤ e− L4r∗κ̃ (R(A0)−R(A∗)) ≤
Since A0 = 0 and κ̃ ≤ κr+r∗ , the result follows.
A.4 TIGHTNESS OF THE ANALYSIS
It is important to note that the κr+r∗ factor that appears in the rank bounds of both Theorems 2.1 and 2.2 is inherent in these algorithms and not an artifact of our analysis. In particular, such lower bounds based on the restricted condition number have been previously shown for the problem of sparse linear regression. More specifically, Foster et al. (2015) showed that there is a family of instances in which the analogues of Greedy and Local Search for sparse optimization require the sparsity to be Ω(s∗κ′) for constant error > 0, where s∗ is the optimal sparsity and κ′ is the sparsity-restricted condition number. These instances can be easily adjusted to give a rank lower bound of Ω(r∗κr+r∗) for constant error > 0, implying that the κ dependence in Theorem 2.1 is tight for Greedy. Furthermore, specifically for Local Search, Axiotis & Sviridenko (2020) additionally
showed that there is a family of instances in which the analogue of Local Search for sparse optimization requires a sparsity of Ω(s∗(κ′)2). Adapting these instances to the setting of rank-constrained convex optimization is less trivial, but we conjecture that it is possible, which would lead to a rank lower bound of Ω(r∗κ2r+r∗) for Local Search.
We present the following lemma, which essentially states that sparse optimization lower bounds for Orthogonal Matching Pursuit (OMP, Pati et al. (1993)) (resp. Orthogonal Matching Pursuit with Replacement (OMPR, Jain et al. (2011))) in which the optimal sparse solution is also a global optimum, immediately carry over (up to constants) to rank-constrained convex optimization lower bounds for Greedy (resp. Local Search). Lemma A.9. Let f ∈ Rn → R and x∗ ∈ Rn be an s∗-sparse vector that is also a global minimizer of f . Also, let f have restricted smoothness parameter β at sparsity level s + s∗ for some s ≥ s∗ and restricted strong convexity parameter α at sparsity level s + s∗. Then we can define the rankconstrained problem, with R : Rn×n → R,
min rank(A)≤s∗
R(A) := f(diag(A)) + β
2 ‖A− diag(A)‖2F , (9)
where diag(A) is a vector containing the diagonal of A. R has rank-restricted smoothness at rank s+s∗ at most 2β and rank-restricted strong convexity at rank s+s∗ at least α. Suppose that we run t iterations of OMP (resp. OMPR) starting from a solution x, to get solution x′, and similarly run t iterations of Greedy (resp. Local Search) starting from solution A = diag(x) (where diag(x) is a diagonal matrix with x on the diagonal) to get solution A′. Then A′ is diagonal and diag(A′) = x′. In other words, in this scenario OMP and Greedy (resp. OMPR and Local Search) are equivalent.
Proof. Note that for any solution Â of R we have R(Â) ≥ f(diag(Â)) ≥ f(x∗), with equality only if Â is diagonal. Furthermore, rank(diag(x∗)) ≤ s∗, meaning that diag(x∗) is an optimal solution of (9). Now, given any diagonal solution A of (9) such that A = diag(x), we claim that one step of either Greedy or Local Search keeps it diagonal. This is because
∇R(A) = diag(∇f(x)) + β 2 (A− diag(A)) = diag(∇f(x)) .
Therefore the largest eigenvalue of ∇R(A) has corresponding eigenvector 1i for some i, which implies that the rank-1 component which will be added is a multiple of 1i1>i . For the same reason the rank-1 component removed by Local Search will be a multiple of 1j1>j for some j. Therefore running Greedy (resp. Local Search) on such an instance is identical to running OMP (resp. OMPR) on the diagonal.
Together with the lower bound instances of Foster et al. (2015) (in which the global minimum property is true), it immediately implies a rank lower bound of Ω(r∗κr+r∗) for getting a solution with constant error for rank-constrained convex optimization. On the other hand, the lower bound instances of Axiotis & Sviridenko (2020) give a quadratic lower bound in κ for OMPR. The above lemma cannot be directly applied since the sparse solutions are not global minima, but we conjecture that a similar proof will give a rank lower bound of Ω(r∗κ2r+r∗) for rank-constrained convex optimization with Local Search.
A.5 ADDENDUM TO SECTION 4

Review summaries:
This paper considers solving rank-constrained convex optimization. This is a fairly general problem that contains several special cases such as matrix completion and robust PCA. This paper presents a local search approach along with an interesting theoretical analysis of their approach. Furthermore, this paper provided extensive simulations to validate their approach. Overall, the paper provided solid justification for their approach.