summaries: This paper considers the problem of minimizing a non-smooth matrix function  $f(X)$ over the spectrahedron $S_n$, i.e., positive semidefinite matrices with trace equal to 1. Such optimization problems often arise as the convex relaxation of various low-rank matrix recovery problems. In general solving these problems with projected subgradient descent requires a projection onto the spectraherdon at every iteration, which is expensive because of the need to compute a full SVD at every step.  Under the assumption that the non-smooth objective function $f(X)$ can be written as the pointwise maximum of a family of smooth functions,  the authors reformulate the original optimization problem as a saddle point problem which can be solved via the projected extragradient method. The key contribution is that under a “generalized strict complementarity” condition and when close to the ground truth, the projection onto $S_n$ can be replaced by a low-rank projection. Therefore, the overall algorithm can be significantly more efficient, because only a low-rank SVD is needed. This paper develops an extragradient method for low-rank and nonsmooth matrix optimization problems. Under certain conditions, the proposed method converges to an optimal solution with rate O(1/t). Numerical experiments are provided to support the theoretical results. In this paper, motivated by a new efficient method in solving smooth low-rank matrix optimization problems, the authors propose to extend the method to nonsmooth problems. However, direct generalization is not possible, as is shown for a specific failure case. Thus, the authors propose to introduce some auxiliary parameters and transform the original nonsmooth problem to a smooth saddle-point problem. The authors then extend to the method of this new problem and derive a new algorithm that is suitable for the saddle-point problem. Theoretically, the authors proved the convergence of their proposed method. And empirical results verify the correctness of their proposed method. This paper studies the problem of minimizing a convex, non-smooth objective function of matrices with low-rank constraints. For the smooth counterpart of the same problem, it has been previously shown that SVDs computed during the iterations of algorithms such as projected sub-gradient descent can be replaced by low-rank SVDs, thus improving the computational complexity. This paper studies the conditions under which such a replacement can be used for non-smooth objective functions.   The main result of the paper consists of describing the class of non-smooth functions for which low-rank SVDs can be used. The paper proves theoretically and empirically that for non-smooth functions that have a saddle point structure (can be written as a sum of a smooth function and a point-wise maximum of affine functions), and when the initialization is within a certain radius w.r.t the ground truth, low0rank SVDs lead to the same solution as full SVDs and a O(1/t) convergence can be achieved.   Further, the paper also shows empirically that using low-rank SVDs on non-smooth objectives leads to the same solution as using full SVDs.  The paper considers a class of nonsmooth and low-rank matrix optimization problems. Major contributions include: (1) Discuss strict complementarity for nonsmooth problems and challenges of low-rank projected subgradient steps. (2) Under certain assumptions, a special type of nonsmooth problems can be reformulated as maximum of smooth ones and thereby a minimax form, and develop a projected extragradient method with low-rank projections. (3) Use the proposed method to solve several classical low-rank approximation problems in applications. 
abstractText: Low-rank and nonsmooth matrix optimization problems capture many fundamental tasks in statistics and machine learning. While significant progress has been made in recent years in developing efficient methods for smooth low-rank optimization problems that avoid maintaining high-rank matrices and computing expensive high-rank SVDs, advances for nonsmooth problems have been slow paced. In this paper we consider standard convex relaxations for such problems. Mainly, we prove that under a natural generalized strict complementarity condition and under the relatively mild assumption that the nonsmooth objective can be written as a maximum of smooth functions, the extragradient method, when initialized with a “warm-start” point, converges to an optimal solution with rate O(1/t) while requiring only two low-rank SVDs per iteration. We give a precise trade-off between the rank of the SVDs required and the radius of the ball in which we need to initialize the method. We support our theoretical results with empirical experiments on several nonsmooth low-rank matrix recovery tasks, demonstrating that using simple initializations, the extragradient method produces exactly the same iterates when full-rank SVDs are replaced with SVDs of rank that matches the rank of the (low-rank) ground-truth matrix to be recovered.
authors.name: Dan Garber
authors.name: Atara Kaplan
references.authors: Emmanuel Abbe
references.title: Community detection and stochastic block models: recent developments
references.venue: The Journal of Machine Learning Research,
references.authors: Praneeth Netrapalli UN Niranjan Sujay Sanghavi Animashree Anandkumar Prateek Jain
references.title: Non-convex robust pca
references.venue: Advances in Neural Information Processing Systems,
references.authors: Amir Beck
references.title: First-order methods in optimization
references.authors: Amir Beck Marc Teboulle
references.title: A fast iterative shrinkage-thresholding algorithm for linear inverse problems
references.venue: SIAM journal on imaging sciences,
references.authors: Dimitri P Bertsekas
references.title: Nonlinear programming
references.venue: Athena Scientific,
references.authors: Afonso S. Bandeira Nicolas Boumal Amit Singer
references.title: Tightness of the maximum likelihood semidefinite relaxation for angular synchronization
references.venue: Mathematical Programming,
references.authors: Sébastien Bubeck
references.title: Convex optimization: Algorithms and complexity
references.venue: Foundations and Trends in Machine Learning,
references.authors: Vasileios Charisopoulos Yudong Chen Damek Davis Mateo Díaz Lijun Ding Dmitriy Drusvyatskiy
references.title: Low-rank matrix recovery with composite optimization: good conditioning and rapid convergence
references.venue: Foundations of Computational Mathematics,
references.authors: PXinyang Yi Dohyung Park Yudong Chen Constantine Caramanis
references.title: Fast algorithms for robust pca via gradient descent
references.venue: Advances in Neural Information Processing Systems,
references.authors: Yudong Chen Huan Xu Constantine Caramanis Sujay Sanghavi
references.title: Robust matrix completion and corrupted columns
references.venue: In Proceedings of the 28th International Conference on Machine Learning
references.authors: Frank H. Clarke
references.title: Optimization and Nonsmooth Analysis
references.authors: Alexandre d’Aspremont Laurent El Ghaoui Michael I Jordan Gert RG Lanckriet
references.title: A direct formulation for sparse pca using semidefinite programming
references.venue: SIAM review,
references.authors: Lijun Ding Jicong Fan Madeleine Udell
references.title: kfw: A frank-wolfe style algorithm with stronger subproblem oracles, 2020
references.authors: Dmitriy Drusvyatskiy Adrian S Lewis
references.title: Error bounds, quadratic growth, and linear convergence of proximal methods
references.venue: Mathematics of Operations Research,
references.authors: Francisco Facchinei Jong-Shi Pang
references.title: Finite-Dimensional Variational Inequalities and Complementarity Problems, volume II
references.authors: Dan Garber
references.title: Linear convergence of frank-wolfe for rank-one matrix recovery without strong convexity
references.venue: arXiv preprint arXiv:1912.01467,
references.authors: Dan Garber
references.title: On the convergence of projected-gradient methods with low-rank projections for smooth convex minimization over trace-norm balls and related problems
references.venue: SIAM Journal on Optimization,
references.authors: Dan Garber
references.title: On the convergence of stochastic gradient descent with low-rank projections for convex low-rank matrix problems
references.venue: Conference on Learning Theory, COLT,
references.authors: Dan Garber Atara Kaplan
references.title: Fast stochastic algorithms for low-rank and nonsmooth matrix problems
references.venue: The 22nd International Conference on Artificial Intelligence and Statistics,
references.authors: Dan Garber Atara Kaplan
references.title: On the efficient implementation of the matrix exponentiated gradient algorithm for low-rank matrix optimization
references.venue: arXiv preprint arXiv:2012.10469,
references.authors: Hui Zou Trevor Hastie Robert Tibshirani
references.title: Sparse principal component analysis
references.venue: Journal of Computational and Graphical Statistics,
references.authors: Olga Klopp Karim Lounici Alexandre B Tsybakov
references.title: Robust matrix completion
references.venue: Probability Theory and Related Fields,
references.authors: G.M. Korpelevich
references.title: The extragradient method for finding saddle points and other problems
references.venue: Matecon, 12:747–756,
references.authors: Srinadh Bhojanapalli Anastasios Kyrillidis Sujay Sanghavi
references.title: Dropping convexity for faster semi-definite optimization
references.venue: In 29th Annual Conference on Learning Theory,
references.authors: Xiao Li Zhihui Zhu Anthony Man-Cho So Rene Vidal
references.title: Nonconvex robust low-rank matrix recovery
references.venue: SIAM Journal on Optimization,
references.authors: Francesco Locatello Alp Yurtsevert Olivier Fercoq Volkan Cevhert
references.title: Stochastic frankwolfe for composite convex minimization
references.venue: Advances in Neural Information Processing Systems,
references.authors: Canyi Lu Jinhui Tang Shuicheng Yan Zhouchen Lin
references.title: Generalized nonconvex nonsmooth low-rank minimization
references.venue: In Proceedings of the IEEE conference on computer vision and pattern recognition,
references.authors: Emmanuel J. Candès Xiaodong Li Yi Ma John Wright
references.title: Robust principal component analysis
references.venue: Journal of the ACM,
references.authors: Adel Javanmard Andrea Montanari Federico Ricci-Tersenghi
references.title: Phase transitions in semidefinite relaxations
references.venue: Proceedings of the National Academy of Sciences,
references.authors: Arkadi Nemirovski
references.title: Prox-method with rate of convergence o(1/t) for variational inequalities with lipschitz continuous monotone operators and smooth convex-concave saddle point problems
references.venue: SIAM J. on Optimization,
references.authors: Gergely Odor Yen-Huan Li Alp Yurtsever Ya-Ping Hsieh Quoc Tran-Dinh Marwa El Halabi Volkan Cevher
references.title: Frank-wolfe works for non-lipschitz continuous gradient objectives: scalable poisson phase retrieval
references.venue: In 2016 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP),
references.authors: Stephen Boyd Neal Parikh Eric Chu Borja Peleato Jonathan Eckstein
references.title: Distributed optimization and statistical learning via the alternating direction method of multipliers
references.venue: Foundations and Trends in Machine Learning,
references.authors: John Wright Arvind Ganesh Shankar Rao Yigang Peng Yi Ma
references.title: Robust principal component analysis: Exact recovery of corrupted low-rank matrices via convex optimization
references.venue: Advances in Neural Information Processing Systems
references.authors: Julien Mairal Francis Bach Jean Ponce Guillermo Sapiro
references.title: Online learning for matrix factorization and sparse coding
references.venue: Journal of Machine Learning Research,
references.authors: Emile Richard Pierre-Andr‘e Savalle Nicolas Vayatis
references.title: Estimation of simultaneously sparse and low rank matrices
references.venue: Proceedings of the 29th International Conference on Machine Learning,
references.authors: Kiran K Thekumparampil Prateek Jain Praneeth Netrapalli Sewoong Oh
references.title: Projection efficient subgradient method and optimal nonsmooth frank-wolfe method
references.venue: In Advances in Neural Information Processing Systems,
references.authors: Nicolas Boumal Vladislav Voroninski Afonso S. Bandeira
references.title: Deterministic guarantees for burer-monteiro factorizations of smooth semidefinite programs
references.venue: Communications on Pure and Applied Mathematics,
references.authors: Cun Mu Yuqian Zhang John Wright Donald Goldfarb
references.title: Scalable robust matrix recovery: Frank-wolfe meets proximal methods
references.venue: SIAM Journal on Scientific Computing,
references.authors: Lijun Ding Yingjie Fei Qiantong Xu Chengrun Yang
references.title: Spectral frank-wolfe algorithm: Strict complementarity and linear convergence
references.authors: Quanming Yao James T Kwok Taifeng Wang Tie-Yan Liu
references.title: Large-scale low-rank matrix learning with nonconvex regularizers
references.venue: IEEE transactions on pattern analysis and machine intelligence,
references.authors: Yiqiao Zhong Nicolas Boumal
references.title: Near-optimal bounds for phase synchronization
references.venue: SIAM Journal on Optimization,
references.authors: Zirui Zhou Anthony Man-Cho So
references.title: A unified approach to error bounds for structured convex optimization problems
references.venue: Mathematical Programming,
sections.heading: 1 Introduction
sections.text: Low-rank and nonsmooth matrix optimization problems have many important applications in statistics, machine learning, and related fields, such as sparse PCA [21, 34], robust PCA [28, 33, 2, 9, 38], phase synchronization [41, 6, 29], community detection and stochastic block models [1]1, low-rank and sparse covariance matrix recovery [35], robust matrix completion [22, 10], and more. For many of these problems, convex relaxations, in which one replaces the nonconvex low-rank constraint with a trace-norm constraint, have been demonstrated in numerous papers to be highly effective both in theory (under suitable assumptions) and empirically (see references above). These convex relaxations can be formulated as the following general nonsmooth optimization problem:
  min X∈Sn g(X), (1)
  where g : Sn → R is convex but nonsmooth, and Sn = {X ∈ Sn | Tr(X) = 1, X 0} is the spectrahedron in Sn, Sn being the space of n× n real symmetric matrices. Problem (1), despite being convex, is notoriously difficult to solve in large scale. The simplest and most general approach applicable to it is the projected subgradient method [3, 7], which requires on
  1in [41, 6, 29] and [1] the authors consider SDPs with linear objective function and affine constraints of the form A(X) = b. By incorporating the linear constraints into the objective function via a `2 penalty term of the form λ‖A(X)− b‖2, λ > 0, we obtain a nonsmooth objective function.
  35th Conference on Neural Information Processing Systems (NeurIPS 2021).
  each iteration to compute a Euclidean projection onto the spectrahedron Sn, which in worst case amounts to O(n3) runtime per iteration. In many applications g(X) follows a composite model, i.e., g(X) = h(X) + w(X), where h(·) is convex and smooth and w(·) is convex and nonsmooth but admits a simple structure (e.g., nonsmooth regularizer). For such composite objectives, without the spectrahedron constraint, proximal methods such as FISTA [4] or splitting methods such as ADMM [32] are often very effective. However, with the spectrahderon constraint, all such methods require on each iteration to apply a subprocedure (e.g., computing the proximal mapping) which in worst case amounts to at least O(n3) runtime. A third type of off-the-shelf methods include those which are based on the conditional gradient method and adapted to nonsmooth problems, see for instance [31, 19, 36, 26]. The advantage of such methods is that no expensive high-rank SVD computations are needed. Instead, only a single leading eigenvector computation (i.e., a rank-one SVD) per iteration is required. However, whenever the number of iterations is not small, these methods still require to store high-rank matrices in memory, even when the optimal solution is low-rank. Thus, to conclude, standard first-order methods for Problem (1) require in worst case Ω(n3) runtime per iteration and/or to store high-rank matrices in memory.
  In the recent works [17, 18] it was established that for smooth objective functions, the high-rank SVD computations required for Euclidean projections onto the spectrahedron in standard gradient methods, can be replaced with low-rank SVDs in the close proximity of a low-rank optimal solution. This is significant since the runtime to compute a rank-r SVD of a given n×n matrix using efficient iterative methods typically scales with rn2 (and further improves when the matrix is sparse), instead of n3 for a full-rank SVD. These results depend on the existence of eigen-gaps in the gradient of the optimal solution, which we refer to as a generalized strict complementarity condition. These results also hinge on a unique property of the Euclidean projection onto the spectrahedron. The projection onto the spectrahedron of a matrix X ∈ Sn, which admits an eigen-decomposition X = ∑n i=1 λiviv > i , is given by
  ΠSn [X] = n∑ i=1 max{0, λi − λ}viv>i , (2)
  where λ ∈ R is the unique scalar satisfying ∑n i=1 max{0, λi−λ} = 1. This operation thus truncates all eigenvalues that are smaller than λ, while leaving the eigenvectors unchanged, thereby returning a matrix with rank equal to the number of eigenvalues greater than λ. Importantly, when the projection of X onto Sn is of rank r, only the first r components in the eigen-decomposition of X are required to compute it in the first place, and thus, only a rank-r SVD of X is required. In other words and simplifying, [17, 18] show that under strict complementary, at the proximity of an optimal solution of rank r, the exact Euclidean projection equals the rank-r truncated projection given by:
  Π̂rSn [X] := ΠSn [ r∑ i=1 λiviv > i ] . (3)
  Extending the results of [17, 18] to the nonsmooth setting is difficult since the smoothness assumption is critical to the analysis. Moreover, while [17, 18] rely on certain eigen-gaps in the gradients at optimal points, for nonsmooth problems, since the subdifferential set is often not a singleton, it is not likely that a similar eigen-gap property holds for all subgradients of an optimal solution.
  In this paper we show that under the mild assumption that Problem (1) can be formulated as a smooth convex-concave saddle-point problem, i.e., the nonsmooth term can be written as a maximum over (possibly infinite number of) smooth convex functions, we can obtain results in the spirit of [17, 18]. Concretely, we show that if a generalized strict complementarity (GSC) assumption holds for a low-rank optimal solution (see Assumption 1 in the sequel), the extragradient method for smooth convex-concave saddle-point problems [23, 30] (see Algorithm 1 below), when initialized in the proximity of the optimal solution, converges with its original convergence rate of O(1/t), while requiring only two low-rank SVDs per iteration2. It is important to recall that while the extragradient method requires two SVDs per iteration, it has the benefit of a fast O(1/t) convergence rate, while simpler saddle-point methods such as mirror-descent-based only achieve a O(1/ √ t) rate [7].
  Our contributions can be summarized as follows: 2note that the extradgradient method computes two projected-gradient steps on each iteration, and thus two SVDs are needed per iteration.
  • We prove that even under (standard) strict complementarity, the projected subgradient method, when initialized with a “warm-start", may produce iterates with rank higher than that of the optimal solution. This phenomena further motivates our saddle-point approach. See Lemma 5.
  • We suggest a generalized strict complementarity (GSC) condition for saddle-point problems and prove that when g(·) — the objective function in Problem (1), admits a highly popular saddle-point structure (one which captures all applications we mentioned in this paper), GSC w.r.t. an optimal solution to Problem (1) implies GSC (with the same parameters) w.r.t. a corresponding optimal solution of the equivalent saddle-point problem (the other direction always holds). See Section 3.
  • Main result: we prove that for a smooth convex-concave saddle-point problem and an optimal solution which satisfies GSC, the extragradient method, when initialized with a “warm-start", converges with its original rate of O(1/t) while requiring only two low-rank SVDs per iteration. Moreover, we prove GSC facilitates a precise and powerful tradeoff: increasing the rank of SVD computations (beyond the rank of the optimal solution) can significantly increase the radius of the ball in which the method needs to be initialized. See Theorem 1.
  • We present extensive numerical evidence that demonstrate both the plausibility of the GSC assumption in various tasks, and more importantly, demonstrate that indeed the extragradient method with simple initialization converges correctly (i.e., produces exactly the same sequences of iterates) when the rank of the SVDs used to compute the (truncated) projections matches the rank of the (low-rank) ground-truth matrix to be recovered, instead of naively using full-rank SVDs (as suggested by (2)). See Section 5.
sections.heading: 1.1 Additional related work
sections.text: Since, as in the works [17, 18] mentioned before which deal with smooth objectives, strict complementarity plays a key role in our analysis, we refer the interested reader to the recent works [16, 39, 13, 20] which also exploit this property for efficient smooth and convex optimization over the spectrahedron. Strict complementarity has also played an instrumental role in two recent and very influential works which used it to prove linear convergence rates for proximal gradient methods [42, 14].
  Besides convex relaxations such as Problem (1), considerable advances have been made in the past several yeas in developing efficient nonconvex methods with global convergence guarantees for low-rank matrix problems. In [37] the authors consider semidefinite programs and prove that under a smooth manifold assumption on the constraints, such methods converge to the optimal global solution. In [24] the authors prove global convergence of factorized nonconvex gradient descent from a ”warm-start” initialization point for non-linear smooth minimization on the positive semidefinite cone. Very recently, [8] has established, under statistical conditions, fast convergence results from “warm-start” initialization of nonconvex first-order methods, when applied to nonsmooth nonconvex matrix recovery problems which are based on the explicit factorization of the low-rank matrix. A result of similar flavor concerning nonsmooth and nonconvex formulation of robust recovery of lowrank matrices from random linear measurements was presented in [25]. Finally, several recent works have considered nonconvex low-rank regularizers which result in nonconvex nonsmooth optimization problems, but guarantee convergence only to a stationary point [27, 40].
sections.heading: 2 Strict complementarity for nonsmooth optimization and difficulty of applying low-rank projected subgradient steps
sections.text: Our analysis of the nonsmooth Problem (1) naturally depends on certain subgradients of an optimal solution which, in many aspects, behave like the gradients of smooth functions. The existence of such a subgradient is guaranteed from the first-order optimality condition for constrained convex minimization problems:
  Lemma 1 (first-order optimality condition, see [3]). Let g : Sn → R be a convex function. Then X∗ ∈ Sn minimizes g over Sn if and only if there exists a subgradient G∗ ∈ ∂g(X∗) such that 〈X−X∗,G∗〉 ≥ 0 for all X ∈ Sn.
  For some G∗ ∈ ∂g(X∗) which satisfies the first-order optimality condition for an optimal solution X∗, if the multiplicity of the smallest eigenvalue equals r∗ = rank(X∗), then it can be shown that the optimal solution satisfies a strict complementarity assumption. The equivalence between a standard strict complementarity assumption on some low-rank optimal solution of a smooth optimization problem over the spectrahedron and an eigen-gap in the gradient of the optimal solution was established in [39]. We generalize this equivalence to also include nonsmooth problems. The proof follows similar arguments and is given in Appendix A.1. Definition 1 (strict complementarity). An optimal solution X∗ ∈ Sn of rank r∗ for Problem (1) satisfies the strict complementarity assumption with parameter δ > 0, if there exists an optimal solution of the dual problem3 (Z∗, s∗) ∈ Sn ×R such that rank(Z∗) = n− r∗, and λn−r∗(Z∗) ≥ δ. Lemma 2. Let X∗ ∈ Sn be a rank-r∗ optimal solution to Problem (1). X∗ satisfies the (standard) strict complementarity assumption with parameter δ > 0 if and only if there exists a subgradient G∗ ∈ ∂g(X∗) such that 〈X−X∗,G∗〉 ≥ 0 for all X ∈ Sn and λn−r∗(G∗)− λn(G∗) ≥ δ.
  Throughout this paper we assume a weaker and more general assumption than strict complementarity, namely generalized strict complementarity (GSC), which we present now. Assumption 1 (generalized strict complementarity). We say an optimal solution X∗ to Problem (1) satisfies the generalized strict complementarity assumption with parameters r, δ, if there exists a subgradient G∗ ∈ ∂g(X∗) such that 〈X−X∗,G∗〉 ≥ 0 for all X ∈ Sn and λn−r(G∗)−λn(G∗) ≥ δ.
  In [17] the author presents several characteristic properties of the gradient of the optimal solution in optimization problems over the spectrahedron. Using the existence of subgradients which satisfy the condition in Lemma 1, we can extend these properties also to the nonsmooth setting. The following lemma shows that GSC with parameters (r, δ) for some δ > 0 (Assumption 1) is a sufficient condition for the optimal solution to be of rank at most r. The proof follows immediately from the proof of the analogous Lemma 7 in [17], by replacing the gradient of the optimal solution with a subgradient for which the first-order optimality condition holds. Lemma 3. Let X∗ be an optimal solution to Problem (1) and write its eigen-decomposition as X∗ = ∑r∗ i=1 λiviv T i . Then, any subgradient G
  ∗ ∈ ∂g(X∗) which satisfies 〈X−X∗,G∗〉 ≥ 0 for all X ∈ Sn, admits an eigen-decomposition such that the set of vectors {vi}r ∗
  i=1 is a set of leading eigenvectors of (−G∗) which corresponds to the eigenvalue λ1(−G∗) = −λn(G∗). Furthermore, there exists at least one such subgradient.
  One motivation for assuming (standard) strict complementarity (Assumption 1 with parameters r = rank(X∗) and δ > 0) is that it guarantees a certain notion of robustness of the problem to small perturbations in the parameters. It is well known (see for instance [3]) that a projected subgradient step from X∗ with respect to a subgradient G∗ ∈ ∂g(X∗) for which the first-order optimality condition holds, returns the optimal solution X∗ itself. This implies that rank (ΠSn [X
  ∗ − ηG∗]) = rank(X∗) (here η is the step-size). Without (standard) strict complementarity however, a small change in the parameters could result in a higher rank matrix. This is captured in the following lemma which is analogous to Lemma 3 in [18], where again the proof is straightforward from the proof in [18] by replacing the gradient of the optimal solution with a subgradient for which the first-order optimality condition holds. Lemma 4. Let X∗ be an optimal solution of rank r∗ to Problem (1). Let G∗ ∈ ∂g(X∗) be a subgradient at X∗ such that 〈X −X∗,G∗〉 ≥ 0 for all X ∈ Sn. Then, λn−r∗(G∗) = λn(G∗) if and only if for any arbitrarily small ζ > 0 it holds that rank ( Π(1+ζ)Sn [X ∗ − ηG∗] ) > r∗, where η > 0, (1 + ζ)Sn = {(1 + ζ)X |X ∈ Sn}, and Π(1+ζ)Sn [·] denotes the Euclidean projection onto the set (1 + ζ)Sn.
sections.heading: 2.1 The challenge of applying low-rank projected subgradient steps
sections.text: We now demonstrate the difficulty of replacing the full-rank SVD computations required in projected subgradient steps over the spectrahedron, with their low-rank SVD counterparts when attempting to solve Problem (1). We prove that a projected subgradient step from a point arbitrarily close to a
  3Denote q(Z, s) = minX∈Sn{g(X) + s(1− Tr(X))− 〈Z,X〉}. The dual problem of Problem (1) can be written as: max{Z 0, s∈R}{q(Z, s) | (Z, s) ∈ dom(q)}.
  low-rank optimal solution — even one that satisfies strict complementarity (Definition 1), may result in a higher rank matrix. The problem on which we demonstrate this phenomena is a well known convex formulation of the sparse PCA problem [12]. The proof is given in Appendix A.2. Lemma 5 (failure of low-rank subgradient descent on sparse PCA). Consider the problem minX∈Sn{g(X) := − 〈 zz> + z⊥z > ⊥,X 〉 + 12k‖X‖1}, where z = (1/ √ k, . . . , 1/ √ k, 0, . . . , 0)> is
  supported on the first k entries, z⊥ = (0, . . . , 0, 1/ √ n− k, . . . , 1/ √ n− k)> is supported on the last n−k entries, and k ≤ n/4. Then, zz> is a rank-one optimal solution for which strict complementarity holds. However, for any η < 23 and any v ∈ R
  n such that ‖v‖ = 1, support(v) ⊆ support(z), and 〈z,v〉2 = 1 − 12‖vv > − zz>‖2F ≥ 1 − 12k2 , it holds that rank ( ΠSn [vv > − ηGvv> ] ) > 1, where Gvv> = −zz> − z⊥z>⊥ + 12k sign(vv >) ∈ ∂g(vv>).
  Note that the subgradient of the `1-norm which we choose for the projected subgradient step simply corresponds to the sign function, which is arguably the most natural choice.
sections.heading: 3 From nonsmooth to saddle-point problems
sections.text: To circumvent the difficulty demonstrated in Lemma 5 in incorporating low-rank SVDs into standard subgradient methods for solving Problem (1), we propose tackling the nonsmooth problem with saddle-point methods.
  We assume the nonsmooth Problem (1) can be written as a maximum of smooth functions, i.e., g(X) = maxy∈K f(X,y), where K ⊂ Y is some compact and convex subset of the finite linear space over the reals Y onto which it is efficient to compute Euclidean projections. We assume f(·,y) is convex for all y ∈ K and f(X, ·) is concave for all X ∈ Sn. That is, we rewrite Problem (1) as the following equivalent saddle-point problem:
  min X∈Sn max y∈K f(X,y). (4)
  Finding an optimal solution to problem (4) is equivalent to finding a saddle-point (X∗,y∗) ∈ Sn ×K such that for all X ∈ Sn and y ∈ K, f(X∗,y) ≤ f(X∗,y∗) ≤ f(X,y∗). We make a standard assumption that f(·, ·) is smooth with respect to all components. That is, we assume there exist βX , βy, βXy, βyX ≥ 0 such that for any X, X̃ ∈ Sn and y, ỹ ∈ K the following four inequalities hold: ‖∇Xf(X,y)−∇Xf(X̃,y)‖F ≤ βX‖X− X̃‖F , ‖∇yf(X,y)− ∇yf(X, ỹ)‖2 ≤ βy‖y − ỹ‖2, ‖∇Xf(X,y) − ∇Xf(X, ỹ)‖F ≤ βXy‖y − ỹ‖2, ‖∇yf(X,y) − ∇yf(X̃,y)‖2 ≤ βyX‖X− X̃‖F , where∇Xf = ∂f∂X and ∇yf = ∂f ∂y .
  We denote by β the full Lipschitz parameter of the gradient, that is for any X, X̃ ∈ Sn and y, ỹ ∈ K,
  ‖(∇Xf(X,y),−∇yf(X,y))− (∇Xf(X̃, ỹ),−∇yf(X̃, ỹ))‖ ≤ β‖(X,Y)− (X̃, ỹ)‖, where ‖ · ‖ denotes the Euclidean norm over the product space Sn × Y. The relationship between the full Lipschitz parameter β and its components βX , βy, βXy, βyX can be written as β = √ 2 max {√
  β2X + β 2 yX , √ β2y + β 2 Xy } . See proof in Appendix B.
  The following lemma highlights a connection between the gradient of a saddle-point of (4) and subgradients of an optimal solution to (1) for which the first order optimality condition holds. One of the connections we will be interested in, is that GSC for Problem (1) implies GSC (with the same parameters) for Problem (4). However, to prove this specific connection we require an additional structural assumption on the objective function g(·). We note that this assumption holds for all applications mentioned in this paper. Assumption 2. g(X) is of the form g(X) = h(X) + maxy∈K y>(A(X)−b), where h(·) is smooth and convex, and A is a linear map. Lemma 6. If (X∗,y∗) is a saddle-point of Problem (4) then X∗ is an optimal solution to Problem (1), ∇Xf(X∗,y∗) ∈ ∂g(X∗), and for all X ∈ Sn it holds that 〈X − X∗,∇Xf(X∗,y∗)〉 ≥ 0. Conversely, under Assumption 2, if X∗ is an optimal solution to Problem (1), and G∗ ∈ ∂g(X∗) which satisfies 〈X −X∗,G∗〉 ≥ 0 for all X ∈ Sn, then there exists y∗ ∈ arg maxy∈K f(X∗,y) such that (X∗,y∗) is a saddle-point of Problem (4), and ∇Xf(X∗,y∗) = G∗.
  The proof is given in Appendix C.1. The connection between the gradient of an optimal solution to the saddle-point problem and a subgradient of a corresponding optimal solution in the equivalent nonsmooth problem established in Lemma 6, naturally leads to the formulation of the following generalized strict complementarity assumption for saddle-point problems.
  Assumption 3 (generalized strict complementarity for saddle-points). We say a saddle-point (X∗,y∗) ∈ Sn×K of Problem (4) with rank(X∗) = r∗ satisfies the generalized strict complementarity assumption with parameters r ≥ r∗, δ > 0, if λn−r(∇Xf(X∗,y∗))− λn(∇Xf(X∗,y∗)) ≥ δ. Remark 1. Note that under Assumption 2, due to Lemma 6, GSC with parameters r, δ for some optimal solution X∗ to Problem (1) implies GSC with parameters r, δ to a corresponding saddle-point (X∗,y∗) of Problem (4). Nevertheless, Assumption 2 is not necessary for proving our convergence results for Problem (4), which are directly stated in terms of Assumption 3.
sections.heading: 4 Projected extragradient method with low-rank projections
sections.text: In this section we formally state and prove our main result: the projected extragradient method for the saddle-point Problem (4), when initialized in the proximity of a saddle-point which satisfies GSC (Assumption 3), converges with its original O(1/t) rate while requiring only two low-rank SVD computations per iteration.
  Algorithm 1 Projected extragradient method for saddle-point problems (see also [23, 30]) Input: sequence of step-sizes {ηt}t≥1 Initialization: (X1,y1) ∈ Sn ×K for t = 1, 2, ... do Zt+1 = ΠSn [Xt − ηt∇Xf(Xt,yt)] wt+1 = ΠK[yt + ηt∇yf(Xt,yt)] Xt+1 = ΠSn [Xt − ηt∇Xf(Zt+1,wt+1)] yt+1 = ΠK[yt + ηt∇yf(Zt+1,wt+1)]
  end for
  Theorem 1 (main theorem). Fix an optimal solution (X∗,y∗) ∈ Sn ×K to Problem (4). Let r̃ denote the multiplicity of λn(∇Xf(X∗,y∗)) and for any r ≥ r̃ define δ(r) = λn−r(∇Xf(X∗,y∗)− λn(∇Xf(X∗,y∗). Let {(Xt,yt)}t≥1 and {(Zt,wt)}t≥2 be the sequences of iterates generated by Algorithm 1 with a fixed step-size η = min { 1
  2 √ β2X+β 2 yX , 1 2 √ β2y+β 2 Xy , 1βX+βXy , 1 βy+βyX
  } . Assume
  the initialization (X1,y1) satisfies ‖(X1,y1)− (X∗,y∗)‖ ≤ R0(r), where
  R0(r) := η
  (1 + √ 2) ( 1 + (2 + √ 2)ηmax{βX , βXy} ) max{√r̃δ(r − r̃ + 1) 2 , δ(r) (1 + 1/ √ r̃) } .
  Then, for all t ≥ 1, the projections ΠSn [Xt − η∇Xf(Xt,yt)] and ΠSn [Xt − η∇Xf(Zt+1,wt+1)] can be replaced with their rank-r truncated counterparts (see (3)) without changing the sequences {(Xt,yt)}t≥1 and {(Zt,wt)}t≥2, and for any T ≥ 0 it holds that
  1
  T T∑ t=1 max y∈K f(Zt+1,y)− 1 T T∑ t=1 min X∈Sn f(X,wt+1)
  ≤ D2 max
  {√ β2X + β 2 yX , √ β2y + β 2 Xy, 1 2 (βX + βXy), 1 2 (βy + βyX) } T ,
  where D := sup(X,y),(Z,w)∈Sn×K ‖(X,y)− (Z,w)‖. Remark 2. Note that Theorem 1 implies that if standard strict complementarity holds for Problem (4), that is Assumption 3 holds with r = r∗ = rank(X∗) and some δ > 0, then only rank-r∗ SVDs are required so that Algorithm 1 converges with the guaranteed convergence rate of O(1/t), when initialized with a “warm-start”. Furthermore, by using SVDs of rank r > r∗, with moderately higher values of r, we can increase the radius of the ball in which Algorithm 1 needs to be initialized quite significantly.
  The complete proof of Theorem 1 is given in Appendix D. Below we give a short sketch of the proof.
  Proof sketch for Theorem 1. The convergence rate of Algorithm 1 is well known in the literature (see Appendix D.1 for completeness). Thus, we focus on the novel part of the low-rank projections. For simplicity consider the case r = r̃. From (2) it can be deduced that for any P ∈ Sn it holds that rank(ΠSn [P]) ≤ r̃ if and only if the condition ∑r̃ i=1 λi(P) ≥ 1 + r̃λr̃+1(P) holds. Denote P∗ = X∗ − η∇Xf(X∗,y∗). From Lemma 6 and Lemma 3 it follows that ∀i ≤ rank(X∗) : λi(P∗) = λi(X∗)− ηλn(∇Xf(X∗,y∗)); ∀i > rank(X∗) : λi(P∗) = −ηλn−i+1(∇Xf(X∗,y∗)).
  Using this and the fact that λn−i+1(∇fX(X∗,y∗)) = λn(∇fX(X∗,y∗)) for all i ≤ r̃ we have, r̃∑ i=1 λi(P ∗) = rank(X∗)∑ i=1 λi(X ∗)− η r̃∑ i=1 λn(∇Xf(X∗,y∗)) = Tr(X∗)− ηr̃λn(∇Xf(X∗,y∗)) = 1− ηr̃λn−r̃(∇Xf(X∗,y∗)) + ηr̃δ(r̃) = 1 + r̃λr̃+1(P∗) + ηr̃δ(r̃). (5) Thus, P∗ = X∗ − η∇Xf(X∗,y∗) not only satisfies the condition that ensures its projection onto Sn is of rank at most r̃, but it satisfies it with a positive slack of ηr̃δ(r̃). Also, for any (X,y) sufficiently close to (X∗,y∗), using the smoothness of f , we have that P = X− η∇Xf(X,y) is close to P∗. Thus, by applying perturbation bounds for the eigenvalues of symmetric matrices to Eq. (5) and using the fact that the positive slack ηr̃δ(r̃) in the RHS of (5) allows to absorb sufficiently small errors, we can establish that for such points (X,y) it holds that rank(ΠSn [P]) ≤ r̃. A similar argument shows that the second primal step in Algorithm 1: Xt+1 = ΠSn [Xt − η∇Xf(Zt+1,wt+1)], also results in matrices of rank at most r̃ when Zt,wt are also in the proximity of X∗,y∗. Finally, we prove a complementary argument, that when initialized in the proximity of (X∗,y∗), the iterates of Algorithm 1 stay in the proximity of (X∗,y∗) throughout all iterations.
  Remark 3. Note that when applying Algorithm 1 towards minimizing a nonsmooth objective of the form g(X) = maxy∈K f(X,y) (as discussed in Section 3), a guarantee of the form 1 T ∑T t=1 maxy∈K f(Zt+1,y) − 1 T ∑T t=1 minX∈Sn f(X,wt+1) ≤ , for some > 0, which is what we get from Theorem 1, implies in particular that mint∈[T ] g(Xt)− g(X∗) ≤ , where X∗ is a minimizer of g(·) over Sn. Remark 4. A downside of considering the saddle-point formulation (4) when attempting to solve Problem (1) that arises from Theorem 1, is that not only do we need a “warm-start” initialization for the original primal matrix variable X, in the saddle-point formulation we need a “warm-start” for the primal-dual pair (X,y). Nevertheless, as we demonstrate extensively in Section 5, it seems that very simple initialization schemes work very well in practice.
sections.heading: 4.1 Efficiently-computable certificates for correctness of low-rank projections
sections.text: Since Theorem 1 only applies in some neighborhood of an optimal solution, it is of interest to have a procedure for verifying if the rank-r truncated projection of a given point indeed equals the exact Euclidean projection. In particular, from a practical point of view, it does not matter whether the conditions of Theorem 1 hold. In practice, as long as the truncated projection (see (3)) equals the exact projection (see (2)), we are guaranteed that Algorithm 1 converges correctly with rate O(1/t), without needing to verify any other condition. Luckily, the expression in (2) which characterizes the structure of the Euclidean projection onto the spectrahedron, yields exactly such a verification procedure. As already noted in [17], for any X ∈ Sn, we have Π̂rSn [X] = ΠSn [X] if an only if the condition ∑r i=1 λi(X) ≥ 1 + r · λr+1(X) holds. Note that verifying this condition simply requires increasing the rank of the SVD computation by one, i.e., computing a rank-(r+ 1) SVD of the matrix to project rather than a rank-r SVD.
sections.heading: 5 Empirical evidence
sections.text: The goal of this section is to bring empirical evidence in support of our theoretical approach. We consider various tasks that take the form of minimizing a composite objective, i.e., the sum of a
  smooth convex function and a nonsmooth convex function, where the nonsmoothness comes from either an `1-norm or `2-norm regularizer / penalty term, over a τ -scaled spectrahedron. In all cases the nonsmooth objective can be written as a saddle-point with function f(X,y) which is linear in y and in particular satisfies Assumption 2.
  The tasks considered include 1. sparse PCA, 2. robust PCA, 3. low-rank and sparse recovery, 4. phase synchronization, and 5. linearly-constrained low-rank estimation, under variety of parameters. Due to lack of space many of the results are deferred to Appendix F.
  For all tasks considered we generate random instances, and examine the sequences of iterates generated by Algorithm 1 {(Xt,yt)}t≥1, {(Zt,wt)}t≥2, when initialized with simple initialization procedures. Out of both sequences generated, we choose our candidate for the optimal solution to be the iterate for which the dual-gap, which is a certificate for optimality, is smallest. See Appendix E.
  In all tasks considered the goal is to recover a ground-truth low-rank matrix M0 ∈ Sn from some noisy observation of it M = M0 + N, where N ∈ Sn is a noise matrix. We measure the signal-to-noise ratio (SNR) as ‖M0‖2F / ‖N‖2F . In all experiments we measure the relative
  initialization error by ∥∥∥Tr(M0)τ X1 −M0∥∥∥2
  F
  / ‖M0‖2F , and similarly we measure the relative re-
  covery error by ∥∥∥Tr(M0)τ X∗ −M0∥∥∥2
  F
  / ‖M0‖2F . Note that in some of the experiments we take
  τ < Tr(M0) to prevent the method from overfitting the noise. In addition, we measure the (standard) strict complementarity parameter which corresponds to the eigen-gap gap(∇Xf(X∗,y∗)) := λn−r(∇Xf(X∗,y∗))− λn(∇Xf(X∗,y∗)), r = rank(M0). In all experiments we use SVDs of rank r = rank(M0) to compute the projections in Algorithm 1 according to the truncated projection given in (3). To certify the correctness of these low-rank projections (that is, that they equal the exact Euclidean projection) we confirm that the inequality∑r i=1 λi(Pj) ≥ τ + r · λr+1(Pj) always holds for P1 = Xt − η∇Xf(Xt,Yt) and P2 = Xt − η∇Xf(Zt+1,Wt+1) (see also Section 4.1). Indeed, we can now already state our main observation from the experiments:
  In all tasks considered and for all random instances generated, throughout all iterations of Algorithm 1, when initialized with a simple “warm-start” strategy and when computing only rank-r truncated projections, r = rank(M0), the truncated projections of P1 = Xt − η∇Xf(Xt,Yt) and P2 = Xt − η∇Xf(Zt+1,Wt+1) equal their exact full-rank counterparts. That is, Algorithm 1, using only rank-r SVDs, computed exactly the same sequences of iterates it would have computed if using full-rank SVDs.
  Aside from the above observation, in the sequel we demonstrate that all models considered indeed satisfy that: 1. the returned solution, denoted (X∗,y∗), is of the same rank as the ground-truth matrix and satisfies the strict complementarity condition with non-negligible parameter (measured by the eigengap λn−r(∇Xf(X∗,y∗))− λn(∇Xf(X∗,y∗))), 2. the recovery error of the returned solution indeed improves significantly over the error of the initialization point.
  Sparse PCA: We consider the sparse PCA problem in a well known convex formulation taken from [12] and its equivalent saddle-point formulation:
  min Tr(X)=1,
  X 0
  〈X,−M〉+ λ‖X‖1 = min Tr(X)=1,
  X 0
  max ‖Y‖∞≤1
  {〈X,−M〉+ λ〈X,Y〉},
  where M = zz> + c2 (N + N >) is a noisy observation of a rank-one matrix zz>, with z being a sparse unit vector. Each entry zi is chosen to be 0 with probability 0.9 and U{1, . . . , 10} with probability 0.1, and then we normalize z to be of unit norm. We test the results obtained when adding different magnitudes of Gaussian or uniform noise. We set the signal-to-noise ratio (SNR) to be a constant. Thus, we set the noise level to c = 2SNR·‖N+N>‖F for our choice of SNR. We initialize the X variable with the rank-one approximation of M. That is, we take X1 = u1u>1 , where u1 is the top eigenvector of M. For the Y variable we initialize it with Y1 = sign(X1) which is a subgradient of ‖X1‖1. We set the step-size to η = 1/(2λ) and we set the number of iterations to T = 1000 and for any set of parameters we average the measurements over 10 i.i.d. runs.
  Robust PCA: We consider the robust PCA problem [28] in the following formulation: min
  Tr(X)=τ, X 0
  ‖X−M‖1 = min Tr(X)=τ,
  X 0
  max ‖Y‖∞≤1
  〈X−M,Y〉,
  where M = rZ0Z>0 + 1 2 (N+N >) is a sparsely-corrupted observation of some rank-r matrix Z0Z>0 . We choose Z0 ∈ Rn×r to be a random unit Frobenius norm matrix. For N ∈ Rn×n, we choose each entry to be 0 with probability 1 − 1/ √ n and otherwise 1 or −1 with equal probability. We initialize the X variable with the projection X1 = Π{Tr(X)=τ, X 0}[M], and the Y variable with Y1 = sign(X1 −M). We vary the rank of Z0, which in turn determines the step-size η. We set the trace bound to τ = 0.95 · Tr(rZ0Z>0 ), and the number of iterations to T = 30, 000. For every set of parameters we average the measurements over 10 i.i.d. runs.
  Linearly-constrained low-rank estimation: Consider the following penalized formulation: min
  Tr(X)=1, X 0
  〈X,−M〉+ λ‖A(X)− b‖2 = min Tr(X)=1,
  X 0
  max ‖y‖2≤1
  〈X,−M〉+ λ〈A(X)− b,y〉,
  where M = z0z>0 + c 2 (N + N >) is a noisy observation of some rank-one matrix z0z>0 such that ‖z0‖2 = 1 and the noise matrix is chosen N ∼ N (0, In). We take A(X) = (〈A1,X〉, . . . , 〈Am,X〉)> with matrices A1, . . . ,Am ∈ Sn of the form Ai = viv>i such that vi ∼ N (0, 1). We take b ∈ Rm such that bi = 〈Ai, z0z>0 〉. We initialize the X variable with the rank-one approximation of M. That is, we take X1 = u1u>1 , where u1 is the top eigenvector of M. The y variable is initialized with y1 = (A(X1) − b)/‖A(X1) − b‖2. We set the number of constraints to m = n, the penalty parameter to λ = 2, and the step-size to η = 1/(2λ). We set the number of iterations in each experiment to T = 2000 and for each value of n we average the measurements over 10 i.i.d. runs.
  Low-rank and sparse matrix recovery: We consider the problem of recovering a simultaneously low-rank and sparse covariance matrix [35], which can be written as the following saddle-point optimization problem:
  min Tr(X)=1,
  X 0
  1 2 ‖X−M‖2F + λ‖X‖1 = min
  Tr(X)=τ, X 0
  max ‖Y‖∞≤1
  1 2 ‖X−M‖2F + λ〈X,Y〉,
  where M = Z0Z0> + c2 (N + N >) is a noisy observation of some low-rank and sparse covariance matrix Z0Z0>. We choose Z0 ∈ Rn×r to be a sparse matrix where each entry Z0i,j is chosen to be
  0 with probability 0.9 and U{1, . . . , 10} with probability 0.1, and then we normalize Z0 to be of unit Frobenius norm. We choose N ∼ N (0.5, In). We test the model with rank(Z0Z>0 ) = 5, 10. We set the signal-to-noise ratio (SNR) to be a constant and set the noise level to c = 2‖Z0Z > 0 ‖F
  SNR·‖N+N>‖F for our choice of SNR. We initialize the X variable with the rank-r approximation of M. That is, we take X1 = Urdiag ( Π∆τ,r [diag(−Λr)] ) U>r , where UrΛrU > r is the rank-r eigen-decomposition of M
  and ∆τ,r = {z ∈ Rr | z ≥ 0, ∑r i=1 zi = τ} is the simplex of radius τ in Rr. For the Y variable we initialize it with Y1 = sign(X1) which is a subgradient of ‖X1‖1. We set the step-size to η = 1, τ = 0.7 · Tr(Z0Z>0 ), and the number of iterations in each experiment to T = 2000. For each value of r and n we average the measurements over 10 i.i.d. runs.
sections.heading: 6 Discussion
sections.text: This work expands upon a line of research that aims to harness the ability of convex relaxations to produce low-rank and high-quality solutions to important low-rank matrix optimization problems, while insisting on methods that, at least locally, store and manipulate only low-rank matrices. Focusing on the challenging case of nonsmooth objective functions and following our evidence for the difficulties of obtaining such a result for subgradient methods (Lemma 5), we consider tackling nonsmooth objectives via saddle-point formulations. We prove that indeed under a generalized strict complementarity condition, a state-of-the-art method for convex-concave saddle-point problems converges locally while storing and manipulating only low-rank matrices. Extensive experiments over several tasks demonstrate that our conceptual approach of utilizing low-rank projections for more efficient optimization is not only of theoretical merit, but indeed seems to work well in practice.
sections.heading: Acknowledgements
sections.text: This research was supported by the ISRAEL SCIENCE FOUNDATION (grant No. 1108/18).
title: Low-Rank Extragradient Method for Nonsmooth and Low-Rank Matrix Optimization Problems