
\section{Introduction}

Low-rank and nonsmooth matrix optimization problems have many important applications in statistics, machine learning, and related fields, such as \textit{sparse PCA} \cite{sparsePCA1, sparsePCA2}, \textit{robust PCA} \cite{robustPCA1, robustPCA2, robustPCA3, robustPCA4, robustPCA5},  \textit{phase synchronization} \cite{phaseSyncronization1, phaseSyncronization2, phaseSyncronization3}, \textit{community detection and stochastic block models} \cite{abbe2017community}\footnote{in \cite{phaseSyncronization1, phaseSyncronization2, phaseSyncronization3} and \cite{abbe2017community} the authors consider SDPs with linear objective function and affine constraints of the form $\mathcal{A}({\mathbf{X}})={\mathbf{b}}$. By incorporating the linear constraints into the objective function via a $\ell_2$ penalty term of the form $\lambda\Vert{\mathcal{A}({\mathbf{X}})-{\mathbf{b}}}\Vert_2$, $\lambda > 0$, we obtain a nonsmooth objective function.},  \textit{low-rank and sparse covariance matrix recovery} \cite{lowRankAndSparse}, \textit{robust matrix completion} \cite{klopp2017robust, chen2011robust}, 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:
\begin{align} \label{nonSmoothProblem}
\min_{{\mathbf{X}}\in{\mathcal{S}_n}} g({\mathbf{X}}),
\end{align} 
where $g:\mathbb{S}^n\rightarrow\mathbb{R}$ is convex but nonsmooth, and $\mathcal{S}_n=\{{\mathbf{X}}\in\mathbb{S}^n\ \vert\ \textrm{Tr}({\mathbf{X}})=1,\ {\mathbf{X}}\succeq0\}$ is the spectrahedron in $\mathbb{S}^n$, $\mathbb{S}^n$ being the space of $n\times n$ real symmetric matrices.

Problem \eqref{nonSmoothProblem}, despite being convex, is notoriously difficult to solve in large scale. The simplest and most general approach applicable to it is the \textit{projected subgradient method}  \cite{beckOptimizationBook, sebastienMD}, which requires on each iteration to compute a Euclidean projection onto the spectrahedron $\mathcal{S}_n$, which in worst case amounts to $O(n^3)$ runtime per iteration. In many applications $g({\mathbf{X}})$ follows a composite model, i.e., $g({\mathbf{X}}) = h({\mathbf{X}}) + w({\mathbf{X}})$, where $h(\cdot)$ is convex and smooth and $w(\cdot)$ 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 \cite{FISTA} or splitting methods such as ADMM  \cite{ADMM} 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(n^3)$ runtime. 
A third type of off-the-shelf methods include those which are based on the \textit{conditional gradient method} and adapted to nonsmooth problems, see for instance \cite{odor2016frank, lowRankAndSparseOurs, NEURIPS2020_8f468c87, locatello2019stochastic}. 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 \eqref{nonSmoothProblem} require in worst case $\Omega(n^3)$ runtime per iteration and/or to store high-rank matrices in memory.

In the recent works \cite{garberLowRankSmooth, garberStochasticLowRank} 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\times n$ matrix using efficient iterative methods typically scales with $rn^2$ (and further improves when the matrix is sparse), instead of $n^3$ 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 \textit{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 ${\mathbf{X}}\in\mathbb{S}^n$, which admits an eigen-decomposition ${\mathbf{X}}=\sum_{i=1}^n\lambda_i{\mathbf{v}}_i{\mathbf{v}}_i^{\top}$, is given by 
\begin{align}\label{eq:euclidProj}
\Pi_{\mathcal{S}_n}[{\mathbf{X}}]=\sum_{i=1}^n\max\{0,\lambda_i-\lambda\}{\mathbf{v}}_i{\mathbf{v}}_i^{\top},
\end{align}
where $\lambda\in\mathbb{R}$ is the unique scalar satisfying  $\sum_{i=1}^n\max\{0,\lambda_i-\lambda\}=1$. This operation thus truncates all eigenvalues that are smaller than $\lambda$, while leaving the eigenvectors unchanged, thereby returning a matrix with rank equal to the number of eigenvalues greater than $\lambda$. Importantly, when the projection of ${\mathbf{X}}$ onto $\mathcal{S}_n$ is of rank $r$, only the first $r$ components in the eigen-decomposition of ${\mathbf{X}}$ are required to compute it in the first place, and thus, only a rank-$r$ SVD of ${\mathbf{X}}$ is required.
In other words and simplifying, \cite{garberLowRankSmooth, garberStochasticLowRank}  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:
\begin{align} \label{truncatedProjection}
\widehat{\Pi}_{\mathcal{S}_n}^r[{\mathbf{X}}]:=\Pi_{\mathcal{S}_n}\left[\sum_{i=1}^r\lambda_i{\mathbf{v}}_i{\mathbf{v}}_i^{\top}\right].
\end{align}

Extending the results of \cite{garberLowRankSmooth, garberStochasticLowRank} to the nonsmooth setting is difficult since the smoothness assumption is critical to the analysis. 
Moreover, while  \cite{garberLowRankSmooth, garberStochasticLowRank} 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 \eqref{nonSmoothProblem} 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 \cite{garberLowRankSmooth, garberStochasticLowRank}. Concretely, we show that if a generalized strict complementarity (GSC) assumption holds for a low-rank optimal solution (see Assumption \ref{ass:strictcompNonSmooth} in the sequel), the \textit{extragradient method} for smooth convex-concave saddle-point problems \cite{extragradientK,Nemirovski} (see Algorithm \ref{alg:EG} 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 iteration\footnote{note that the extradgradient method computes two projected-gradient steps on each iteration, and thus two SVDs are needed per iteration.}. 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/\sqrt{t})$ rate \cite{sebastienMD}.

Our contributions can be summarized  as follows:
\begin{itemize}

\item
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 \ref{lemma:negativeExampleNonSmooth}.

\item
We suggest a generalized strict complementarity (GSC) condition for saddle-point problems and prove that when $g(\cdot)$ --- the objective function in Problem \eqref{nonSmoothProblem}, 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  \eqref{nonSmoothProblem} 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 \ref{sec:smooth2Saddle}.


\item
\textbf{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 \ref{theroem:allPutTogether}.

\item 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 \eqref{eq:euclidProj}). See Section \ref{sec:expr}.
\end{itemize}



\subsection{Additional related work}
Since, as in the works  \cite{garberLowRankSmooth, garberStochasticLowRank} 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 \cite{garber2019linear, spectralFrankWolfe, ding2020kfw, garber2020efficient} 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 \cite{zhou2017unified, drusvyatskiy2018error}.

Besides convex relaxations such as Problem \eqref{nonSmoothProblem}, considerable advances have been made in the past several yeas in developing efficient \textit{nonconvex} methods with global convergence guarantees for low-rank matrix problems. 
In \cite{nonConvexFactorizedAssumptions} 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 \cite{nonconvexFactorizedGlobal} 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, \cite{charisopoulos2021low} 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 low-rank matrices from random linear measurements was presented in \cite{li2020nonconvex}.  
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 \cite{lu2014generalized, yao2018large}. 



\section{Strict Complementarity for Nonsmooth Optimization and Difficulty of Applying Low-Rank Projected Subgradient Steps}
\label{sec:strictComp}

Our analysis of the nonsmooth Problem \eqref{nonSmoothProblem} 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:
\begin{lemma}[first-order optimality condition, see \cite{beckOptimizationBook}] \label{optCondition}
Let $g:\mathbb{S}^n\rightarrow\mathbb{R}$ be a convex function. Then ${\mathbf{X}}^*\in{\mathcal{S}_n}$ minimizes $g$ over ${\mathcal{S}_n}$ if and only if there exists a subgradient ${\mathbf{G}}^*\in\partial g({\mathbf{X}}^*)$ such that $\langle {\mathbf{X}}-{\mathbf{X}}^*,{\mathbf{G}}^*\rangle\ge0$ for all ${\mathbf{X}}\in{\mathcal{S}_n}$.
\end{lemma}


For some ${\mathbf{G}}^*\in\partial g({\mathbf{X}}^*)$ which satisfies the first-order optimality condition for an optimal solution ${\mathbf{X}}^*$, if the multiplicity of the smallest eigenvalue equals $r^*=\textrm{rank}({\mathbf{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 \textit{smooth} optimization problem over the spectrahedron and an eigen-gap in the gradient of the optimal solution was established in \cite{spectralFrankWolfe}. We generalize this equivalence to also include nonsmooth problems. The proof follows similar arguments and is given in \cref{sec:proofKKT}.

\begin{definition}[strict complementarity] \label{def:strict}
An optimal solution ${\mathbf{X}}^*\in{\mathcal{S}_n}$ of rank $r^*$ for Problem \eqref{nonSmoothProblem} satisfies the strict complementarity assumption with parameter $\delta >0$, if there exists an optimal solution of the dual problem\footnote{Denote $q({\mathbf{Z}},s)=\min_{{\mathbf{X}}\in\mathbb{S}^n} \lbrace g({\mathbf{X}})+s(1-\textrm{Tr}({\mathbf{X}}))-\langle{\mathbf{Z}},{\mathbf{X}}\rangle\rbrace$.
The dual problem of Problem \eqref{nonSmoothProblem} can be written as:
$ \max_{\lbrace{\mathbf{Z}}\succeq0,\ s\in\mathbb{R}\rbrace} \lbrace q({\mathbf{Z}},s) \ \vert \ ({\mathbf{Z}},s)\in \textrm{dom}(q)\rbrace$.}
 $({\mathbf{Z}}^*,s^*)\in\mathbb{S}^n\times\mathbb{R}$ such that $\textrm{rank}({\mathbf{Z}}^*)=n-r^*$, and $\lambda_{n-r^*}({\mathbf{Z}}^*)\ge\delta$. 
\end{definition}

\begin{lemma} \label{lemma:kkt}
Let ${\mathbf{X}}^*\in{\mathcal{S}_n}$ be a rank-$r^*$ optimal solution to Problem \eqref{nonSmoothProblem}. ${\mathbf{X}}^*$ satisfies the (standard) strict complementarity assumption with parameter $\delta>0$ if and only if there exists a subgradient ${\mathbf{G}}^*\in\partial g({\mathbf{X}}^*)$ such that $\langle {\mathbf{X}}-{\mathbf{X}}^*,{\mathbf{G}}^*\rangle\ge0$ for all ${\mathbf{X}}\in\mathcal{S}_n$ and $\lambda_{n-r^*}({\mathbf{G}}^*)-\lambda_{n}({\mathbf{G}}^*)\ge\delta$.
\end{lemma}


Throughout this paper we assume a weaker and more general assumption than strict complementarity, namely generalized strict complementarity (GSC), which we present now.

\begin{assumption}[generalized strict complementarity]\label{ass:strictcompNonSmooth}
We say an optimal solution ${\mathbf{X}}^*$ to Problem \eqref{nonSmoothProblem} satisfies the generalized strict complementarity assumption with parameters $r,\delta$, if there exists a subgradient ${\mathbf{G}}^*\in\partial g({\mathbf{X}}^*)$ such that $\langle {\mathbf{X}}-{\mathbf{X}}^*,{\mathbf{G}}^*\rangle\ge0$ for all ${\mathbf{X}}\in\mathcal{S}_n$ and $\lambda_{n-r}({\mathbf{G}}^*)-\lambda_{n}({\mathbf{G}}^*) \geq \delta$.
\end{assumption}


In \cite{garberLowRankSmooth} 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 \cref{optCondition}, we can extend these properties also to the nonsmooth setting.
The following lemma shows that GSC with parameters $(r,\delta)$ for some $\delta>0$ (\cref{ass:strictcompNonSmooth}) 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 \cite{garberLowRankSmooth}, by replacing the gradient of the optimal solution with a subgradient for which the first-order optimality condition holds. 
\begin{lemma} \label{lemma:Xopt-subgradEigen}
Let ${\mathbf{X}}^*$ be an optimal solution to Problem \eqref{nonSmoothProblem} and write its eigen-decomposition as ${\mathbf{X}}^*=\sum_{i=1}^{r^*}{\lambda_i{\mathbf{v}}_i{\mathbf{v}}_i^T}$. Then, any subgradient ${\mathbf{G}}^*\in\partial g({\mathbf{X}}^*)$ which satisfies $\langle {\mathbf{X}}-{\mathbf{X}}^*,{\mathbf{G}}^*\rangle\ge0$ for all ${\mathbf{X}}\in\mathcal{S}_n$, admits an eigen-decomposition such that the set of vectors $\{{\mathbf{v}}_i\}_{i=1}^{r^*}$ is a set of leading eigenvectors of $(-{\mathbf{G}}^*)$ which corresponds to the eigenvalue $\lambda_1(-{\mathbf{G}}^*)=-\lambda_n({\mathbf{G}}^*)$. Furthermore, there exists at least one such subgradient. 
\end{lemma}

One motivation for assuming (standard) strict complementarity (\cref{ass:strictcompNonSmooth} with parameters $r=\textrm{rank}({\mathbf{X}}^*)$ and $\delta>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 \cite{beckOptimizationBook}) that a projected subgradient step from ${\mathbf{X}}^*$ with respect to a subgradient ${\mathbf{G}}^*\in\partial g({\mathbf{X}}^*)$ for which the first-order optimality condition holds, returns the optimal solution ${\mathbf{X}}^*$ itself. This implies that $\textrm{rank}\left( \Pi_{{\mathcal{S}_n}} [ {\mathbf{X}}^* - \eta{\mathbf{G}}^* ]\right)=\textrm{rank}({\mathbf{X}}^*)$ (here $\eta$ 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 \cite{garberStochasticLowRank}, where again the proof is straightforward from the proof in  \cite{garberStochasticLowRank} by replacing the gradient of the optimal solution with a subgradient for which the first-order optimality condition holds. 
\begin{lemma}
Let ${\mathbf{X}}^*$ be an optimal solution of rank $r^*$ to Problem \eqref{nonSmoothProblem}. Let ${\mathbf{G}}^*\in\partial g({\mathbf{X}}^*)$ be a subgradient at ${\mathbf{X}}^*$ such that $\langle {\mathbf{X}}-{\mathbf{X}}^*,{\mathbf{G}}^*\rangle\ge0$ for all ${\mathbf{X}}\in\mathcal{S}_n$. Then, $\lambda_{n-r^*}({\mathbf{G}}^*)=\lambda_{n}({\mathbf{G}}^*)$ if and only if for any arbitrarily small $\zeta>0$ it holds that 
$\textrm{rank}\left( \Pi_{(1+\zeta)\mathcal{S}_n} [ {\mathbf{X}}^* - \eta{\mathbf{G}}^* ]\right)>r^*$,
where $\eta>0$, $(1+\zeta)\mathcal{S}_n = \{(1+\zeta){\mathbf{X}} \ \vert\ {\mathbf{X}}\in\mathcal{S}_n\}$, and $\Pi_{(1+\zeta)\mathcal{S}_n}{[\cdot]}$ denotes the Euclidean projection onto the set $(1+\zeta)\mathcal{S}_n$.
\end{lemma}

\subsection{The challenge of applying low-rank projected subgradient steps}

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 \eqref{nonSmoothProblem}. We prove that a projected subgradient step from a point arbitrarily close to a low-rank optimal solution --- even one that satisfies strict complementarity (Definition \ref{def:strict}), may result in a higher rank matrix.  The problem on which we demonstrate this phenomena is a well known convex formulation of the \textit{sparse PCA} problem \cite{d2007direct}.


\begin{lemma}[failure of low-rank subgradient descent on sparse PCA] \label{lemma:negativeExampleNonSmooth}
Consider the problem
\begin{align*}
\min_{{\mathbf{X}}\in{\mathcal{S}_n}}\{g({\mathbf{X}}):=-\left\langle{\mathbf{z}}\z^{\top}+{\mathbf{z}}_{\perp}{\mathbf{z}}_{\perp}^{\top},{\mathbf{X}}\right\rangle+\frac{1}{2k}\Vert{\mathbf{X}}\Vert_1\},
\end{align*}
where ${\mathbf{z}}=(1/\sqrt{k},\ldots,1/\sqrt{k},0,\ldots,0)^{\top}$ is supported on the first $k$ entries,  ${\mathbf{z}}_{\perp}=(0,\ldots,0,\\ 1/\sqrt{n-k},\ldots,1/\sqrt{n-k})^{\top}$ is supported on the last $n-k$ entries, and $k\le n/4$. Then, ${\mathbf{z}}\z^{\top}$ is a rank-one optimal solution for which  strict complementarity holds.
However, for any $\eta<\frac{2}{3}$ and any ${\mathbf{v}}\in\mathbb{R}^n$ such that $\Vert{\mathbf{v}}\Vert=1$, $\textrm{support}({\mathbf{v}})\subseteq\textrm{support}({\mathbf{z}})$, and $\langle{\mathbf{z}},{\mathbf{v}}\rangle^2=1 - \frac{1}{2}\Vert{{\mathbf{v}}\v^{\top}-{\mathbf{z}}\z^{\top}}\Vert_F^2 \ge1-\frac{1}{2k^2}$, it holds that
\begin{align*}
\textrm{rank}\left(\Pi_{{\mathcal{S}_n}}[{\mathbf{v}}\v^{\top}-\eta{\mathbf{G}}_{{\mathbf{v}}\v^{\top}}]\right)>1,
\end{align*}
where ${\mathbf{G}}_{{\mathbf{v}}\v^{\top}}=-{\mathbf{z}}\z^{\top}-{\mathbf{z}}_{\perp}{\mathbf{z}}_{\perp}^{\top}+\frac{1}{2k}\textrm{sign}({\mathbf{v}}\v^{\top})\in\partial g({\mathbf{v}}\v^{\top})$.
\end{lemma}
Note that the subgradient of the $\ell_1$-norm which we choose for the projected subgradient step simply corresponds to the sign function, which is arguably the most natural choice. 


\begin{proof}

${\mathbf{z}}\z^{\top}$ is a rank-one optimal solution for this problem since for the subgradient $k{\mathbf{z}}\z^{\top}+2k{\mathbf{z}}_{\perp}{\mathbf{z}}_{\perp}^{\top}\in\partial \left(\Vert{\mathbf{z}}\z^{\top}\Vert_1\right)$ the first-order optimality condition holds. Indeed, for all ${\mathbf{X}}\in{\mathcal{S}_n}$ 
\begin{align} \label{eq:ExampleFirstOrderOptCondHolds}
\langle{\mathbf{X}}-{\mathbf{z}}\z^{\top},-{\mathbf{z}}\z^{\top}-{\mathbf{z}}_{\perp}{\mathbf{z}}_{\perp}^{\top}+\frac{1}{2}{\mathbf{z}}\z^{\top}+{\mathbf{z}}_{\perp}{\mathbf{z}}_{\perp}^{\top}\rangle = \langle{\mathbf{X}}-{\mathbf{z}}\z^{\top},-\frac{1}{2}{\mathbf{z}}\z^{\top}\rangle\ge0.
\end{align}

For the subgradient $-\frac{1}{2}{\mathbf{z}}\z^{\top}\in\partial g({\mathbf{z}}\z^{\top})$ there is a gap $\lambda_{n-1}(-\frac{1}{2}{\mathbf{z}}\z^{\top})-\lambda_{n}(-\frac{1}{2}{\mathbf{z}}\z^{\top})=\frac{1}{2}>0$, and as we showed in \eqref{eq:ExampleFirstOrderOptCondHolds} the first order optimality condition holds for $-\frac{1}{2}{\mathbf{z}}\z^{\top}$. Thus, by \cref{lemma:kkt} the optimal solution ${\mathbf{z}}\z^{\top}$ satisfies standard strict complementarity. 

We will show that the projection onto the spectrahedron of a subgradient step from ${\mathbf{v}}\v^{\top}$ with respect to the natural subgradient of the $\ell_1$-norm $\textrm{sign}({\mathbf{v}}\v^{\top})\in\partial \left(\Vert{\mathbf{v}}\v^{\top}\Vert_1\right)$ returns a rank-2 solution.

It holds that
\begin{align*}
1-\frac{1}{2k^2}\le\langle{\mathbf{z}}\z^{\top},{\mathbf{v}}\v^{\top}\rangle=\frac{1}{2}\left(\Vert{\mathbf{z}}\z^{\top}\Vert_F^2+\Vert{\mathbf{v}}\v^{\top}\Vert_F^2-\Vert{\mathbf{v}}\v^{\top}-{\mathbf{z}}\z^{\top}\Vert_F^2\right) = 1- \frac{1}{2}\Vert{\mathbf{v}}\v^{\top}-{\mathbf{z}}\z^{\top}\Vert_F^2,
\end{align*}
and equivalently
\begin{align*}
\sum_{i=1}^k\sum_{j=1}^k\left(\frac{1}{k}-({\mathbf{v}}\v^{\top})_{ij}\right)^2=\Vert{\mathbf{v}}\v^{\top}-{\mathbf{z}}\z^{\top}\Vert_F^2\le\frac{1}{k^2}.
\end{align*}
Therefore, for every $i,j\in\lbrace 1,\ldots,k\rbrace$ it holds that
\begin{align*}
\left\vert({\mathbf{v}}\v^{\top})_{ij}-\frac{1}{k}\right\vert\le\frac{1}{k},
\end{align*}
which implies that $0\le({\mathbf{v}}\v^{\top})_{ij}\le\frac{2}{k}$. Therefore, $k{\mathbf{z}}\z^{\top}=\textrm{sign}({\mathbf{v}}\v^{\top})\in\partial \left(\Vert{\mathbf{v}}\v^{\top}\Vert_1\right)$.

Taking a projected subgradient step from ${\mathbf{v}}\v^{\top}$ with respect to the subgradient $-{\mathbf{z}}\z^{\top}-{\mathbf{z}}_{\perp}{\mathbf{z}}_{\perp}^{\top}+\frac{1}{2}{\mathbf{z}}\z^{\top}\in\partial g({\mathbf{v}}\v^{\top})$ has the form
\begin{align*}
&\Pi_{{\mathcal{S}_n}}\left[{\mathbf{v}}\v^{\top}-\eta\left(-{\mathbf{z}}\z^{\top}-{\mathbf{z}}_{\perp}{\mathbf{z}}_{\perp}^{\top}+\frac{1}{2}{\mathbf{z}}\z^{\top}\right)\right]
 = \Pi_{{\mathcal{S}_n}}\left[{\mathbf{v}}\v^{\top}+\frac{\eta}{2}{\mathbf{z}}\z^{\top}+\eta{\mathbf{z}}_{\perp}{\mathbf{z}}_{\perp}^{\top}\right].
\end{align*}

Since $\textrm{support}({\mathbf{v}})\subseteq\textrm{support}({\mathbf{z}})$ it holds that $\left({\mathbf{v}}\v^{\top}+\frac{\eta}{2}{\mathbf{z}}\z^{\top}\right)\perp{\mathbf{z}}_{\perp}{\mathbf{z}}_{\perp}^{\top}$.
${\mathbf{v}}\v^{\top}+\frac{\eta}{2}{\mathbf{z}}\z^{\top}$ is a rank-2 matrix and so we can denote the eigen-decomposition of ${\mathbf{v}}\v^{\top}+\frac{\eta}{2}{\mathbf{z}}\z^{\top}+\eta{\mathbf{z}}_{\perp}{\mathbf{z}}_{\perp}^{\top}$ as ${\mathbf{v}}\v^{\top}+\frac{\eta}{2}{\mathbf{z}}\z^{\top}+\eta{\mathbf{z}}_{\perp}{\mathbf{z}}_{\perp}^{\top}=\lambda_1{\mathbf{v}}_1{\mathbf{v}}_1^{\top}+\lambda_2{\mathbf{v}}_2{\mathbf{v}}_2^{\top}+\eta{\mathbf{z}}_{\perp}{\mathbf{z}}_{\perp}^{\top}$, where $\lambda_1\ge\lambda_2$. 
Thus, invoking \eqref{eq:euclidProj} to calculate the projection we need to find the scalar $\lambda\in\mathbb{R}$ for which the following holds.
\begin{align*}
\max\left\lbrace\lambda_1-\lambda,0\right\rbrace+\max\left\lbrace\lambda_2-\lambda,0\right\rbrace+\max\lbrace\eta-\lambda,0\rbrace
+\sum_{i=4}^n\max \lbrace0-\lambda,0\rbrace = 1.
\end{align*}

$\lambda_1$ is the largest eigenvalue of ${\mathbf{v}}\v^{\top}+\frac{\eta}{2}{\mathbf{z}}\z^{\top}+\eta{\mathbf{z}}_{\perp}{\mathbf{z}}_{\perp}^{\top}$ since, under our assumption that $\eta < 2/3$, we have that
\[\lambda_1\ge \frac{1}{2}(\lambda_1+\lambda_2)=\frac{1}{2}\textrm{Tr}({\mathbf{v}}\v^{\top}+\frac{\eta}{2}{\mathbf{z}}\z^{\top})=\frac{1}{2}+\frac{\eta}{4}>\eta.\]
Therefore, $\lambda<\lambda_1\le \lambda_1+\lambda_2=\textrm{Tr}({\mathbf{v}}\v^{\top}+\frac{\eta}{2}{\mathbf{z}}\z^{\top})=1+\frac{\eta}{2}$.

In addition, $\max\lbrace\lambda_2,\eta\rbrace\ge\eta$. Therefore, 
\[\lambda_1-\max\lbrace\lambda_2,\eta\rbrace \le \lambda_1 + \lambda_2 - \eta = 1+\frac{\eta}{2}-\eta <1,\]
and so we must have that $\lambda<\max\lbrace\lambda_2,\eta\rbrace\le\lambda_1$. 

This implies that both $\max\left\lbrace\lambda_1-\lambda,0\right\rbrace>0$
and $\max\lbrace\max\lbrace\lambda_2,\eta\rbrace-\lambda,0\rbrace>0$. Thus, using  \eqref{eq:euclidProj} we conclude that 
\[\textrm{rank}\left(\Pi_{{\mathcal{S}_n}}\left[{\mathbf{v}}\v^{\top}-\eta\left(-{\mathbf{z}}\z^{\top}-{\mathbf{z}}_{\perp}{\mathbf{z}}_{\perp}^{\top}+\frac{1}{2}{\mathbf{z}}\z^{\top}\right)\right]\right)\ge2.\]

\end{proof}


\section{From Nonsmooth to Saddle-Point Problems}\label{sec:smooth2Saddle}


To circumvent the difficulty demonstrated in \cref{lemma:negativeExampleNonSmooth} in incorporating low-rank SVDs into standard subgradient methods for solving Problem \eqref{nonSmoothProblem}, we propose tackling the nonsmooth problem with saddle-point methods.


We assume the nonsmooth Problem \eqref{nonSmoothProblem} can be written as a maximum of smooth functions, i.e., $g({\mathbf{X}})=\max_{{\mathbf{y}}\in\mathcal{K}}f({\mathbf{X}},{\mathbf{y}})$, where $\mathcal{K}\subset\mathbb{Y}$ is some compact and convex subset of the finite linear space over the reals $\mathbb{Y}$ onto which it is efficient to compute Euclidean projections. We assume $f(\cdot,{\mathbf{y}})$ is convex for all ${\mathbf{y}}\in\mathcal{K}$ and $f({\mathbf{X}},\cdot)$ is concave for all ${\mathbf{X}}\in{\mathcal{S}_n}$. That is, we rewrite Problem \eqref{nonSmoothProblem} as the following equivalent saddle-point  problem:
\begin{align} \label{problem1}
 \min_{{\mathbf{X}}\in\mathcal{S}_n}\max_{{\mathbf{y}}\in\mathcal{K}}{f({\mathbf{X}},{\mathbf{y}})}.
\end{align}


Finding an optimal solution to problem \eqref{problem1} is equivalent to finding a saddle-point $({\mathbf{X}}^*,{\mathbf{y}}^*)\in{{\mathcal{S}_n}\times\mathcal{K}}$ such that for all ${\mathbf{X}}\in{\mathcal{S}_n}$ and ${\mathbf{y}}\in\mathcal{K}$, 
\begin{align*}
f({\mathbf{X}}^*,{\mathbf{y}}) \le f({\mathbf{X}}^*,{\mathbf{y}}^*) \le f({\mathbf{X}},{\mathbf{y}}^*).
\end{align*}

We make a standard assumption that $f(\cdot,\cdot)$ is smooth with respect to all the components. That is, we assume there exist $\beta_X,\beta_y,\beta_{Xy},\beta_{yX}\ge0$ such that for any ${\mathbf{X}},\tilde{{\mathbf{X}}}\in{\mathcal{S}_n}$ and ${\mathbf{y}},\tilde{{\mathbf{y}}}\in\mathcal{K}$ the following four inequalities hold:
\begin{align*}
& \Vert\nabla_{{\mathbf{X}}}f({\mathbf{X}},{\mathbf{y}})-\nabla_{{\mathbf{X}}}f(\tilde{{\mathbf{X}}},{\mathbf{y}})\Vert_{F} \le\beta_{X}\Vert{\mathbf{X}}-\tilde{{\mathbf{X}}}\Vert_{F}, \\
&\Vert\nabla_{{\mathbf{y}}}f({\mathbf{X}},{\mathbf{y}})-\nabla_{{\mathbf{y}}}f({\mathbf{X}},\tilde{{\mathbf{y}}})\Vert_{2} \le\beta_{y}\Vert{\mathbf{y}}-\tilde{{\mathbf{y}}}\Vert_{2}, \\
& \Vert\nabla_{{\mathbf{X}}}f({\mathbf{X}},{\mathbf{y}})-\nabla_{{\mathbf{X}}}f({\mathbf{X}},\tilde{{\mathbf{y}}})\Vert_{F} \le\beta_{Xy}\Vert{\mathbf{y}}-\tilde{{\mathbf{y}}}\Vert_{2}, \\
& \Vert\nabla_{{\mathbf{y}}}f({\mathbf{X}},{\mathbf{y}})-\nabla_{{\mathbf{y}}}f(\tilde{{\mathbf{X}}},{\mathbf{y}})\Vert_{2} \le\beta_{yX}\Vert{\mathbf{X}}-\tilde{{\mathbf{X}}}\Vert_{F},
\end{align*}
where $\nabla_{{\mathbf{X}}}f=\frac{\partial f}{\partial {\mathbf{X}}}$ and $\nabla_{{\mathbf{y}}}f=\frac{\partial f}{\partial {\mathbf{y}}}$.

We denote by $\beta$ the full Lipschitz parameter of the gradient, that is for any ${\mathbf{X}},\tilde{{\mathbf{X}}}\in{\mathcal{S}_n}$ and ${\mathbf{y}},\tilde{{\mathbf{y}}}\in\mathcal{K}$, 
\begin{align*}
\Vert(\nabla_{{\mathbf{X}}}f({\mathbf{X}},{\mathbf{y}}),-\nabla_{{\mathbf{y}}}f({\mathbf{X}},{\mathbf{y}}))-(\nabla_{{\mathbf{X}}}f(\tilde{{\mathbf{X}}},\tilde{{\mathbf{y}}}),-\nabla_{{\mathbf{y}}}f(\tilde{{\mathbf{X}}},\tilde{{\mathbf{y}}}))\Vert \le\beta\Vert({\mathbf{X}},{\mathbf{Y}})-(\tilde{{\mathbf{X}}},\tilde{{\mathbf{y}}})\Vert,
\end{align*}
where $\Vert\cdot\Vert$ denotes the Euclidean norm over the product space ${\mathcal{S}_n}\times\mathbb{Y}$.


To establish the relationship between $\beta$ and $\beta_X,\beta_y,\beta_{Xy},\beta_{yX}$, we can see that for all ${\mathbf{X}},\tilde{{\mathbf{X}}}\in{\mathcal{S}_n}$ and all ${\mathbf{y}},\tilde{{\mathbf{y}}}\in\mathcal{K}$
\begin{align*}
& \Vert(\nabla_{{\mathbf{X}}}f({\mathbf{X}},{\mathbf{y}}),-\nabla_{{\mathbf{y}}}f({\mathbf{X}},{\mathbf{y}}))-(\nabla_{{\mathbf{X}}}f(\tilde{{\mathbf{X}}},\tilde{{\mathbf{y}}}),-\nabla_{{\mathbf{y}}}f(\tilde{{\mathbf{X}}},\tilde{{\mathbf{y}}}))\Vert^2 
\\ & = \Vert\nabla_{{\mathbf{X}}}f({\mathbf{X}},{\mathbf{y}})-\nabla_{{\mathbf{X}}}f(\tilde{{\mathbf{X}}},{\mathbf{v}})\Vert_F^2 + \Vert\nabla_{{\mathbf{y}}}f({\mathbf{X}},{\mathbf{y}})-\nabla_{{\mathbf{y}}}f(\tilde{{\mathbf{X}}},\tilde{{\mathbf{y}}})\Vert_2^2
\\ & \le 2\Vert\nabla_{{\mathbf{X}}}f({\mathbf{X}},{\mathbf{y}})-\nabla_{{\mathbf{X}}}f(\tilde{{\mathbf{X}}},{\mathbf{y}})\Vert_F^2 + 2\Vert\nabla_{{\mathbf{X}}}f(\tilde{{\mathbf{X}}},{\mathbf{y}})-\nabla_{{\mathbf{X}}}f(\tilde{{\mathbf{X}}},\tilde{{\mathbf{y}}})\Vert_F^2 \\ &  \ \ \ + 2\Vert\nabla_{{\mathbf{y}}}f({\mathbf{X}},{\mathbf{y}})-\nabla_{{\mathbf{y}}}f(\tilde{{\mathbf{X}}},{\mathbf{y}})\Vert_2^2 + 2\Vert\nabla_{{\mathbf{y}}}f(\tilde{{\mathbf{X}}},{\mathbf{y}})-\nabla_{{\mathbf{y}}}f(\tilde{{\mathbf{X}}},\tilde{{\mathbf{y}}})\Vert_2^2
\\ & \le 2(\beta_{X}^2+\beta_{yX}^2)\Vert{\mathbf{X}}-\tilde{{\mathbf{X}}}\Vert_F^2 + 2(\beta_{y}^2+\beta_{Xy}^2)\Vert{\mathbf{y}}-\tilde{{\mathbf{y}}}\Vert_2^2
\\ & \le 2\max\lbrace\beta_{X}^2+\beta_{yX}^2,\beta_{y}^2+\beta_{Xy}^2\rbrace\Vert({\mathbf{X}},{\mathbf{y}})-(\tilde{{\mathbf{X}}},\tilde{{\mathbf{y}}})\Vert^2.
\end{align*}

Therefore, $\beta= \sqrt{2}\max\left\lbrace\sqrt{\beta_{X}^2+\beta_{yX}^2},\sqrt{\beta_{y}^2+\beta_{Xy}^2}\right\rbrace$.


The following lemma highlights a connection between the gradient of a saddle-point of \eqref{problem1} and subgradients of an optimal solution to \eqref{nonSmoothProblem} for which the first order optimality condition holds. One of the connections we will be interested in, is that GSC for Problem \eqref{nonSmoothProblem} implies GSC (with the same parameters) for Problem \eqref{problem1}. However, to prove this specific connection we require an additional structural assumption on the objective function $g(\cdot)$. We note that this assumption holds for all applications mentioned in this paper.
\begin{assumption}\label{ass:struct}
$g({\mathbf{X}})$ is of the form $g({\mathbf{X}}) = h({\mathbf{X}}) + \max_{{\mathbf{y}}\in\mathcal{K}}{\mathbf{y}}^{\top}(\mathcal{A}({\mathbf{X}})-{\mathbf{b}})$, where $h(\cdot)$ is smooth and convex, and $\mathcal{A}$ is a linear map.
\end{assumption}

\begin{lemma} \label{lemma:connectionSubgradientNonSmoothAndSaddlePoint}
If $({\mathbf{X}}^*,{\mathbf{y}}^*)$ is a saddle-point of Problem \eqref{problem1} then ${\mathbf{X}}^*$ is an optimal solution to Problem \eqref{nonSmoothProblem},  $\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*)\in\partial g({\mathbf{X}}^*)$,  and for all ${\mathbf{X}}\in{\mathcal{S}_n}$ it holds that $\langle{\mathbf{X}}-{\mathbf{X}}^*,\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*)\rangle\ge0$.
Conversely, under \cref{ass:struct}, if ${\mathbf{X}}^*$ is an optimal solution to Problem \eqref{nonSmoothProblem}, and ${\mathbf{G}}^*\in\partial{}g({\mathbf{X}}^*)$ which satisfies $\langle{\mathbf{X}}-{\mathbf{X}}^*,{\mathbf{G}}^*\rangle\ge0$  for all ${\mathbf{X}}\in{\mathcal{S}_n}$, then there exists ${\mathbf{y}}^*\in\argmax_{{\mathbf{y}}\in\mathcal{K}}f({\mathbf{X}}^*,{\mathbf{y}})$ such that $({\mathbf{X}}^*,{\mathbf{y}}^*)$ is a saddle-point of Problem \eqref{problem1}, and $\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*)={\mathbf{G}}^*$.
\end{lemma}

The proof is given in \cref{sec:proofSubgradEquiv}. 
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 \cref{lemma:connectionSubgradientNonSmoothAndSaddlePoint}, naturally
leads to the formulation of the following generalized strict complementarity assumption for saddle-point problems.

\begin{assumption}[generalized strict complementarity for saddle-points]\label{ass:strictcompSaddlePoint}
We say a saddle-point $({\mathbf{X}}^*,{\mathbf{y}}^*)\in{\mathcal{S}_n}\times\mathcal{K}$ of Problem \eqref{problem1} with $\textrm{rank}({\mathbf{X}}^*)=r^*$ satisfies the generalized strict complementarity assumption with parameters $r\geq r^*, \delta>0$, if $\lambda_{n-r}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))-\lambda_{n}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))\ge\delta$.
\end{assumption}


\begin{remark}
Note that under \cref{ass:struct}, due to Lemma \ref{lemma:connectionSubgradientNonSmoothAndSaddlePoint},  GSC with parameters $r,\delta$ for some optimal solution ${\mathbf{X}}^*$ to Problem \eqref{nonSmoothProblem} implies GSC with parameters $r,\delta$ to a corresponding saddle-point $({\mathbf{X}}^*,{\mathbf{y}}^*)$ of Problem  \eqref{problem1}. Nevertheless, \cref{ass:struct} is not necessary for proving our convergence results for Problem \eqref{problem1}, which are directly stated in terms of Assumption \ref{ass:strictcompSaddlePoint}.
\end{remark}


\section{Projected Extragradient Method with Low-Rank Projections}

In this section we formally state and prove our main result: the projected extragradient method for the saddle-point Problem \eqref{problem1}, when initialized in the proximity of a saddle-point which satisfies GSC (\cref{ass:strictcompSaddlePoint}), converges with its original $O(1/t)$ rate while requiring only two low-rank SVD computations per iteration.

\begin{algorithm}[H]
	\caption{Projected extragradient method for saddle-point problems (see also \cite{extragradientK,Nemirovski})}\label{alg:EG}
	\begin{algorithmic}
		\STATE \textbf{Input:} sequence of step-sizes $\{\eta_t\}_{t\geq 1}$ 
		\STATE \textbf{Initialization:} $({\mathbf{X}}_1,{\mathbf{y}}_1)\in{\mathcal{S}_n}\times\mathcal{K}$
		\FOR{$t = 1,2,...$} 
            \STATE ${\mathbf{Z}}_{t+1}=\Pi_{{\mathcal{S}_n}}[{\mathbf{X}}_t-\eta_t\nabla_{{\mathbf{X}}}f({\mathbf{X}}_t,{\mathbf{y}}_t)]$   
			\STATE ${\mathbf{w}}_{t+1}=\Pi_{\mathcal{K}}[{\mathbf{y}}_t+\eta_t\nabla_{{\mathbf{y}}}f({\mathbf{X}}_t,{\mathbf{y}}_t)]$
			\STATE ${\mathbf{X}}_{t+1}=\Pi_{{\mathcal{S}_n}}[{\mathbf{X}}_t-\eta_t\nabla_{{\mathbf{X}}}f({\mathbf{Z}}_{t+1},{\mathbf{w}}_{t+1})]$   
			\STATE ${\mathbf{y}}_{t+1}=\Pi_{\mathcal{K}}[{\mathbf{y}}_t+\eta_t\nabla_{{\mathbf{y}}}f({\mathbf{Z}}_{t+1},{\mathbf{w}}_{t+1})]$
        \ENDFOR
	\end{algorithmic}
\end{algorithm}

First, in the following lemma we state the standard convergence result of the projected extragradient method, which is a well known result\footnote{\cite{sebastienMD,Nemirovski} prove this result with respect to the ergodic series. A small adjustment of the proof proves the same with respect to the minimum and maximum iterates.}. For completeness we include the proof in \cref{appendixConvergenceProof}.



\begin{lemma} \label{lemma:EGconvergenceRate}
Let $\lbrace({\mathbf{X}}_t,{\mathbf{y}}_t)\rbrace_{t\ge1}$ and $\lbrace({\mathbf{Z}}_{t},{\mathbf{w}}_t)\rbrace_{t\ge2}$ be the sequences generated by \cref{alg:EG} with a fixed step-size $\eta_t=\eta\le\min\left\lbrace\frac{1}{\beta_{X}+\beta_{Xy}},\frac{1}{\beta_{y}+\beta_{yX}},\frac{1}{\beta_{X}+\beta_{yX}},\frac{1}{\beta_{y}+\beta_{Xy}}\right\rbrace$ then
\begin{align*}
\frac{1}{T}\sum_{t=1}^T \max_{{\mathbf{y}}\in\mathcal{K}} f({\mathbf{Z}}_{t+1},{\mathbf{y}}) - \frac{1}{T}\sum_{t=1}^T \min_{{\mathbf{X}}\in{\mathcal{S}_n}} f({\mathbf{X}},{\mathbf{w}}_{t+1}) & \le \frac{D^2}{2\eta T},
\end{align*}
where $D:=\sup_{({\mathbf{X}},{\mathbf{y}}),(\tilde{{\mathbf{X}}},\tilde{{\mathbf{y}}})\in{{\mathcal{S}_n}\times\mathcal{K}}}\Vert({\mathbf{X}},{\mathbf{y}})-(\tilde{{\mathbf{X}}},\tilde{{\mathbf{y}}})\Vert$. In particular,
\begin{enumerate}
\item $\min_{t\in[T]}\left(\max_{{\mathbf{y}}\in\mathcal{K}} f({\mathbf{Z}}_{t+1},{\mathbf{y}}) - \min_{{\mathbf{X}}\in{\mathcal{S}_n}} f({\mathbf{X}},{\mathbf{w}}_{t+1})\right)
\le \frac{D^2}{2\eta T}$.
\item $\max_{{\mathbf{y}}\in\mathcal{K}} f\left(\frac{1}{T}\sum_{t=1}^T {\mathbf{Z}}_{t+1},{\mathbf{y}}\right) -  \min_{{\mathbf{X}}\in{\mathcal{S}_n}} f\left({\mathbf{X}},\frac{1}{T}\sum_{t=1}^T{\mathbf{w}}_{t+1}\right) \le \frac{D^2}{2\eta T}$.
\end{enumerate}
\end{lemma}



We can now state our main theorem.

\begin{theorem}[main theorem] \label{theroem:allPutTogether}
Fix an optimal solution $({\mathbf{X}}^*,{\mathbf{y}}^*)\in{{\mathcal{S}_n}\times\mathcal{K}}$ to Problem \eqref{problem1}. Let $\tilde{r}$ denote the multiplicity of $\lambda_{n}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))$ and for any $r\ge \tilde{r}$ define $\delta(r)=\lambda_{n-r}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*)-\lambda_{n}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*)$. Let $\lbrace({\mathbf{X}}_t,{\mathbf{y}}_t)\rbrace_{t\ge1}$ and $\lbrace({\mathbf{Z}}_{t},{\mathbf{w}}_t)\rbrace_{t\ge2}$ be the sequences of iterates generated by \cref{alg:EG} with a fixed step-size \[ \eta=\min\Bigg\lbrace\frac{1}{2\sqrt{\beta_X^2+\beta_{yX}^2}},\frac{1}{2\sqrt{\beta_y^2+\beta_{Xy}^2}},\frac{1}{\beta_{X}+\beta_{Xy}},\frac{1}{\beta_y+\beta_{yX}}\Bigg\rbrace. \] 
Assume the initialization $({\mathbf{X}}_1,{\mathbf{y}}_1)$ satisfies $\Vert({\mathbf{X}}_1,{\mathbf{y}}_1)-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert\le R_0(r)$, where
\begin{align*}
R_0(r):= \frac{\eta}{(1+\sqrt{2})\left(1+(2+\sqrt{2})\eta\max\lbrace \beta_{X},\beta_{Xy}\rbrace\right)}\max\Bigg\{\frac{\sqrt{\tilde{r}}\delta(r-\tilde{r}+1)}{2},\frac{\delta(r)}{(1+1/\sqrt{\tilde{r}})}\Bigg\}.
\end{align*}
Then, for all $t\ge1$, the projections $\Pi_{{\mathcal{S}_n}}[{\mathbf{X}}_t-\eta\nabla_{{\mathbf{X}}}f({\mathbf{X}}_t,{\mathbf{y}}_t)]$ and $\Pi_{{\mathcal{S}_n}}[{\mathbf{X}}_t-\eta\nabla_{{\mathbf{X}}}f({\mathbf{Z}}_{t+1},{\mathbf{w}}_{t+1})]$ can be replaced with their rank-r truncated counterparts (see \eqref{truncatedProjection}) without changing the sequences $\lbrace({\mathbf{X}}_t,{\mathbf{y}}_t)\rbrace_{t\ge1}$ and $\lbrace({\mathbf{Z}}_{t},{\mathbf{w}}_t)\rbrace_{t\ge2}$, and for any $T\ge0$ it holds that
\begin{align*}
& \frac{1}{T}\sum_{t=1}^T \max_{{\mathbf{y}}\in\mathcal{K}} f({\mathbf{Z}}_{t+1},{\mathbf{y}}) - \frac{1}{T}\sum_{t=1}^T \min_{{\mathbf{X}}\in{\mathcal{S}_n}} f({\mathbf{X}},{\mathbf{w}}_{t+1})
\\ & \le \frac{D^2\max{\left\lbrace\sqrt{\beta_X^2+\beta_{yX}^2},\sqrt{\beta_y^2+\beta_{Xy}^2},\frac{1}{2}(\beta_{X}+\beta_{Xy}),\frac{1}{2}(\beta_y+\beta_{yX})\right\rbrace}}{T},
\end{align*} 
where $D:=\sup_{({\mathbf{X}},{\mathbf{y}}),({\mathbf{Z}},{\mathbf{w}})\in{{\mathcal{S}_n}\times\mathcal{K}}}\Vert({\mathbf{X}},{\mathbf{y}})-({\mathbf{Z}},{\mathbf{w}})\Vert$.
\end{theorem}

\begin{remark}
Note that \cref{theroem:allPutTogether} implies that if standard strict complementarity holds for Problem \eqref{problem1}, that is \cref{ass:strictcompSaddlePoint} holds with $r=r^*=\textrm{rank}({\mathbf{X}}^*)$ and some $\delta>0$, then only rank-$r^*$ SVDs are required  so that Algorithm \ref{alg:EG} 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 \ref{alg:EG} needs to be initialized quite significantly.
\end{remark}


To prove \cref{theroem:allPutTogether} we first prove two technical lemmas. We begin by proving that the iterates of Algorithm \ref{alg:EG} always remain inside a ball of a certain radius around an optimal solution.



\begin{lemma} \label{lemma:convergenceOfXYZWsequences}
Let $\lbrace({\mathbf{X}}_t,{\mathbf{y}}_t)\rbrace_{t\ge1}$ and $\lbrace({\mathbf{Z}}_{t},{\mathbf{w}}_t)\rbrace_{t\ge2}$ be the sequences generated by \cref{alg:EG} with a step-size $\eta_t\le\frac{1}{\beta}$, and let $({\mathbf{X}}^*,{\mathbf{y}}^*)$ be some optimal solution to Problem \eqref{problem1}. Then for all $t\ge1$ it holds that 
\begin{align*}
\Vert({\mathbf{X}}_{t+1},{\mathbf{y}}_{t+1})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert & \le \Vert({\mathbf{X}}_{t},{\mathbf{y}}_t)-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert,
\\ \Vert({\mathbf{Z}}_{t+1},{\mathbf{w}}_{t+1})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert & \le \Bigg(1+\frac{1}{\sqrt{1-\eta_{t}^2\beta^2}}\Bigg)\Vert({\mathbf{X}}_t,{\mathbf{y}}_t)-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert.
\end{align*}
\end{lemma}

\begin{proof}

A known inequality of the EG algorithm (see for example Lemma 12.1.10 in \cite{nonexpansiveProof}) is
\begin{align} \label{ineq:EGconvergenceOfSeries}
& \Vert({\mathbf{X}}_{t+1},{\mathbf{y}}_{t+1})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert^2 \nonumber
\\& \le \Vert({\mathbf{X}}_{t},{\mathbf{y}}_t)-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert^2 - (1-\eta_t^2\beta^2)\Vert({\mathbf{X}}_t,{\mathbf{y}}_t)-({\mathbf{Z}}_{t+1},{\mathbf{w}}_{t+1})\Vert^2.
\end{align}

Since $\eta_t^2\beta^2\le1$ it follows that
\begin{align*}
\Vert({\mathbf{X}}_{t+1},{\mathbf{y}}_{t+1})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert \le \Vert({\mathbf{X}}_{t},{\mathbf{y}}_t)-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert.
\end{align*}

In addition, using \eqref{ineq:EGconvergenceOfSeries}
\begin{align*}
\Vert({\mathbf{X}}_t,{\mathbf{y}}_t)-({\mathbf{Z}}_{t+1},{\mathbf{w}}_{t+1})\Vert & \le \sqrt{(1-\eta_t^2\beta^2)^{-1}}\Vert({\mathbf{X}}_{t},{\mathbf{y}}_t)-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert.
\end{align*}

Therefore,
\begin{align*}
\Vert({\mathbf{Z}}_{t+1},{\mathbf{w}}_{t+1})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert & \le \Vert({\mathbf{Z}}_{t+1},{\mathbf{w}}_{t+1})-({\mathbf{X}}_{t},{\mathbf{y}}_{t})\Vert+\Vert({\mathbf{X}}_{t},{\mathbf{y}}_{t})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert
\\ & \le \sqrt{(1-\eta_{t}^2\beta^2)^{-1}}\Vert({\mathbf{X}}_{t},{\mathbf{y}}_{t})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert \\
& \ \ \ +\Vert({\mathbf{X}}_{t},{\mathbf{y}}_{t})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert
\\ & = \Bigg(1+\frac{1}{\sqrt{1-\eta_{t}^2\beta^2}}\Bigg)\Vert({\mathbf{X}}_{t},{\mathbf{y}}_{t})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert.
\end{align*}

\end{proof}




We now prove that when close enough to a low-rank saddle-point of Problem \eqref{problem1}, under an assumption of an eigen-gap in the gradient of the saddle-point, both projections onto the spectrahedron that are necessary in each iteration of \cref{alg:EG}, result in low-rank matrices.

\begin{lemma}
\label{lemma:radiusOfLowRankProjections}
Let $({\mathbf{X}}^*,{\mathbf{y}}^*)$ be an optimal solution to Problem \eqref{problem1}. Let $\tilde{r}$ denote the multiplicity of $\lambda_{n}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))$ and for any $r\ge \tilde{r}$ denote $\delta(r):=\lambda_{n-r}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))-\lambda_{n}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))$.
Then, for any $\eta\ge0$ and $({\mathbf{X}},{\mathbf{y}})\in\mathcal{S}_n\times\mathcal{K}$, if
\begin{align*}
& \Vert({\mathbf{X}},{\mathbf{y}})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert
\\ & \le \frac{\eta}{1+\sqrt{2}\eta\max\lbrace \beta_{X},\beta_{Xy}\rbrace\left(1+\frac{1}{\sqrt{1-\eta^2\beta^2}}\right)}\max\left\lbrace\frac{\sqrt{\tilde{r}}\delta(r-\tilde{r}+1)}{2},\frac{\delta(r)}{(1+1/\sqrt{\tilde{r}})}\right\rbrace
\end{align*}
then $\textrm{rank}\left( \Pi_{\mathcal{S}_n}[{\mathbf{X}}-\eta \nabla_{{\mathbf{X}}}f({\mathbf{X}},{\mathbf{y}})]\right)\le r$ and $\textrm{rank}\left( \Pi_{\mathcal{S}_n}[{\mathbf{X}}-\eta \nabla_{{\mathbf{X}}}f({\mathbf{Z}}_+,{\mathbf{w}}_+)]\right)\le r$ where ${\mathbf{Z}}_+=\Pi_{{\mathcal{S}_n}}[{\mathbf{X}}-\eta\nabla_{{\mathbf{X}}}f({\mathbf{X}},{\mathbf{y}})]$ and ${\mathbf{w}}_+=\Pi_{\mathcal{K}}[{\mathbf{y}}-\eta\nabla_{{\mathbf{y}}}f({\mathbf{X}},{\mathbf{y}})]$. 
\end{lemma}

\begin{proof}
Denote ${\mathbf{P}}^*={\mathbf{X}}^*-\eta\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*)$. By \cref{lemma:connectionSubgradientNonSmoothAndSaddlePoint}, $\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*)$ is a subgradient of the corresponding nonsmooth objective $g({\mathbf{X}})=\max_{{\mathbf{y}}\in\mathcal{K}}f({\mathbf{X}},{\mathbf{y}})$ at the point ${\mathbf{X}}^*$. Moreover, this subgradient also satisfies the first-order optimality condition. Hence, invoking \cref{lemma:Xopt-subgradEigen} with this subgradient  we have that
\begin{align} \label{eq:eigsOfPstar}
 \forall i\le \textrm{rank}({\mathbf{X}}^*) & :\ \lambda_i({\mathbf{P}}^*) = \lambda_i({\mathbf{X}}^*)-\eta\lambda_n(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*)); \nonumber \\
\forall i>\textrm{rank}({\mathbf{X}}^*) & :\ \lambda_i({\mathbf{P}}^*) = -\eta\lambda_{n-i+1}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*)). 
\end{align}

Therefore, using \eqref{eq:eigsOfPstar} and the fact that $\lambda_{n-i+1}(\nabla{}f_{{\mathbf{X}}}({\mathbf{X}}^*,{\mathbf{y}}^*)) = \lambda_{n}(\nabla{}f_{{\mathbf{X}}}({\mathbf{X}}^*,{\mathbf{y}}^*))$ for all $i\leq \tilde{r}$ we have,
\begin{align} \label{ineq:sumOfXeigsInProof}
\sum_{i=1}^{\tilde{r}} {\lambda_i({\mathbf{P}}^*)} & = \sum_{i=1}^{\tilde{r}} {\lambda_i({\mathbf{X}}^*-\eta\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))}
= \sum_{i=1}^{\tilde{r}} \lambda_i({\mathbf{X}}^*) -\eta\sum_{i=1}^{\tilde{r}} \lambda_{n-i+1}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*)) \nonumber
\\ & = \sum_{i=1}^{\textrm{rank}({\mathbf{X}}^*)}\lambda_i({\mathbf{X}}^*) -\eta \sum_{i=1}^{\tilde{r}}\lambda_n(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*)) = 1 -\eta \tilde{r}\lambda_n(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*)).
\end{align}

Let ${\mathbf{P}}\in\mathbb{S}^n$.  From the structure of the Euclidean projection onto the spectrahedron (see Eq. \eqref{eq:euclidProj}), it follows that
a sufficient condition so that $\textrm{rank}\left( \Pi_{\mathcal{S}_n}[{\mathbf{P}}]\right)\le r$ is that $\sum_{i=1}^{r} {\lambda_i({\mathbf{P}})}-r\lambda_{r+1}({\mathbf{P}}) \ge 1$. We will bound the LHS of this inequality.

First, it holds that
\begin{align} \label{ineq:sumTildeRbound}
\sum_{i=1}^{\tilde{r}} {\lambda_i({\mathbf{P}})} & \underset{(a)}{\ge} \sum_{i=1}^{\tilde{r}} {\lambda_i({\mathbf{P}}^*)} - \sum_{i=1}^{\tilde{r}} {\lambda_i({\mathbf{P}}^*-{\mathbf{P}})} 
\ge \sum_{i=1}^{\tilde{r}} {\lambda_i({\mathbf{P}}^*)} - \sqrt{\tilde{r}\sum_{i=1}^{\tilde{r}} {\lambda_i^2({\mathbf{P}}-{\mathbf{P}}^*)}} \nonumber
\\ & \ge \sum_{i=1}^{\tilde{r}} {\lambda_i({\mathbf{P}}^*)} - \sqrt{\tilde{r}\sum_{i=1}^{n} {\lambda_i^2({\mathbf{P}}-{\mathbf{P}}^*)}} 
\ge \sum_{i=1}^{\tilde{r}} {\lambda_i({\mathbf{P}}^*)} - \sqrt{\tilde{r}}\Vert{\mathbf{P}}-{\mathbf{P}}^*\Vert_F \nonumber
\\ & \underset{(b)}{\ge}  1 -\eta \tilde{r}\lambda_n(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))- \sqrt{\tilde{r}}\Vert{\mathbf{P}}-{\mathbf{P}}^*\Vert_F,
\end{align}
where (a) holds from Ky Fan’s inequality for eigenvalues and (b) holds from \eqref{ineq:sumOfXeigsInProof}.

Now, for any $r\ge \tilde{r}$ using Weyl's inequality and \eqref{eq:eigsOfPstar}
\begin{align} \label{ineq:Rplus1bound1}
\lambda_{r+1}({\mathbf{P}}) & \le \lambda_{r+1}({\mathbf{P}}^*)+\lambda_{1}({\mathbf{P}}-{\mathbf{P}}^*) \le \lambda_{r+1}({\mathbf{P}}^*)+\Vert{\mathbf{P}}-{\mathbf{P}}^*\Vert_F \nonumber
\\ & = -\eta\lambda_{n-r}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))+\Vert{\mathbf{P}}-{\mathbf{P}}^*\Vert_F.
\end{align}

Thus, combining \eqref{ineq:sumTildeRbound} and \eqref{ineq:Rplus1bound1} we obtain
\begin{align} \label{ineq:opt1_Pbound}
\sum_{i=1}^{r} {\lambda_i({\mathbf{P}})}-r\lambda_{r+1}({\mathbf{P}})  \ge \sum_{i=1}^{\tilde{r}} {\lambda_i({\mathbf{P}})}-\tilde{r}\lambda_{r+1}({\mathbf{P}})
\ge 1 +\eta \tilde{r}\delta(r) - (\tilde{r}+\sqrt{\tilde{r}})\Vert{\mathbf{P}}-{\mathbf{P}}^*\Vert_F.
\end{align} 

Alternatively, if $r\ge2\tilde{r}-1$ then using the general Weyl inequality and \eqref{eq:eigsOfPstar} we obtain
\begin{align} \label{ineq:Rplus1bound2}
\lambda_{r+1}({\mathbf{P}}) & \le \lambda_{r-\tilde{r}+2}({\mathbf{P}}^*)+\lambda_{\tilde{r}}({\mathbf{P}}-{\mathbf{P}}^*) = \lambda_{r-\tilde{r}+2}({\mathbf{P}}^*)+\sqrt{\lambda^2_{\tilde{r}}({\mathbf{P}}-{\mathbf{P}}^*)} \nonumber
\\ & \le \lambda_{r-\tilde{r}+2}({\mathbf{P}}^*)+\frac{1}{\sqrt{\tilde{r}}}\Vert{\mathbf{P}}-{\mathbf{P}}^*\Vert_F
= -\eta\lambda_{n-r+\tilde{r}-1}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))+\frac{1}{\sqrt{\tilde{r}}}\Vert{\mathbf{P}}-{\mathbf{P}}^*\Vert_F.
\end{align} 

Thus, combining \eqref{ineq:sumTildeRbound} and \eqref{ineq:Rplus1bound2} we obtain
\begin{align} \label{ineq:opt2_Pbound}
& \sum_{i=1}^{r} {\lambda_i({\mathbf{P}})}-r\lambda_{r+1}({\mathbf{P}})  \ge \sum_{i=1}^{\tilde{r}} {\lambda_i({\mathbf{P}})}-\tilde{r}\lambda_{r+1}({\mathbf{P}})
\ge 1 +\eta \tilde{r}\delta(r-\tilde{r}+1) - 2\sqrt{\tilde{r}}\Vert{\mathbf{P}}-{\mathbf{P}}^*\Vert_F.
\end{align} 

Now we are left with bounding $\Vert{\mathbf{P}}-{\mathbf{P}}^*\Vert_F$. Note that by the smoothness of $f$, for any $({\mathbf{X}},{\mathbf{y}})\in{\mathcal{S}_n}\times\mathcal{K}$ it holds that
\begin{align} \label{ineq:betaSmoothnessInProof}
& \Vert\nabla_{{\mathbf{X}}}f({\mathbf{X}},{\mathbf{y}})-\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert_F  \nonumber
\\ & \le \Vert\nabla_{{\mathbf{X}}}f({\mathbf{X}},{\mathbf{y}})-\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}})\Vert_F+\Vert\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}})-\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert_F \nonumber
\\ & \le \beta_{X}\Vert{\mathbf{X}}-{\mathbf{X}}^*\Vert_F+\beta_{Xy}\Vert{\mathbf{y}}-{\mathbf{y}}^*\Vert_2.
\end{align}

Taking ${\mathbf{P}}={\mathbf{X}}-\eta\nabla_{{\mathbf{X}}}f({\mathbf{X}},{\mathbf{y}})$ we get
\begin{align*}
\Vert{\mathbf{P}}-{\mathbf{P}}^*\Vert_F & = \Vert{\mathbf{X}}-\eta\nabla_{{\mathbf{X}}}f({\mathbf{X}},{\mathbf{y}})-{\mathbf{X}}^*+\eta\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert_F
\\ & \le \Vert{\mathbf{X}}-{\mathbf{X}}^*\Vert_F+\eta\Vert\nabla_{{\mathbf{X}}}f({\mathbf{X}},{\mathbf{y}})-\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert_F
\\ & \le \Vert({\mathbf{X}},{\mathbf{y}})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert+\eta\Vert\nabla_{{\mathbf{X}}}f({\mathbf{X}},{\mathbf{y}})-\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert_F
\\ & \le \Vert({\mathbf{X}},{\mathbf{y}})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert+\eta\beta_{X}\Vert{\mathbf{X}}-{\mathbf{X}}^*\Vert_F+\eta\beta_{Xy}\Vert{\mathbf{y}}-{\mathbf{y}}^*\Vert_2,
\end{align*}
where the last inequality holds from \eqref{ineq:betaSmoothnessInProof}.

For any $a,b\ge0$ it holds that 
\begin{align*}
a\Vert{\mathbf{X}}-{\mathbf{X}}^*\Vert_F+b\Vert{\mathbf{y}}-{\mathbf{y}}^*\Vert_2 & \le \max\lbrace a,b\rbrace\left(\Vert{\mathbf{X}}-{\mathbf{X}}^*\Vert_F+\Vert{\mathbf{y}}-{\mathbf{y}}^*\Vert_2\right)
\\ & \le \sqrt{2}\max\lbrace a,b\rbrace \Vert({\mathbf{X}},{\mathbf{y}})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert.
\end{align*}

Thus, by taking $a=\eta\beta_{X}$ and $b=\eta\beta_{Xy}$ we obtain
\begin{align} \label{ineq:p1NormBound}
\Vert{\mathbf{P}}-{\mathbf{P}}^*\Vert_F & \le \left(1+\sqrt{2}\eta\max\lbrace\beta_{X},\beta_{Xy}\rbrace\right)\Vert({\mathbf{X}},{\mathbf{y}})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert.
\end{align}

Therefore, plugging \eqref{ineq:p1NormBound} into \eqref{ineq:opt1_Pbound} we obtain that the condition $\sum_{i=1}^{r} {\lambda_i({\mathbf{P}})}-r\lambda_{r+1}({\mathbf{P}})\ge1$ holds if
\begin{align*} 
\Vert({\mathbf{X}},{\mathbf{y}})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert \le \frac{\eta\delta(r)}{(1+1/\sqrt{\tilde{r}})\left(1+\sqrt{2}\eta\max\lbrace\beta_{X},\beta_{Xy}\rbrace\right)}.
\end{align*}

Alternatively, plugging \eqref{ineq:p1NormBound} into \eqref{ineq:opt2_Pbound} we obtain that if $r\ge2\tilde{r}-1$ then the condition $\sum_{i=1}^{r} {\lambda_i({\mathbf{P}})}-r\lambda_{r+1}({\mathbf{P}})\ge1$ holds if
\begin{align*}
\Vert({\mathbf{X}},{\mathbf{y}})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert \le \frac{\eta\sqrt{\tilde{r}}\delta(r-\tilde{r}+1)}{2\left(1+\sqrt{2}\eta\max\lbrace\beta_{X},\beta_{Xy}\rbrace\right)}.
\end{align*}

Note that $\delta(r-\tilde{r}+1)>0$ only if $r\ge2\tilde{r}-1$. Therefore, we can combine the last two inequalities to conclude that for any $r\ge\tilde{r}$ if
\begin{align} \label{ineq:XseriesRadius}
\Vert({\mathbf{X}},{\mathbf{y}})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert \le \frac{\eta}{1+\sqrt{2}\eta\max\lbrace\beta_{X},\beta_{Xy}\rbrace}\max\Bigg\{\frac{\sqrt{\tilde{r}}\delta(r-\tilde{r}+1)}{2},\frac{\delta(r)}{(1+1/\sqrt{\tilde{r}})}\Bigg\}
\end{align}
then $\textrm{rank}(\Pi_{{\mathcal{S}_n}}[{\mathbf{X}}-\nabla_{{\mathbf{X}}}f({\mathbf{X}},{\mathbf{y}})])\le r$.

Similarly, taking ${\mathbf{P}}={\mathbf{X}}-\eta\nabla_{{\mathbf{X}}}f({\mathbf{Z}}_+,{\mathbf{w}}_+)$ we get
\begin{align*} 
\Vert{\mathbf{P}}-{\mathbf{P}}^*\Vert_F & = \Vert{\mathbf{X}}-\eta\nabla_{{\mathbf{X}}}f({\mathbf{Z}}_+,{\mathbf{w}}_+)-{\mathbf{X}}^*+\eta\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert_F 
\\ & \le \Vert{\mathbf{X}}-{\mathbf{X}}^*\Vert_F+\eta\Vert\nabla_{{\mathbf{X}}}f({\mathbf{Z}}_+,{\mathbf{w}}_+)-\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert_F
\\ & \le \Vert{\mathbf{X}}-{\mathbf{X}}^*\Vert_F+\eta\beta_{X}\Vert{\mathbf{Z}}_+-{\mathbf{X}}^*\Vert_F+\eta\beta_{Xy}\Vert{\mathbf{w}}_+-{\mathbf{y}}^*\Vert_2,
\end{align*}
where the last inequality holds from \eqref{ineq:betaSmoothnessInProof}.

For any $a,b,c\ge0$ it holds that 
\begin{align*}
& a\Vert{\mathbf{X}}-{\mathbf{X}}^*\Vert_F+b\Vert{\mathbf{Z}}_+-{\mathbf{X}}^*\Vert_F+c\Vert{\mathbf{w}}_+-{\mathbf{y}}^*\Vert_2 
\\ & \le a\Vert{\mathbf{X}}-{\mathbf{X}}^*\Vert_F+\max\lbrace b,c\rbrace\left(\Vert{\mathbf{Z}}_+-{\mathbf{X}}^*\Vert_F+\Vert{\mathbf{w}}_+-{\mathbf{y}}^*\Vert_2\right)
\\ & \le a\Vert{\mathbf{X}}-{\mathbf{X}}^*\Vert_F+\sqrt{2}\max\lbrace b,c\rbrace \Vert({\mathbf{Z}}_+,{\mathbf{w}}_+)-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert
\\ & \le a\Vert({\mathbf{X}},{\mathbf{y}})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert+\sqrt{2}\max\lbrace b,c\rbrace \Vert({\mathbf{Z}}_+,{\mathbf{w}}_+)-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert
\\ & \le \left(a+\sqrt{2}\max\lbrace b,c\rbrace\left(1+\frac{1}{\sqrt{1-\eta^2\beta^2}}\right)\right)\Vert({\mathbf{X}},{\mathbf{y}})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert,
\end{align*}
where the second to last inequality holds from \cref{lemma:convergenceOfXYZWsequences}.

Thus, by taking $a=1$, $b=\eta\beta_{X}$, and $c=\eta\beta_{Xy}$ we obtain
\begin{align} \label{ineq:p2NormBound}
\Vert{\mathbf{P}}-{\mathbf{P}}^*\Vert_F & \le \left(1+\sqrt{2}\eta\max\lbrace \beta_{X},\beta_{Xy}\rbrace\left(1+\frac{1}{\sqrt{1-\eta^2\beta^2}}\right)\right)\Vert({\mathbf{X}},{\mathbf{y}})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert.
\end{align}

Therefore, plugging \eqref{ineq:p2NormBound} into \eqref{ineq:opt1_Pbound} we obtain that the condition $\sum_{i=1}^{r} {\lambda_i({\mathbf{P}})}-r\lambda_{r+1}({\mathbf{P}})\ge1$ holds if
\begin{align*}
\Vert({\mathbf{X}},{\mathbf{y}})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert \le \frac{\eta\delta(r)}{(1+1/\sqrt{r^*})\left(1+\sqrt{2}\eta\max\lbrace \beta_{X},\beta_{Xy}\rbrace\left(1+\frac{1}{\sqrt{1-\eta^2\beta^2}}\right)\right)}.
\end{align*}

Alternatively, plugging \eqref{ineq:p2NormBound} into \eqref{ineq:opt2_Pbound} we obtain that if $r\ge2\tilde{r}-1$ then the condition $\sum_{i=1}^{r} {\lambda_i({\mathbf{P}})}-r\lambda_{r+1}({\mathbf{P}})\ge1$ holds if
\begin{align*}
\Vert({\mathbf{X}},{\mathbf{y}})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert \le \frac{\eta\sqrt{\tilde{r}}\delta(r-\tilde{r}+1)}{1+\sqrt{2}\eta\max\lbrace \beta_{X},\beta_{Xy}\rbrace\left(1+\frac{1}{\sqrt{1-\eta^2\beta^2}}\right)}.
\end{align*}

Note that $\delta(r-\tilde{r}+1)>0$ only if $r\ge2\tilde{r}-1$. Therefore, we can combine the last two inequalities to conclude that for any $r\ge\tilde{r}$ if
\begin{align} \label{ineq:ZseriesRadius}
& \Vert({\mathbf{X}},{\mathbf{y}})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert \nonumber
\\ & \le \frac{\eta}{1+\sqrt{2}\eta\max\lbrace \beta_{X},\beta_{Xy}\rbrace\left(1+\frac{1}{\sqrt{1-\eta^2\beta^2}}\right)}\max\left\lbrace\frac{\sqrt{\tilde{r}}\delta(r-\tilde{r}+1)}{2},\frac{\delta(r)}{(1+1/\sqrt{\tilde{r}})}\right\rbrace
\end{align}
then $\textrm{rank}(\Pi_{{\mathcal{S}_n}}[{\mathbf{X}}-\nabla_{{\mathbf{X}}}f({\mathbf{Z}}_+,{\mathbf{w}}_+)])\le r$.

Taking the minimum between \eqref{ineq:XseriesRadius} and \eqref{ineq:ZseriesRadius} gives us the bound on the radius in the lemma.

\end{proof}




Now we can prove \cref{theroem:allPutTogether}.

\begin{proof}[Proof of \cref{theroem:allPutTogether}]
We will prove by induction that for all $t\ge 1$ it holds that $\Vert({\mathbf{X}}_t,{\mathbf{y}}_t)-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert \le R_0(r)$ and $\Vert({\mathbf{Z}}_{t},{\mathbf{w}}_{t})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert \le \left(1+\sqrt{2}\right)R_0(r)$, thus implying through  \cref{lemma:radiusOfLowRankProjections} that all projections $\Pi_{{\mathcal{S}_n}}[{\mathbf{X}}_t-\eta\nabla_{{\mathbf{X}}}f({\mathbf{X}}_t,{\mathbf{y}}_t)]$ and $\Pi_{{\mathcal{S}_n}}[{\mathbf{X}}_t-\eta\nabla_{{\mathbf{X}}}f({\mathbf{Z}}_{t+1},{\mathbf{w}}_{t+1})]$ can be replaced with their rank-r truncated counterparts given in \eqref{truncatedProjection}, without any change to the result. 

The initialization $\Vert({\mathbf{X}}_1,{\mathbf{y}}_1)-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert\le R_0(r)$ holds trivially.
Now, by \cref{lemma:convergenceOfXYZWsequences}, using recursion, we have that for all $t\ge1$,
\begin{align*}
\Vert({\mathbf{X}}_{t+1},{\mathbf{y}}_{t+1})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert \le \Vert({\mathbf{X}}_{t},{\mathbf{y}}_t)-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert \le \cdots \le \Vert({\mathbf{X}}_{1},{\mathbf{y}}_1)-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert\le R_0(r),
\end{align*}
and for $\beta= \sqrt{2}\max\left\lbrace\sqrt{\beta_{X}^2+\beta_{yX}^2},\sqrt{\beta_{y}^2+\beta_{Xy}^2}\right\rbrace$ we have that,
\begin{align*}
\Vert({\mathbf{Z}}_{t+1},{\mathbf{w}}_{t+1})-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert & \le \Bigg(1+\frac{1}{\sqrt{1-\eta_{t}^2\beta^2}}\Bigg)\Vert({\mathbf{X}}_t,{\mathbf{y}}_t)-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert
\\ & \le \Bigg(1+\frac{1}{\sqrt{1-\eta_{t}^2\beta^2}}\Bigg)\Vert({\mathbf{X}}_1,{\mathbf{y}}_1)-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert
\\ & \le (1+\sqrt{2})\Vert({\mathbf{X}}_1,{\mathbf{y}}_1)-({\mathbf{X}}^*,{\mathbf{y}}^*)\Vert\le (1+\sqrt{2})R_0(r).
\end{align*}

Therefore, under the assumptions of the theorem, \cref{alg:EG} can be run using only rank-r truncated projections, while maintaining its original convergence rate stated in \cref{lemma:EGconvergenceRate}.
\end{proof}



\begin{remark}
A downside of considering the saddle-point formulation \eqref{problem1} when attempting to solve Problem \eqref{nonSmoothProblem} that arises from Theorem \ref{theroem:allPutTogether}, is that not only do we need a ``warm-start'' initialization for the original primal matrix variable ${\mathbf{X}}$, in the saddle-point formulation we need a ``warm-start''  for the primal-dual pair $({\mathbf{X}},{\mathbf{y}})$. Nevertheless, as we demonstrate extensively in Section \ref{sec:expr}, it seems that very simple initialization schemes work very well in practice.
\end{remark}


\subsection{Back to nonsmooth problems}

\begin{corollary} \label{cor:nonsmoothExtragradient}
Fix an optimal solution ${\mathbf{X}}^*\in{\mathcal{S}_n}$ to Problem \eqref{nonSmoothProblem} and assume \cref{ass:struct} holds. Let ${\mathbf{G}}^*\in\partial g({\mathbf{X}}^*)$ which satisfies that $\langle{\mathbf{X}}-{\mathbf{X}}^*,{\mathbf{G}}^*\rangle\ge0$  for all ${\mathbf{X}}\in{\mathcal{S}_n}$. Let $\tilde{r}$ denote the multiplicity of $\lambda_{n}({\mathbf{G}}^*)$ and for any $r\ge \tilde{r}$ define $\delta(r):=\lambda_{n-r}({\mathbf{G}}^*)-\lambda_{n}({\mathbf{G}}^*)$.
Define $f$ as in Problem \eqref{problem1} and let $\lbrace({\mathbf{X}}_t,{\mathbf{y}}_t)\rbrace_{t\ge1}$ and $\lbrace({\mathbf{Z}}_{t},{\mathbf{w}}_t)\rbrace_{t\ge2}$ be the sequences of iterates generated by \cref{alg:EG} where $\eta$ and $R_0(r)$ are as defined in \cref{theroem:allPutTogether}.
Then, for all $t\ge1$ the projections $\Pi_{{\mathcal{S}_n}}[{\mathbf{X}}_t-\eta\nabla_{{\mathbf{X}}}f({\mathbf{X}}_t,{\mathbf{y}}_t)]$ and $\Pi_{{\mathcal{S}_n}}[{\mathbf{X}}_t-\eta\nabla_{{\mathbf{X}}}f({\mathbf{Z}}_{t+1},{\mathbf{w}}_{t+1})]$ can be replaced with rank-r truncated projections \eqref{truncatedProjection} without changing the sequences $\lbrace({\mathbf{X}}_t,{\mathbf{y}}_t)\rbrace_{t\ge1}$ and $\lbrace({\mathbf{Z}}_{t},{\mathbf{w}}_t)\rbrace_{t\ge2}$, and for any $T\ge0$ it holds that
\begin{align*}
\min_{t\in[T]}g({\mathbf{Z}}_{t+1})-g({\mathbf{X}}^*) 
& \le \frac{D^2\max{\left\lbrace\sqrt{\beta_X^2+\beta_{yX}^2},\sqrt{\beta_y^2+\beta_{Xy}^2},\frac{1}{2}(\beta_{X}+\beta_{Xy}),\frac{1}{2}(\beta_y+\beta_{yX})\right\rbrace}}{T},
\end{align*} 
where $D:=\sup_{({\mathbf{X}},{\mathbf{y}}),({\mathbf{Z}},{\mathbf{w}})\in{{\mathcal{S}_n}\times\mathcal{K}}}\Vert({\mathbf{X}},{\mathbf{y}})-({\mathbf{Z}},{\mathbf{w}})\Vert$.
\end{corollary}

\begin{proof}
Since \cref{ass:struct} holds, invoking \cref{lemma:connectionSubgradientNonSmoothAndSaddlePoint} we obtain that there exists a point ${\mathbf{y}}^*\in\argmax_{{\mathbf{y}}\in\mathcal{K}}f({\mathbf{X}}^*,{\mathbf{y}})$ such that $({\mathbf{X}}^*,{\mathbf{y}}^*)$ is a saddle-point of Problem \eqref{problem1}, and $\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*)={\mathbf{G}}^*$.
Therefore, the assumptions of \cref{theroem:allPutTogether} hold and so by \cref{theroem:allPutTogether} we get that
\begin{align} \label{ineq:boundFromTheorem1}
& \frac{1}{T}\sum_{t=1}^T \max_{{\mathbf{y}}\in\mathcal{K}} f({\mathbf{Z}}_{t+1},{\mathbf{y}}) - \frac{1}{T}\sum_{t=1}^T \min_{{\mathbf{X}}\in{\mathcal{S}_n}} f({\mathbf{X}},{\mathbf{w}}_{t+1}) \nonumber
\\ & \le \frac{D^2\max{\left\lbrace\sqrt{\beta_X^2+\beta_{yX}^2},\sqrt{\beta_y^2+\beta_{Xy}^2},\frac{1}{2}(\beta_{X}+\beta_{Xy}),\frac{1}{2}(\beta_y+\beta_{yX})\right\rbrace}}{T}.
\end{align} 
From the definition of $g$ it holds that
\begin{align} \label{ineq:replaceMaxTerm}
\min_{t\in[T]}g({\mathbf{Z}}_{t+1})=\min_{t\in[T]}\max_{{\mathbf{y}}\in\mathcal{K}}f({\mathbf{Z}}_{t+1},{\mathbf{y}}) \le \frac{1}{T}\sum_{t=1}^T\max_{{\mathbf{y}}\in\mathcal{K}}f({\mathbf{Z}}_{t+1},{\mathbf{y}})
\end{align}
and
\begin{align} \label{ineq:replaceMinTerm}
\frac{1}{T}\sum_{t=1}^T \min_{{\mathbf{X}}\in{\mathcal{S}_n}}f({\mathbf{X}},{\mathbf{w}}_{t+1}) \le \frac{1}{T}\sum_{t=1}^T f({\mathbf{X}}^*,{\mathbf{w}}_{t+1}) \le \max_{{\mathbf{y}}\in\mathcal{K}}f\left({\mathbf{X}}^*,{\mathbf{y}}\right)=g({\mathbf{X}}^*).
\end{align}
Plugging \eqref{ineq:replaceMaxTerm} and \eqref{ineq:replaceMinTerm} into the RHS of \eqref{ineq:boundFromTheorem1} we obtain the required result. 



\end{proof}



\subsection{Efficiently-computable certificates for correctness of  low-rank projections}\label{sec:certificate}

Since Theorem \ref{theroem:allPutTogether} 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 \ref{theroem:allPutTogether} hold. In practice, as long as the truncated projection (see \eqref{truncatedProjection}) equals the exact projection (see \eqref{eq:euclidProj}), we are guaranteed that Algorithm \ref{alg:EG} converges correctly with rate $O(1/t)$, without needing to verify any other condition.
Luckily, the expression in \eqref{eq:euclidProj} which characterizes the structure of the Euclidean projection onto the spectrahedron, yields exactly such a verification procedure. As already noted in  \cite{garberLowRankSmooth}, for any ${\mathbf{X}}\in\mathbb{S}^n$, we have $\widehat{\Pi}_{\mathcal{S}_n}^r[{\mathbf{X}}]=\Pi_{{\mathbf{S}}_n}[{\mathbf{X}}]$ if an only if the condition
\begin{align*}
\sum_{i=1}^r\lambda_i({\mathbf{X}})\ge 1+r\cdot\lambda_{r+1}({\mathbf{X}})
\end{align*}
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.


\section{Empirical Evidence}\label{sec:expr}


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 $\ell_1$-norm or $\ell_2$-norm regularizer / penalty term, over a $\tau$-scaled spectrahedron. In all cases the nonsmooth objective can be written as a saddle-point with function $f({\mathbf{X}},{\mathbf{y}})$ which is linear in ${\mathbf{y}}$ and in particular satisfies Assumption \ref{ass:struct}. 

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.

For all tasks considered we generate random instances, and examine the sequences of iterates generated by Algorithm \ref{alg:EG} $\lbrace({\mathbf{X}}_t,{\mathbf{y}}_t)\rbrace_{t\ge1}$, $\lbrace({\mathbf{Z}}_{t},{\mathbf{w}}_t)\rbrace_{t\ge2}$, 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 \cref{appendixDualGapCalculation}.

\begin{comment}
For a given candidate  $(\widehat{{\mathbf{Z}}},\widehat{{\mathbf{w}}})$ for an approximated saddle-point, the dual-gap is the RHS of the following Eq. \eqref{ineq:dualGapBound}:
\begin{align} \label{ineq:dualGapBound}
& g(\widehat{{\mathbf{Z}}})-g^* \leq \max_{{\mathbf{y}}\in\mathcal{K}}f(\widehat{{\mathbf{Z}}},{\mathbf{y}})-\min_{\substack{\textrm{Tr}({\mathbf{X}})=\tau,\\ {\mathbf{X}}\succeq 0}}f({\mathbf{X}},\widehat{{\mathbf{w}}}) \nonumber
\\ & \le \max_{\substack{\textrm{Tr}({\mathbf{X}})=\tau,\\ {\mathbf{X}}\succeq 0}}\langle\widehat{{\mathbf{Z}}}-{\mathbf{X}},\nabla_{{\mathbf{X}}}f(\widehat{{\mathbf{Z}}},\widehat{{\mathbf{w}}})\rangle - \min_{{\mathbf{y}}\in\mathcal{K}} \langle\widehat{{\mathbf{w}}}-{\mathbf{y}},\nabla_{{\mathbf{y}}}f(\widehat{{\mathbf{Z}}},\widehat{{\mathbf{w}}})\rangle.
\end{align} 
We prove  the correctness of the last inequality in \cref{appendixDualGapCalculation}
\end{comment}


In all tasks considered the goal is to recover a ground-truth low-rank matrix ${\mathbf{M}}_0\in\mathbb{S}^n$ from some noisy observation of it ${\mathbf{M}}={\mathbf{M}}_0+{\mathbf{N}}$, where ${\mathbf{N}}\in\mathbb{S}^{n}$ is a noise matrix. We measure the signal-to-noise ratio (SNR) as $\Vert{\mathbf{M}}_0\Vert_F^2\big/\Vert{\mathbf{N}}\Vert_F^2$. 
In all experiments we measure the relative initialization error by $\left\Vert\frac{\textrm{Tr}({\mathbf{M}}_0)}{\tau}{\mathbf{X}}_1-{\mathbf{M}}_0\right\Vert_F^2\Big{/}\left\Vert{\mathbf{M}}_0\right\Vert_F^2$, and similarly we measure the relative recovery error by $\left\Vert\frac{\textrm{Tr}({\mathbf{M}}_0)}{\tau}{\mathbf{X}}^*-{\mathbf{M}}_0\right\Vert_F^2\Big{/}\left\Vert{\mathbf{M}}_0\right\Vert_F^2$. Note that in some of the experiments we  take $\tau<\textrm{Tr}({\mathbf{M}}_0)$ to prevent the method from overfitting the noise. In addition, we measure the (standard) strict complementarity parameter which corresponds to the eigen-gap $\textrm{gap}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*)):=\lambda_{n-r}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))-\lambda_{n}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))$, $r = \textrm{rank}({\mathbf{M}}_0)$.

In all experiments we use SVDs of rank $r=\textrm{rank}({\mathbf{M}}_0)$ to compute the projections in Algorithm \ref{alg:EG} according to the truncated projection given in \eqref{truncatedProjection}. To certify the correctness of these low-rank projections (that is, that they equal the exact Euclidean projection) we confirm that the inequality
\begin{align*} 
\sum_{i=1}^r\lambda_i({\mathbf{P}}_j)\ge \tau+r\cdot\lambda_{r+1}({\mathbf{P}}_j)
\end{align*}
always holds for ${\mathbf{P}}_1={\mathbf{X}}_t-\eta\nabla_{{\mathbf{X}}}f({\mathbf{X}}_t,{\mathbf{Y}}_t)$ and ${\mathbf{P}}_2={\mathbf{X}}_t-\eta\nabla_{{\mathbf{X}}}f({\mathbf{Z}}_{t+1},{\mathbf{W}}_{t+1})$ (see also Section \ref{sec:certificate}). Indeed, we can now already state our main observation from the experiments:

\framebox{\parbox{\dimexpr\linewidth-2\fboxsep-2\fboxrule}{\itshape%
In all tasks considered and for all random instances generated, throughout all iterations of Algorithm \ref{alg:EG}, when initialized with a simple ``warm-start'' strategy and when computing only rank-$r$ truncated projections, $r=\textrm{rank}({\mathbf{M}}_0)$, the truncated projections of ${\mathbf{P}}_1={\mathbf{X}}_t-\eta\nabla_{{\mathbf{X}}}f({\mathbf{X}}_t,{\mathbf{Y}}_t)$ and ${\mathbf{P}}_2={\mathbf{X}}_t-\eta\nabla_{{\mathbf{X}}}f({\mathbf{Z}}_{t+1},{\mathbf{W}}_{t+1})$ equal their exact full-rank counterparts. That is, Algorithm \ref{alg:EG}, 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 $({\mathbf{X}}^*,{\mathbf{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 $\lambda_{n-r}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))-\lambda_{n}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))$), 2. the recovery error of the returned solution indeed improves significantly over the error of the initialization point.

 
\subsection{Sparse PCA}

We consider the sparse PCA problem in a well known convex formulation taken from \cite{d2007direct} and its equivalent saddle-point formulation:
\begin{align*}
\min_{\substack{\textrm{Tr}({\mathbf{X}}) =1,\\ {\mathbf{X}}\succeq 0}} \langle{{\mathbf{X}},-{\mathbf{M}}}\rangle + \lambda\Vert{{\mathbf{X}}}\Vert_1 =
\min_{\substack{\textrm{Tr}({\mathbf{X}}) =1,\\ {\mathbf{X}}\succeq 0}} \max_{\Vert{\mathbf{Y}}\Vert_{\infty}\le1}\lbrace\langle{\mathbf{X}},-{\mathbf{M}}\rangle+\lambda\langle{\mathbf{X}},{\mathbf{Y}}\rangle\rbrace,
\end{align*}
where ${\mathbf{M}} = {\mathbf{z}}\z^{\top} + \frac{c}{2}({\mathbf{N}}+{\mathbf{N}}^{\top})$ is a noisy observation of a rank-one matrix ${\mathbf{z}}\z^{\top}$, with ${\mathbf{z}}$ being a sparse unit vector. Each entry ${\mathbf{z}}_i$ is chosen to be $0$ with probability $0.9$ and $U\lbrace1,\ldots,10\rbrace$ with probability $0.1$, and then we normalize ${\mathbf{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=\frac{2}{\textrm{SNR}\cdot\Vert{\mathbf{N}}+{\mathbf{N}}^{\top}\Vert_F}$ for our choice of $\textrm{SNR}$.

We initialize the ${\mathbf{X}}$ variable with the rank-one approximation of ${\mathbf{M}}$. That is, we take ${\mathbf{X}}_1={\mathbf{u}}_1{\mathbf{u}}_1^{\top}$, where ${\mathbf{u}}_1$ is the top eigenvector of ${\mathbf{M}}$. For the ${\mathbf{Y}}$ variable we initialize it with ${\mathbf{Y}}_1=\textrm{sign}({\mathbf{X}}_1)$ which is a subgradient of $\Vert{\mathbf{X}}_1\Vert_1$.

We set the step-size to $\eta=1/(2\lambda)$ 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.



\begin{table*}[!htb]\renewcommand{\arraystretch}{1.3}
{\footnotesize
\begin{center}
  \begin{tabular}{| l | c | c | c | c |} \hline 
  dimension (n) & $100$ & $200$ & $400$ & $600$ \\ \hline
  \multicolumn{5}{|c|}{$\downarrow$ ${\mathbf{N}}\sim U[0,1]$, $\textrm{SNR}=1$ $\downarrow$}\\ \hline
 
  $\lambda$ & $0.008$ & $0.004$ & $0.002$ & $0.0013$ \\ \hline
  initialization error & $0.5997$ & $0.6009$ & $0.5990$ & $0.6002$ \\ \hline
    recovery error & $0.0054$ & $0.0040$ & $0.0035$ & $0.0043$ \\ \hline
    dual gap & $4.1\times{10}^{-5}$ & $7.9\times{10}^{-5}$ & $4.9\times{10}^{-5}$ & $3.4\times{10}^{-6}$ \\ \hline
  $\lambda_{n-1}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))-\lambda_{n}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))$ & $0.8840$ & $0.8898$ & $0.8938$ & $0.8777$ \\ \hline 
  
  \multicolumn{5}{|c|}{$\downarrow$ ${\mathbf{N}}\sim U[0,1]$, $\textrm{SNR}=0.05$ $\downarrow$}\\ \hline
 
  $\lambda$ & $0.04$ & $0.02$ & $0.01$ & $0.0067$ \\ \hline
  initialization error & $1.7456$ & $1.7494$ & $1.7566$ & $1.7625$ \\ \hline
  recovery error& $0.0425$ & $0.0244$ & $0.0149$ & $0.0100$ \\ \hline
    dual gap & $2.0\times{10}^{-9}$ & $5.8\times{10}^{-6}$ & $4.5\times{10}^{-4}$ & $0.0018$ \\ \hline
  $\lambda_{n-1}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))-\lambda_{n}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))$ & $0.7092$ & $0.7854$ & $0.8340$ & $0.8622$ \\ \hline 
  
  \multicolumn{5}{|c|}{$\downarrow$ ${\mathbf{N}}\sim\mathcal{N}(0.5,{\mathbf{I}}_n)$, $\textrm{SNR}=1$ $\downarrow$}\\ \hline
 
  $\lambda$ & $0.006$ & $0.003$ & $0.0015$ & $0.001$ \\ \hline
  initialization error & $0.1584$ & $0.1464$ & $0.1443$ & $0.1411$ \\ \hline
  recovery error & $0.0059$ & $0.0033$ & $0.0019$ & $0.0015$ \\ \hline
    dual gap & $8.6\times{10}^{-4}$ & $0.0031$ & $0.0053$ & $0.0060$ \\ \hline
 
  $\lambda_{n-1}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))-\lambda_{n}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))$ & $0.8406$ & $0.8869$ & $0.9178$ & $0.9331$ \\ \hline 
  
  \multicolumn{5}{|c|}{$\downarrow$ ${\mathbf{N}}\sim\mathcal{N}(0.5,{\mathbf{I}}_n)$, $\textrm{SNR}=0.05$ $\downarrow$}\\ \hline
   
  $\lambda$ & $0.04$ & $0.02$ & $0.01$ & $0.005$ \\ \hline
  initialization error & $1.6701$ & $1.6620$ & $1.6542$ & $1.6610$ \\ \hline
  recovery error & $0.0502$ & $0.0234$ & $0.0137$ & $0.0109$ \\ \hline
    dual gap & $1.9\times{10}^{-5}$ & $0.0041$ & $0.0534$ & $0.0409$ \\ \hline
  $\lambda_{n-1}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))-\lambda_{n}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))$ & $0.2200$ & $0.4076$ & $0.5460$ & $0.6788$ \\ \hline
   \end{tabular}
  \caption{Numerical results for the sparse PCA problem.
  }\label{table:sparsePCAuniformBigNoise}
\end{center}
}
\vskip -0.2in
\end{table*}\renewcommand{\arraystretch}{1}






\subsection{Low-rank and sparse matrix recovery}

We consider the problem of recovering a simultaneously low-rank and sparse covariance matrix \cite{lowRankAndSparse}, which can be written as the following saddle-point optimization problem:
\begin{align*}
\min_{\substack{\textrm{Tr}({\mathbf{X}}) =1,\\ {\mathbf{X}}\succeq 0}}\frac{1}{2}\Vert{{\mathbf{X}}-{\mathbf{M}}}\Vert_F^2 + \lambda\Vert{{\mathbf{X}}}\Vert_1 =
\min_{\substack{\textrm{Tr}({\mathbf{X}}) =\tau,\\ {\mathbf{X}}\succeq 0}} \max_{\Vert{\mathbf{Y}}\Vert_{\infty}\le1}\frac{1}{2}\Vert{\mathbf{X}}-{\mathbf{M}}\Vert_F^2+\lambda\langle{\mathbf{X}},{\mathbf{Y}}\rangle,
\end{align*} 
where ${\mathbf{M}} = {{\mathbf{Z}}_0}{{\mathbf{Z}}_0}^{\top} + \frac{c}{2}({\mathbf{N}}+{\mathbf{N}}^{\top})$ is a noisy observation of some low-rank and sparse covariance matrix ${{\mathbf{Z}}_0}{{\mathbf{Z}}_0}^{\top}$. We choose ${\mathbf{Z}}_0\in\mathbb{R}^{n\times r}$ to be a sparse matrix where each entry ${{\mathbf{Z}}_0}_{i,j}$ is chosen to be $0$ with probability $0.9$ and $U\lbrace1,\ldots,10\rbrace$ with probability $0.1$, and then we normalize ${\mathbf{Z}}_0$ to be of unit Frobenius norm. We choose ${\mathbf{N}}\sim\mathcal{N}(0.5,{\mathbf{I}}_n)$.

We test the model with $\textrm{rank}({\mathbf{Z}}_0{\mathbf{Z}}_0^{\top})=1,5,10$.
We set the signal-to-noise ratio (SNR) to be a constant and set the noise level to $c=\frac{2\Vert{\mathbf{Z}}_0{\mathbf{Z}}_0^{\top}\Vert_F}{\textrm{SNR}\cdot\Vert{\mathbf{N}}+{\mathbf{N}}^{\top}\Vert_F}$ for our choice of $\textrm{SNR}$. 

We initialize the ${\mathbf{X}}$ variable with the rank-r approximation of ${\mathbf{M}}$. That is, we take ${\mathbf{X}}_1={\mathbf{U}}_r\textrm{diag}\left(\Pi_{\Delta_{\tau,r}}[\textrm{diag}(-\Lambda_r)]\right){\mathbf{U}}_r^{\top}$, where ${\mathbf{U}}_r\Lambda_r{\mathbf{U}}_r^{\top}$ is the rank-r eigen-decomposition of ${\mathbf{M}}$ and $\Delta_{\tau,r}=\lbrace{\mathbf{z}}\in\mathbb{R}^r~|~{\mathbf{z}}\ge0,\ \sum_{i=1}^r {\mathbf{z}}_i=\tau\rbrace$ is the simplex of radius $\tau$ in $\mathbb{R}^r$. For the ${\mathbf{Y}}$ variable we initialize it with ${\mathbf{Y}}_1=\textrm{sign}({\mathbf{X}}_1)$ which is a subgradient of $\Vert{\mathbf{X}}_1\Vert_1$.  

We set the step-size to $\eta=1$, $\tau=0.7\cdot\textrm{Tr}({\mathbf{Z}}_0{\mathbf{Z}}_0^{\top})$, and the number of iterations in each experiment to $T=2000$. For each value of $r$ and $n$ we average the measurements over over $10$ i.i.d. runs.

\begin{table*}[!htb]\renewcommand{\arraystretch}{1.3}
{\footnotesize
\begin{center}
  \begin{tabular}{| l | c | c | c | c | c | c | c | c | c |} \hline 
  dimension (n) & $100$ & $200$ & $400$ & $600$ \\ \hline
  
  \multicolumn{5}{|c|}{$\downarrow$ $r=\textrm{rank}({\mathbf{Z}}_0{\mathbf{Z}}_0^{\top})=1$, $\textrm{SNR}=0.48$ $\downarrow$}\\ \hline
  $\lambda$ & $0.0012$ & $0.0035$ & $0.0016$ & $0.001$ \\ \hline
  initialization error & $0.4562$ & $0.4471$ & $0.4507$ & $0.4450$ \\ \hline
    recovery error & $0.0364$ & $0.0193$ & $0.0160$ & $0.0168$ \\ \hline
    dual gap & $0.0083$ & $0.0086$ & $0.0020$ & $4.2\times{10}^{-4}$ \\ \hline
  $\lambda_{n-r}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))-\lambda_{n}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))$ & $0.0628$ & $0.1439$ & $0.1258$ & $0.1069$ \\ \hline

\multicolumn{5}{|c|}{$\downarrow$ $r=\textrm{rank}({\mathbf{Z}}_0{\mathbf{Z}}_0^{\top})=5$, $\textrm{SNR}=2.4$ $\downarrow$}\\ \hline
  $\lambda$ & $0.0012$ & $0.0006$ & $0.0003$ & $0.0002$ \\ \hline
  initialization error & $0.2132$ & $0.2103$ & $0.1983$ & $0.1907$ \\ \hline
    recovery error & $0.0641$ & $0.0478$ & $0.0349$ & $0.0274$ \\ \hline
    dual gap & $9.0\times{10}^{-4}$ & $4.3\times{10}^{-4}$ & $1.4\times{10}^{-4}$ & $7.3\times{10}^{-5}$ \\ \hline
  $\lambda_{n-r}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))-\lambda_{n}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))$ & $0.0148$ & $0.0200$ & $0.0257$ & $0.0277$ \\ \hline 

\multicolumn{5}{|c|}{$\downarrow$ $r=\textrm{rank}({\mathbf{Z}}_0{\mathbf{Z}}_0^{\top})=10$, $\textrm{SNR}=4.8$ $\downarrow$}\\ \hline
  $\lambda$ & $0.0007$ & $0.0004$ & $0.0002$ & $0.0001$ \\ \hline
  initialization error & $0.1855$ & $0.1661$ & $0.1527$ & $0.1473$ \\ \hline
    recovery error & $0.0702$ & $0.0403$ & $0.0268$ & $0.0356$ \\ \hline
    dual gap & $4.9\times{10}^{-4}$ & $6.6\times{10}^{-4}$ & $4.2\times{10}^{-4}$ & $3.4\times{10}^{-5}$ \\ \hline
  $\lambda_{n-r}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))-\lambda_{n}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))$ & $0.0072$ & $0.0142$ & $0.0187$ & $0.0160$ \\ \hline
   \end{tabular}
  \caption{Numerical results for the low-rank and sparse matrix recovery problem.
  }\label{table:lowRank&SparseRank10}
\end{center}
}
\vskip -0.2in
\end{table*}\renewcommand{\arraystretch}{1} 




\subsection{Robust PCA}

We consider the robust PCA problem \cite{robustPCA1} in the following formulation:
\begin{align*}
\min_{\substack{\textrm{Tr}({\mathbf{X}}) =\tau,\\ {\mathbf{X}}\succeq 0}}\Vert{{\mathbf{X}}-{\mathbf{M}}}\Vert_1 =
\min_{\substack{\textrm{Tr}({\mathbf{X}}) =\tau,\\ {\mathbf{X}}\succeq 0}} \max_{\Vert{\mathbf{Y}}\Vert_{\infty}\le1}\langle{\mathbf{X}}-{\mathbf{M}},{\mathbf{Y}}\rangle,
\end{align*}
where ${\mathbf{M}} = r{\mathbf{Z}}_0{\mathbf{Z}}_0^{\top} + \frac{1}{2}({\mathbf{N}}+{\mathbf{N}}^{\top})$ is a sparsely-corrupted observation of some rank-r  matrix ${\mathbf{Z}}_0{\mathbf{Z}}_0^{\top}$. We choose ${\mathbf{Z}}_0\in\mathbb{R}^{n\times r}$ to be a random unit Frobenius norm matrix. For ${\mathbf{N}}\in\mathbb{R}^{n\times n}$, we choose each entry to be $0$ with probability $1-1/\sqrt{n}$ and otherwise $1$ or $-1$ with equal probability.  
 

We initialize the ${\mathbf{X}}$ variable with the projection ${\mathbf{X}}_1=\Pi_{\lbrace\textrm{Tr}({\mathbf{X}})=\tau,\ {\mathbf{X}}\succeq 0\rbrace}[{\mathbf{M}}]$, and the ${\mathbf{Y}}$ variable with ${\mathbf{Y}}_1=\textrm{sign}({\mathbf{X}}_1-{\mathbf{M}})$. 


We test the model with $\textrm{rank}({\mathbf{Z}}_0{\mathbf{Z}}_0^{\top})=1,5,10$. For $\textrm{rank}({\mathbf{Z}}_0{\mathbf{Z}}_0^{\top})=1$ we set the step-size to $\eta=n/10$ and for $\textrm{rank}({\mathbf{Z}}_0{\mathbf{Z}}_0^{\top})=5,10$ we set it to $\eta=1$. We set the trace bound to $\tau=0.95\cdot\textrm{Tr}(r{\mathbf{Z}}_0{\mathbf{Z}}_0^{\top})$. For every set of parameters we average the measurements over $10$ i.i.d. runs. 



\begin{table*}[!htb]\renewcommand{\arraystretch}{1.3}
{\footnotesize
\begin{center}
  \begin{tabular}{| l | c | c | c | c | c | c | c | c | c |} \hline 
  dimension (n) & $100$ & $200$ & $400$ & $600$ \\ \hline
  
 \multicolumn{5}{|c|}{$\downarrow$ $r=\textrm{rank}({\mathbf{Z}}_0{\mathbf{Z}}_0^{\top})=1$, $T=3000$ $\downarrow$}\\ \hline 
SNR & $0.0021$ & $7.2\times{10}^{-4}$ & $2.5\times{10}^{-4}$ & $1.3\times{10}^{-4}$ \\ \hline
  initialization error & $1.3511$ & $1.3430$ & $1.2889$ & $1.2606$ \\ \hline
    recovery error & $0.0084$ & $0.0107$ & $0.0109$ & $0.0107$ \\ \hline
    dual gap & $0.0016$ & $0.0029$ & $0.0044$ & $0.0069$ \\ \hline
  $\lambda_{n-r}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))-\lambda_{n}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))$ & $15.5944$ & $41.2139$ & $85.8117$ & $140.5349$ \\ \hline 

\multicolumn{5}{|c|}{$\downarrow$ $r=\textrm{rank}({\mathbf{Z}}_0{\mathbf{Z}}_0^{\top})=5$, $T=20,000$ $\downarrow$}\\ \hline
SNR & $0.0110$ & $0.0038$ & $0.0013$ & $6.9\times{10}^{-4}$ \\ \hline
  initialization error & $1.5501$ & $1.5527$ & $1.5221$ & $1.4833$ \\ \hline
    recovery error & $0.0092$ & $0.0092$ & $0.0087$ & $0.0075$  \\ \hline
    dual gap & $0.0084$ & $0.0390$ & $0.1866$ & $0.4721$ \\ \hline
  $\lambda_{n-r}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))-\lambda_{n}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))$ & $7.6734$ & $26.2132$ & $66.1113$ & $108.7215$  \\ \hline 

\multicolumn{5}{|c|}{$\downarrow$ $r=\textrm{rank}({\mathbf{Z}}_0{\mathbf{Z}}_0^{\top})=10$, $T=30,000$ $\downarrow$}\\ \hline
  SNR & $0.0229$ & $0.0077$ & $0.0026$ & $0.0014$ \\ \hline
  initialization error & $1.5729$ & $1.6485$ & $1.6317$ & $1.5949$ \\ \hline
  recovery error & $0.0079$ & $0.0081$ & $0.0073$ & $0.0065$  \\ \hline
  dual gap & $0.0139$ & $0.0338$ & $0.1533$ & $0.3561$ \\ \hline
  $\lambda_{n-r}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))-\lambda_{n}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))$ & $1.7945$ & $16.9890$ & $48.9799$ & $82.2727$  \\ \hline 
   \end{tabular}
  \caption{Numerical results for the robust PCA problem.
  }\label{table:robustPCArank10}
\end{center}
}
\vskip -0.2in
\end{table*}\renewcommand{\arraystretch}{1}



\subsection{Phase synchronization}

We consider the phase synchronization problem (see for instance \cite{phaseSyncronization1}) which can be written as:
\begin{align} \label{phaseSync}
\max_{\substack{{\mathbf{z}}\in\mathbb{C}^n ,\\ \vert z_j\vert=1\ \forall j\in[n]}}{\mathbf{z}}^*{\mathbf{M}}{\mathbf{z}},
\end{align}
where ${\mathbf{M}} = {\mathbf{z}}_0{\mathbf{z}}_0^{*} + c{\mathbf{N}}$ is a noisy observation of some rank-one matrix ${\mathbf{z}}_0{\mathbf{z}}_0^{*}$ such that ${\mathbf{z}}_0\in\mathbb{C}^n$ and ${{\mathbf{z}}_0}_j=e^{i\theta_j}$ where $\theta_j\in[0,2\pi]$. We follow the statistical model in \cite{phaseSyncronization1} where the noise matrix ${\mathbf{N}}\in\mathbb{C}^{n\times{}n}$ is chosen such that every entry is
\begin{align*}
{\mathbf{N}}_{jk}=\Bigg\lbrace\begin{array}{lc}\mathcal{N}(0,1)+i\mathcal{N}(0,1) & j<k
\\ \overline{{\mathbf{N}}}_{kj} & j>k 
\\ 0 & j=k\end{array}.
\end{align*}
It is known that for a large $n$ and $c=\mathcal{O}\left(\sqrt{\frac{n}{\log{}n}}\right)$, with high probability the SDP relaxation of \eqref{phaseSync} is able to recover the original signal (see \cite{phaseSyncronization1}).

We solve a penalized version of the SDP relaxation of \eqref{phaseSync} which can be written as the following saddle-point optimization problem:
\begin{align*}
\min_{\substack{\textrm{Tr}({\mathbf{X}}) =n,\\ {\mathbf{X}}\succeq 0}}\langle{{\mathbf{X}},-{\mathbf{M}}}\rangle + \lambda\Vert{\textrm{diag}({\mathbf{X}})-\overrightarrow{\mathbf{1}}}\Vert_2 =
\min_{\substack{\textrm{Tr}({\mathbf{X}})=n,\\ {\mathbf{X}}\succeq 0}} \max_{\Vert{\mathbf{y}}\Vert_2\le1}\langle{\mathbf{X}},-{\mathbf{M}}\rangle+\lambda\langle\textrm{diag}({\mathbf{X}})-\overrightarrow{\mathbf{1}},{\mathbf{y}}\rangle,
\end{align*}
where $\overrightarrow{\mathbf{1}}$ is the all-ones vector.

While the phase synchronization problem is formulated over the complex numbers, extending our results to handle this model is straightforward.

We initialize the ${\mathbf{X}}$ variable with the rank-one approximation of ${\mathbf{M}}$. That is, we take ${\mathbf{X}}_1=n{\mathbf{u}}_1{\mathbf{u}}_1^{*}$, where ${\mathbf{u}}_1$ is the top eigenvector of ${\mathbf{M}}$. For the ${\mathbf{y}}$ variable we initialize it with ${\mathbf{y}}_1=(\textrm{diag}({\mathbf{X}}_1)-\overrightarrow{\mathbf{1}})/\Vert\textrm{diag}({\mathbf{X}}_1)-\overrightarrow{\mathbf{1}}\Vert_2$.

We set the noise level to $c=0.18\sqrt{n}$. We set the number of iterations in each experiment to $T=10,000$ and for each choice of $n$ we average the measurements over $10$ i.i.d. runs.

\begin{table*}[!htb]\renewcommand{\arraystretch}{1.3}
{\footnotesize
\begin{center}
  \begin{tabular}{| l | c | c | c | c | c | c | c | c | c |} \hline 
  dimension (n) & $100$ & $200$ & $400$ & $600$ \\ \hline
  SNR & $0.1553$ & $0.0775$ & $0.0387$ & $0.0258$ \\ \hline
$\lambda$ & $200$ & $600$ & $1600$ & $2800$ \\ \hline
$\eta$ & $1/400$ & $1/800$ & $1/1800$ & $1/1800$ \\ \hline
  initialization error & $0.1270$ & $0.1255$ & $0.1284$ & $0.1323$ \\ \hline
  recovery error & $0.0698$ & $0.0659$ & $0.0683$ & $0.0719$  \\ \hline
  dual gap & $7.8\times10^{-8}$ & $3.9\times10^{-5}$ & $0.1553$ & $0.5112$ \\ \hline
  $\lambda_{n-1}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))-\lambda_{n}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))$ & $39.8591$ & $78.9982$ & $150.3524$ & $217.06$  \\ \hline 
  $\Vert\textrm{diag}({\mathbf{X}}^*)-\overrightarrow{\mathbf{1}}\Vert_2$ & $3.2\times10^{-10}$ & $2.1\times10^{-8}$ & $5.1\times10^{-7}$ & $3.7\times10^{-7}$ \\ \hline
   \end{tabular}
  \caption{Numerical results for the phase synchronization problem.
  }\label{table:phaseSync}
\end{center}
}
\vskip -0.2in
\end{table*}\renewcommand{\arraystretch}{1}









\subsection{Linearly constrained low-rank matrix estimation}

Consider the following penalized formulation:
\begin{align*}
\min_{\substack{\textrm{Tr}({\mathbf{X}})=1,\\ {\mathbf{X}}\succeq 0}} \langle{\mathbf{X}},-{\mathbf{M}}\rangle+\lambda\Vert\mathcal{A}({\mathbf{X}})-{\mathbf{b}}\Vert_2 = \min_{\substack{\textrm{Tr}({\mathbf{X}})=1,\\ {\mathbf{X}}\succeq 0}} \max_{\Vert{\mathbf{y}}\Vert_2\le1}\langle{\mathbf{X}},-{\mathbf{M}}\rangle+\lambda\langle\mathcal{A}({\mathbf{X}})-{\mathbf{b}},{\mathbf{y}}\rangle,
\end{align*}

where ${\mathbf{M}} = {\mathbf{z}}_0{\mathbf{z}}_0^{\top} + \frac{c}{2}({\mathbf{N}}+{\mathbf{N}}^{\top})$ is the noisy observation of some rank-one matrix ${\mathbf{z}}_0{\mathbf{z}}_0^{\top}$ such that $\Vert{\mathbf{z}}_0\Vert_2=1$ and the noise matrix is chosen ${\mathbf{N}}\sim\mathcal{N}(0,{\mathbf{I}}_n)$. We take $\mathcal{A}({\mathbf{X}})=(\langle{\mathbf{A}}_1,{\mathbf{X}}\rangle,\ldots,\langle{\mathbf{A}}_m,{\mathbf{X}}\rangle)^{\top}$ with matrices ${\mathbf{A}}_1,\ldots,{\mathbf{A}}_m\in\mathbb{S}^n$ of the form ${\mathbf{A}}_i={\mathbf{v}}_i{\mathbf{v}}_i^{\top}$ such that ${\mathbf{v}}_i\sim\mathcal{N}(0,1)$. We take ${\mathbf{b}}\in\mathbb{R}^m$ such that $b_i=\langle{\mathbf{A}}_i,{\mathbf{z}}_0{\mathbf{z}}_0^{\top}\rangle$.

We initialize the ${\mathbf{X}}$ variable with the rank-one approximation of ${\mathbf{M}}$. That is, we take ${\mathbf{X}}_1={\mathbf{u}}_1{\mathbf{u}}_1^{\top}$, where ${\mathbf{u}}_1$ is the top eigenvector of ${\mathbf{M}}$. The ${\mathbf{y}}$ variable is initialized with ${\mathbf{y}}_1=(\mathcal{A}({\mathbf{X}}_1)-{\mathbf{b}})/\Vert\mathcal{A}({\mathbf{X}}_1)-{\mathbf{b}}\Vert_2$.

We set the number of constraints to $m=n$, the penalty parameter to $\lambda=2$, and the step-size to $\eta=1/(2\lambda)$. 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.


\begin{table*}[!htb]\renewcommand{\arraystretch}{1.3}
{\footnotesize
\begin{center}
  \begin{tabular}{| l | c | c | c | c | c | c | c | c | c |} \hline 
  dimension (n) & $100$ & $200$ & $400$ & $600$ \\ \hline
  $\textrm{SNR}$ & $0.15$ & $0.075$ & $0.04$ & $0.027$ \\ \hline
  initialization error & $0.1219$ & $0.1324$ & $0.1242$ & $0.1228$ \\ \hline
  recovery error & $0.0437$ & $0.0617$ & $0.0685$ & $0.0735$  \\ \hline
  dual gap & $5.3\times10^{-11}$ & $5.0\times10^{-12}$ & $8.5\times10^{-12}$ & $2.3\times10^{-11}$ \\ \hline
  $\lambda_{n-1}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))-\lambda_{n}(\nabla_{{\mathbf{X}}}f({\mathbf{X}}^*,{\mathbf{y}}^*))$ & $0.2941$ & $0.3409$ & $0.4690$ & $0.5069$  \\ \hline
  $\Vert\mathcal{A}({\mathbf{X}}^*)-{\mathbf{b}}\Vert_2$ & $0.0080$ & $0.0082$ & $0.0079$ & $0.0073$ \\ \hline
   \end{tabular}
  \caption{Numerical results for the linearly constrained low-rank matrix estimation problem. 
  }\label{table:linearlyConstrained}
\end{center}
}
\vskip -0.2in
\end{table*}\renewcommand{\arraystretch}{1}


\section{Discussion}
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  (\cref{lemma:negativeExampleNonSmooth}), 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.

\section*{Acknowledgements}
This research was supported by the ISRAEL SCIENCE FOUNDATION (grant No. 1108/18).

\bibliographystyle{plain}
