text:
\section{Introduction}

In this paper we consider the following optimization problem:
\begin{align} \label{eq:Model}
\min_{{\mathbf{X}}\in\mathcal{S}_n}{f({\mathbf{X}})},
\end{align}
where $\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$ - the space
of $n\times n$ real symmetric matrices, and $f:\mathbb{R}^{n\times n}\rightarrow\mathbb{R}$ is convex. For the sake of simplicity, throughout this paper we assume the gradient vector of $f(\cdot)$ is a symmetric matrix over $\mathbb{S}^n$, i.e., $\nabla{}f({\mathbf{X}})\in\mathbb{S}^n$ for all ${\mathbf{X}}\in\mathbb{S}^n$\footnote{In particular this holds for some of the most popular objective functions which take the form: $f({\mathbf{X}}) = g(\mathcal{A}{\mathbf{X}})+h({\mathbf{X}})$, where $\mathcal{A}:\mathbb{R}^{n\times n}\rightarrow\mathbb{R}^m$ is a linear map of the form $\mathcal{A}{\mathbf{X}} = (\textrm{Tr}({\mathbf{A}}_1^{\top}{\mathbf{X}}),\dots,\textrm{Tr}({\mathbf{A}}_m^{\top}{\mathbf{X}}))^{\top}$ for ${\mathbf{A}}_1,\dots,{\mathbf{A}}_m\in\mathbb{R}^{n\times n}$ (note that since ${\mathbf{X}}$ is symmetric, in case ${\mathbf{A}}_i$ is not symmetric we can always replace it with the symmetric matrix $({\mathbf{A}}_i+{\mathbf{A}}_i^{\top})/2$ without changing the result $\mathcal{A}{\mathbf{X}}$) and $g(\cdot)$ is differentiable over $\mathbb{R}^m$, and $h(\cdot)$ is of the form $h({\mathbf{X}})=w(\vec{\lambda}({\mathbf{X}}))$ for a symmetric and differentiable $w:\mathbb{R}^n\rightarrow\mathbb{R}$, where $\vec{\lambda}({\mathbf{X}})$ denotes the vector of eigenvalues of ${\mathbf{X}}$ (see \cite{Lewis1996}).}.
We also assume $f(\cdot)$ is $\beta$-smooth over $\mathcal{S}_n$ in the following typical sense:
$\Vert \nabla f({\mathbf{X}})-\nabla f({\mathbf{Y}})\Vert_2\le \beta\Vert {\mathbf{X}}-{\mathbf{Y}}\Vert_*$ for all ${\mathbf{X}},{\mathbf{Y}}\in\mathcal{S}_n$, where $\Vert{\cdot}\Vert_2$ is the spectral norm for matrices  (largest singular value) and $\Vert{\cdot}\Vert_*$ is the nuclear norm for matrices (sum of singular values).

The model \eqref{eq:Model} is a general convex framework for both deterministic and stochastic low-rank matrix recovery problems. These problems have many important modern applications in machine learning, statistics, signal processing, which include, among others, the well known tasks of matrix completion \cite{MatrixCompletion1,MatrixCompletion2,MatrixCompletion3,MatrixCompletion4}, robust PCA \cite{robustPCA1,robustPCA2,robustPCA3,robustPCA4,robustPCA5}, and phase retrieval \cite{phaseRetrieval1,phaseRetrieval2,phaseRetrieval3}.

We consider Problem \eqref{eq:Model} both in the deterministic first-order setting where the full gradient of $f(\cdot)$ is available, and in the stochastic first-order setting where only unbiased estimates for the gradient vector are available. In both settings, in terms of gradient complexity, the methods  of choice for solving Problem \eqref{eq:Model} in large scale are the non-Eulcidean proximal gradient methods that are based on Bregman distances w.r.t. the (negative) von Neumann entropy (entropy w.r.t. the eigenvalues) \cite{beckMD,sebastienMD,beckOptimizationBook, BregmanMatrices}.  This is because for Problem \eqref{eq:Model},  these methods are almost dimension independent \cite{sebastienMD}. In particular, as opposed to Euclidean methods whose complexities depend on the Lipschitz parameter of the gradient in Euclidean norm when using full gradients, or the Euclidean norm of the stochastic gradients and their variance in Euclidean norm when using stochastic gradients, proximal gradient methods based on the von Neumann entropy, measure these quantities w.r.t. to the spectral norm, which can often lead to convergence rates with significantly improved dependence on the dimension - a crucial consideration in high dimensional settings \footnote{For example, in the stochastic setting this can improve a factor of $\sqrt{n}$ in the rate of Euclidean SGD to only $\sqrt{\log{n}}$ for stochastic mirror descent with von Neumman entropy, or a factor of $n$ in the rate of Eucldiean projected gradient descent to $\log{n}$ for proximal gradient with Bregman distance w.r.t. the von Neumann entropy \cite{sebastienMD}.}. The family of proximal gradient methods with the von Neumann entropy for Problem \eqref{eq:Model} is also referred to in the literature as the Matrix Exponentiated Gradient (MEG) method (see for instance \cite{BregmanMatrices}), which we also adopt in this work. See also Algorithm \ref{alg:exactMEG} in the sequel for a full description of the method (both with deterministic and stochastic gradients).

Despite the improved gradient complexity of the MEG method over its Euclidean counterparts, a clear caveat is the high computational cost of executing each iteration, since each iteration requires in worst case a full-rank SVD computation of an $n\times n$ matrix  (to compute matrix logarithm and matrix exponential) which amounts to $O(n^3)$ runtime and is clearly not practical for high dimensional problems.

While the high computational cost of high-rank SVD computations is also present in Euclidean proximal gradient methods for Problem \eqref{eq:Model}, in the recent works \cite{garberLowRankSmooth, garberStochasticLowRank} it was established that under a natural strict complementarity condition (see discussion in the sequel) and at least in a certain ball around an optimal solution, Euclidean proximal gradient methods (both deterministic and stochastic) which only apply low-rank SVD computations (with rank that matches that of the optimal solution), provably converge to an optimal solution with rates similar to their full-SVD-based counterparts, resulting in computationally efficient methods, at least from a ``warm-start" initialization. This was also supported by empirical evidence brought in   \cite{garberLowRankSmooth, garberStochasticLowRank} which demonstrated the correct convergence of these methods using only low-rank computations in practice.

Thus, given that von Neumann entropy-based gradient methods are superior (in terms of gradient complexity) to their Euclidean counterparts for Problem \eqref{eq:Model}, and given the progress in \cite{garberLowRankSmooth, garberStochasticLowRank} on SVD-efficient Euclidean gradient methods, it is natural to ask whether similar results could also be obtained for the matrix exponentiated gradient method, and by that, to obtain the best of both worlds.

In this paper we provide a positive resolution to this question. We propose a variant of the MEG method which utilizes only low-rank SVD computations, and prove that under the strict complementary condition, it converges from a ``warm-start" initialization to an optimal solution, with similar rates to its full-SVD-based counterparts ($O(1/t)$ for deterministic gradients and $O(1/\sqrt{t})$ for stochastic gradients). Aside from our local convergence guarantees  we also suggest a recipe for computing certificates which can validate the correct convergence of our proposed low-rank MEG variant throughout the run. Finally, we present numerical simulations which supports both of our theory and overall methodology of incorporating low-rank updates to MEG and computing certificates for convergence.

As discussed, our provably local convergence guarantees hold under a strict complementarity condition which is also central to the analyses in \cite{garberLowRankSmooth, garberStochasticLowRank}, and also in the recent related work \cite{spectralFrankWolfe}. This condition is given in the following Assumption \ref{ass:strictcomp}. While, similarly to \cite{garberLowRankSmooth, garberStochasticLowRank}, we state this condition in terms of a certain eigen-gap in the gradient vector at optimal solutions which is natural for our analysis, in \cite{spectralFrankWolfe} it was showed that it is indeed equivalent to strict complementarity for Problem \eqref{eq:Model}. \footnote{This means that there exists a corresponding optimal solution to the dual problem ${\mathbf{Z}}^*\succeq 0$, such that $\textrm{rank}({\mathbf{X}}^*) + \textrm{rank}({\mathbf{Z}}^*)=n$.}

\begin{assumption}[strict complementarity]\label{ass:strictcomp}
An optimal solution ${\mathbf{X}}^*$  of rank $r^*$ for Problem \eqref{eq:Model} satisfies the strict complementarity assumption with parameter $\delta >0$ if $\lambda_{n-r^*}(\nabla f({\mathbf{X}}^*))-\lambda_{n}(\nabla f({\mathbf{X}}^*)) \geq \delta$.
\end{assumption}

Importantly, a positive eigen-gap between $\lambda_{n-r}(\nabla{}f({\mathbf{X}}^*))$ and $\lambda_n(\nabla{}f({\mathbf{X}}^*))$ for some optimal solution ${\mathbf{X}}^*$, as considered in Assumption \ref{ass:strictcomp}, is a sufficient condition (though not necessary)  for  ${\mathbf{X}}^*$ to have rank at most $r$. This is captured by the following lemma.

\begin{lemma}[Lemma 7 in \cite{garberLowRankSmooth}] \label{lemma:opt-grad}
Let ${\mathbf{X}}^*\in\mathcal{X^*}$ be any optimal solution, such that $\textrm{rank}({\mathbf{X}}^*)=r^*$, and write its eigen-decomposition as ${\mathbf{X}}^*=\sum_{i=1}^{r^*}{\lambda_i{\mathbf{v}}_i{\mathbf{v}}_i^{\top}}$. Then, the gradient vector $\nabla f({\mathbf{X}}^*)$ admits an eigen-decomposition such that the set of vectors $\{{\mathbf{v}}_i\}_{i=1}^{r^*}$ is a set of top eigen-vectors of $(-\nabla f({\mathbf{X}}^*))$ which corresponds to the eigenvalue $\lambda_1(-\nabla f({\mathbf{X}}^*))=-\lambda_n(\nabla f({\mathbf{X}}^*))$.
\end{lemma}

A motivation for the plausibility of Assumption \ref{ass:strictcomp} is that it implies the robustness of the low-rank of optimal solutions to Problem \eqref{eq:Model}, to small perturbations in the problem's parameters. See detailed discussions in \cite{garberLowRankSmooth, garberStochasticLowRank, spectralFrankWolfe}. This assumption is also supported by empirical experiments conducted in \cite{garberLowRankSmooth, garberStochasticLowRank, spectralFrankWolfe} and also in this current work.

Finally, while Problem \eqref{eq:Model} has received significant attention within the optimization, statistics, and machine learning communities in recent years with many attempts to provide efficient algorithms for large scale instances under a variety of assumptions, our current interest is to study efficient implementations of the classical proximal gradient methods with the von Neumann entropy because, at least in terms of gradient complexity, these are considered the methods of choice among first order methods. We refer the interested reader to \cite{garberLowRankSmooth, garberStochasticLowRank, spectralFrankWolfe} for comprehensive discussions about alternative methods/approaches and more related work.

The rest of this paper is organized as follows.
\begin{itemize}
\item
In Section \ref{sec:perm} we bring preliminaries including a fully detailed review of the matrix exponentiated gradient method. Importantly, in this section we also discuss the fundamental challenges of extending the results in \cite{garberLowRankSmooth, garberStochasticLowRank} which consider a Euclidean setup, to our von Neumann entropy-based Bregman distance setup.

\item
In Section \ref{sec:lowrankMEG} we present in full detail our low-rank MEG variant and give basic convergence results. These convergence results depend on how well the low-rank sequences approximates the steps of the standard MEG method. In this section we also present an efficient way to compute certificates that can verify that these approximations are well-behaved and by that certify to correct convergence of the method.

\item
In Section \ref{sec:main_tools} we present our main technical tool towards proving our main results: that in a certain ball around an optimal low-rank solution, under the strict complementarity assumption, our low-rank MEG variant indeed approximates sufficiently well the steps of the standard MEG method.

\item
In Section \ref{section:smoothCase} we prove our first main result:  local convergence of the low-rank MEG method with \textit{deterministic} updates with rate $O(1/t)$, under strict complementarity.

\item
In Section \ref{section:stochasticCase} we prove our second main result:  local convergence of the low-rank MEG method with \textit{stochastic} updates with rate $O(1/\sqrt{t})$, under strict complementarity.

\item
In Section \ref{section:experiments} we bring numerical experiments in support of our theoretical findings.
\end{itemize}

\section{Preliminaries}\label{sec:perm}

\subsection{Notation}

Throughout this work we use the following notations.
For real matrices we let $\Vert\cdot\Vert_2$ denote the spectral norm (i.e., the largest singular value), we let $\Vert\cdot\Vert_*$ denote the nuclear norm (i.e., the sum of all singular values),
and we let we let $\Vert\cdot\Vert_F$ denote the Frobenius norm.
For any ${\mathbf{X}},{\mathbf{Y}}\in\mathbb{S}^{n}$, we denote the standard matrix inner product as $\langle{\mathbf{X}},{\mathbf{Y}}\rangle={\mathbf{X}}\bullet{\mathbf{Y}}=\textrm{Tr}({\mathbf{X}}{\mathbf{Y}})$.
For a real symmetric matrix ${\mathbf{X}}\in\mathbb{S}^{n}$, we let
$\lambda_i({\mathbf{X}})$ denote its $i$th largest eigenvalue.
We denote $\mathbb{S}^n_{++}$ to be the set of all symmetric $n\times n$ positive-definite matrices.
For a matrix ${\mathbf{X}}\in\mathbb{S}^n_{++}$ with an eigen-decomposition ${\mathbf{X}}={\mathbf{V}} \mathrm{diag}(\lambda_1,\ldots,\lambda_n){\mathbf{V}}^{\top}$ we define the matrix logarithm as $\log({\mathbf{X}})={\mathbf{V}} \mathrm{diag}(\log(\lambda_1),\ldots,\log(\lambda_n)){\mathbf{V}}^{\top}$ and the matrix exponent as $\exp({\mathbf{X}})={\mathbf{V}} \mathrm{diag}(\exp(\lambda_1),\ldots,\exp(\lambda_n)){\mathbf{V}}^{\top}$.

\subsection{Bregman Distances, the von Neumann Entropy and the Matrix Exponentiated Gradient Algorithm}

\begin{definition}[Bregman distance]
Let $\omega$ be a real-valued, proper and $\alpha$-strongly convex function over a nonempty, closed and convex subset of the parameter domain that is continuously differentiable.
Then the Bregman distance is defined by: $B_{\omega}({\mathbf{x}},{\mathbf{y}})=\omega({\mathbf{x}})-\omega({\mathbf{y}})-\left<{\mathbf{x}}-{\mathbf{y}},\nabla\omega({\mathbf{y}})\right>$.
\end{definition}

For the problem under consideration in this paper, Problem \eqref{eq:Model}, where the parameter domain is the set of symmetric positive definite matrices, the standard strongly convex function that is considered is based on the von Neumann entropy, and is given by
\begin{align*
\omega({\mathbf{X}})=\textrm{Tr}({\mathbf{X}}\log({\mathbf{X}})-{\mathbf{X}}).
\end{align*}

The Bergman distance corresponding to the von Neumann entropy, has the form of
\begin{align*
B_{\omega}({\mathbf{X}},{\mathbf{Y}})& = \textrm{Tr}({\mathbf{X}}\log({\mathbf{X}})-{\mathbf{X}}\log({\mathbf{Y}})-{\mathbf{X}}+{\mathbf{Y}}).
\end{align*}

In this paper, we are interested in matrices in $\mathcal{S}_n$, for which the trace is equal to $1$. Therefore, for any ${\mathbf{X}},{\mathbf{Y}}\in\mathcal{S}_n$, the Bregman distance reduces to
\begin{align} \label{def:bregmanInSpectrahedron}
B_{\omega}({\mathbf{X}},{\mathbf{Y}})& = \textrm{Tr}({\mathbf{X}}\log({\mathbf{X}})-{\mathbf{X}}\log({\mathbf{Y}})).
\end{align}
In order for this definition to include symmetric positive semi-definite matrices, we use the convention $0\log(0):=0$.

For simplicity, from now on we will denote the Bergman distance corresponding to the von Neumann entropy as $B({\mathbf{X}},{\mathbf{Y}}):=B_{\omega}({\mathbf{X}},{\mathbf{Y}})$.

It is known that  the von Neumann entropy is $1$-strongly convex with respect to the nuclear-norm over $\mathcal{S}_n$. This implies that
\begin{align} \label{ineq:strongConvexityOfBregmanDistance}
B({\mathbf{X}},{\mathbf{Y}})\ge\frac{1}{2}\Vert {\mathbf{X}}-{\mathbf{Y}}\Vert_*^2.
\end{align}

Another important property of Bregman distances is the following three point identity (see \cite{beckMD}):
\begin{align} \label{lemma:threePointLemma}
B_{\omega}({\mathbf{X}},{\mathbf{Y}})+B_{\omega}({\mathbf{Y}},{\mathbf{Z}})-B_{\omega}({\mathbf{X}},{\mathbf{Z}})=\langle \nabla\omega({\mathbf{Z}})-\nabla\omega({\mathbf{Y}}),{\mathbf{X}}-{\mathbf{Y}}\rangle.
\end{align}

Bregman distances are central to Mirror-Decent algorithms, where the update step for some ${\mathbf{Z}}\in\mathcal{S}_n$ can be written, in our case, as:
\begin{align} \label{eq:MDGeneralUpdate}
 {\mathbf{Z}}_{+}&=\argmin_{{\mathbf{X}}\in\mathcal{S}_n}\left\{\langle\eta\nabla f({\mathbf{Z}}),{\mathbf{X}}\rangle+B({\mathbf{X}},{\mathbf{Z}})\right\} \nonumber \\
 &=\argmin_{{\mathbf{X}}\in\mathcal{S}_n}\left\{\langle\eta\nabla f({\mathbf{Z}})-\nabla\omega({\mathbf{Z}}),{\mathbf{X}}\rangle+\omega({\mathbf{X}})\right\},
\end{align}
where $\eta > 0$ is the step-size.

In the case of the von Neumann entropy, this update is equivalent to an update of the following form (see \cite{BregmanMatrices}):
\begin{align} \label{eq:MDvonNeumannUpdate}
{\mathbf{Z}}_{+}=\frac{1}{b}\exp(\log({\mathbf{Z}})-\eta\nabla f({\mathbf{Z}})),
\end{align}
where $b:=\textrm{Tr}(\exp(\log({\mathbf{Z}})-\eta\nabla f({\mathbf{Z}})))$.

For stochastic optimization, the gradient is replaced with an unbiased estimator of the gradient, $\widehat{\nabla}\in\mathbb{S}^n$, resulting in the following update:
  \begin{align} \label{eq:MDvonNeumannStochasticUpdate}
{\mathbf{Z}}_{+}=\frac{1}{b}\exp(\log({\mathbf{Z}})-\eta\widehat{\nabla}),
\end{align}
where $b:=\textrm{Tr}(\exp(\log({\mathbf{Z}})-\eta\widehat{\nabla}))$.

The update rules \eqref{eq:MDvonNeumannUpdate} and \eqref{eq:MDvonNeumannStochasticUpdate} are simply the Matrix Exponentiated Gradient (MEG) method and its stochastic counterpart \cite{BregmanMatrices}, which are given in complete form in Algorithm \ref{alg:exactMEG}.

\begin{algorithm}
	\caption{Matrix Exponentiated Gradient for Problem \eqref{eq:Model}}\label{alg:MEG}
	\label{alg:exactMEG}
	\begin{algorithmic}
		\STATE \textbf{Input:} $\{\eta_t\}_{t\geq 1}$ - sequence of (positive) step-sizes
		\STATE ${\mathbf{Z}}_1$ - an arbitrary point in $\mathcal{S}_{n}\cap\mathbb{S}^n_{++}$
		\FOR{$t = 1,2,...$}
            \STATE \ $ {\mathbf{Y}}_t=\exp(\log ({\mathbf{Z}}_{t})-\eta_{t}\nabla f({\mathbf{Z}}_{t}))$ (for deterministic optimization)
            \\ \OR
            \STATE \ $ {\mathbf{Y}}_t=\exp(\log ({\mathbf{Z}}_{t})-\eta_{t}\widehat{\nabla}_t)$ (for stochastic optimization)
			\STATE $b_{t}=\textrm{Tr}({\mathbf{Y}}_t)$
			\STATE $ {\mathbf{Z}}_{t+1} = {\mathbf{Y}}_{t}/b_t$
        \ENDFOR
	\end{algorithmic}
\end{algorithm}

The following result is standard and states the convergence rates of Algorithm \ref{alg:exactMEG} with either full or stochastic gradients (see for example \cite{beckOptimizationBook,sebastienMD}).
\begin{theorem}[Convergence of Algorithm \ref{alg:exactMEG}]
When used with deterministic gradients and fixed step-size $\eta = 1/\beta$, the iterates $\{{\mathbf{X}}_t\}_{t\geq 1}$ of Algorithm \ref{alg:exactMEG} satisfy: $f({\mathbf{X}}_t) -f^* = O\left({\frac{\beta\log{n}}{t}}\right)$ for all $t\geq 1$. Alternatively, when used with stochastic gradients and step-size $\eta_t = \frac{\sqrt{\log{n}}}{G\sqrt{t}}$, the ergodic sequence $\{\bar{{\mathbf{X}}}_t\}_{t\geq 1}$, where $\bar{{\mathbf{X}}}_t := \frac{1}{t}\sum_{i=1}^t{\mathbf{X}}_t$, satisfy: $\mathbb{E}[f(\bar{{\mathbf{X}}}_t) - f^*] = O\left({\frac{G\sqrt{\log{n}}}{\sqrt{t}}}\right)$ for all $t\geq 1$, where $G \geq \max_{t}\Vert{\widehat{\nabla}_t}\Vert_2$.
\end{theorem}

It is important to note that from a computational perspective, the most expensive step in Algorithm \ref{alg:exactMEG} is the computation of the matrix exponential in the update of the matrix variable ${\mathbf{Y}}_t$ on each iteration. Computing a matrix exponent requires in worst case a full-rank SVD of the input matrix which takes $O(n^3)$ runtime.

\subsection{The challenge of low-rank MEG updates}\label{sec:challenges}

Since our work is mostly inspired by the recent works \cite{garberLowRankSmooth, garberStochasticLowRank} which developed similar local convergence results for Problem \eqref{eq:Model} using Euclidean gradient methods that only apply low-rank SVD computations, we first discuss the differences between the Euclidean and our non-Euclidean Bregman setups, and we emphasize and clarify the difficulty of extending the results \cite{garberLowRankSmooth, garberStochasticLowRank} to the von Neumann entropy-based Bregman setup.

We first recall the Euclidean projected gradient mapping for the spectrahedron $\mathcal{S}_n$, given by $\Pi_{\mathcal{S}_n}[{\mathbf{X}}-\eta\nabla{}f({\mathbf{X}})]$ for some ${\mathbf{X}}\in\mathcal{S}_n$, where $\Pi_{\mathcal{S}_n}[\cdot]$ denotes the Euclidean projection onto $\mathcal{S}_n$. Let ${\mathbf{Y}}=\sum_{i=1}^n\lambda_i{\mathbf{v}}_i{\mathbf{v}}_i^{\top}$ denote the eigen-decomposition of the matrix ${\mathbf{Y}}={\mathbf{X}}-\eta\nabla{}f({\mathbf{X}})$. Its Euclidean projection onto the spectrahedron is given by
\begin{align}\label{eq:euclidProj}
\Pi_{\mathcal{S}_n}[{\mathbf{Y}}]=\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$.

We can see that the Euclidean projection has a thresholding effect on the eigenvalues, i.e., all eigenvalues below or equal to $\lambda$ become zero after the projection operation. This unique property is central to the analysis in \cite{garberLowRankSmooth, garberStochasticLowRank} and is used to establish that in the proximity of an optimal solution which satisfies an eigen-gap assumption such as the strict complementarity condition (Assumption \ref{ass:strictcomp}), the Euclidean projected gradient mapping always results in a low-rank matrix, and thus, as evident from \eqref{eq:euclidProj}, can be computed exactly using only a low-rank SVD computation.

This unique property and consequence of the Euclidean projected gradient mapping does not hold anymore for the von Neumann entropy-based MEG mapping given in \eqref{eq:MDvonNeumannUpdate}. In particular we see that almost everywhere in $\mathcal{S}_n$ (except for the points for which this mapping results in a low-rank matrix, such as low-rank optimal solutions), due to the matrix logarithm and exponent, this mapping is defined only for full-rank matrices and also results in full-rank matrices. Thus, in stark contrast to the Euclidean case, we cannot expect that in any proximity of a low-rank optimal solution that this mapping could be computed exactly via a low-rank SVD. Thus, our last hope is that the MEG mapping \eqref{eq:MDvonNeumannUpdate} could be sufficiently approximated (in a suitable sense) using only a low-rank SVD computation.

Thus, when attempting to derive results in the spirit of  \cite{garberLowRankSmooth, garberStochasticLowRank} for the MEG method, our challenge is two-folded: i. we need to suggest a way to approximate the steps of MEG, which generally requires to store and manipulate full-rank matrices, using only low-rank SVD computations, and ii. we need to be able to establish that the resulting approximation errors can indeed be properly controlled to guarantee the convergence of the new approximated updates.

\section{The Approximated MEG Method with Low-Rank Updates}\label{sec:lowrankMEG}

Our approach to replacing the exact MEG updates in Algorithm \ref{alg:MEG} with approximated low-rank updates is fairly straightforward. Instead of computing the mapping ${\mathbf{Y}} = \exp(\log{\mathbf{X}}-\eta\nabla{}f({\mathbf{X}}))$ for some ${\mathbf{X}}\in\mathcal{S}_n$ (e.g., in case of deterministic updates), we compute explicitly only its best rank-$r$ approximation ${\mathbf{Y}}^r$ for some parameter $r$, which is simply given by the top $r$ components in the eigen-decomposition of ${\mathbf{Y}}$, and then slightly perturb it in a way that sets all eigenvalues in the places $r+1$,\dots,$n$ to some small fixed value, so that we end up with a positive definite matrix. We first show that such consecutive updates could indeed be performed by computing only a rank-$r$ SVD and that the iterates of the proposed method could be stored in memory in the form of low-rank matrices.

Concretely, focusing on deterministic updates (the stochastic case follows the same reasoning with the obvious changes), we define the low-rank MEG update to some positive-definite matrix ${\mathbf{X}}\in\mathcal{S}_n$, as
\begin{align} \label{eq:LowRankMDvonNeumannUpdate}
{\mathbf{X}}_+:=(1-\varepsilon)\frac{{\mathbf{Y}}^r}{a}+\frac{\varepsilon}{n-r}({\mathbf{I}}-{\mathbf{V}}^r{{\mathbf{V}}^r}^{\top})
\end{align}
for some $\varepsilon\in[0,1]$, where ${\mathbf{Y}}:=\exp(\log({\mathbf{X}})-\eta\nabla f({\mathbf{X}}))$ admits the  eigen-decomposition ${\mathbf{Y}}={\mathbf{V}}\Lambda{\mathbf{V}}^{\top}$, and ${\mathbf{Y}}^r={\mathbf{V}}^r\Lambda^r{{\mathbf{V}}^r}^{\top}$ is its rank-$r$ approximation taken to be the top $r$ components in the eigen-decomposition and $a:=\sum_{i=1}^{r} \lambda_i({\mathbf{Y}})$ is its trace.

It is important to note that in case the matrix $\log({\mathbf{X}})$ is readily given, the matrices ${\mathbf{Y}}^r$ and ${\mathbf{V}}^r$ can indeed be calculated in an efficient manner  by computing the top-$r$ components in the eigen-decomposition of $\log({\mathbf{X}})-\eta\nabla f({\mathbf{X}})$ and then applying the matrix exponent to the resulting matrix to obtain ${\mathbf{Y}}^r$ \footnote{Note that applying the matrix exponent to a matrix given by its eigen-decomposition simply requires to exponentiate the eigenvalues since the eigen-vectors remain unchanged.}.

Now suppose ${\mathbf{X}}$  itself has a similar structure to that in \eqref{eq:LowRankMDvonNeumannUpdate} that is,
${\mathbf{X}}=(1-\varepsilon_{-})\tilde{{\mathbf{X}}}+\frac{\varepsilon_{-}}{n-r}({\mathbf{I}}-{\mathbf{V}}_{-}^r{{\mathbf{V}}_{-}^r}^{\top})$,
where $\tilde{{\mathbf{X}}}\in\mathcal{S}_n$ is such that $\textrm{rank}(\tilde{{\mathbf{X}}})=r$ and $\varepsilon_{-}\in[0,1]$, and the eigen-decomposition of $\tilde{{\mathbf{X}}}$ is $\tilde{{\mathbf{X}}}={\mathbf{V}}_{-}^r\Lambda_{-}^r{{\mathbf{V}}_{-}^r}^{\top}$, where ${\mathbf{V}}_-^r\in\mathbb{R}^{n\times r}, \Lambda_-^r\in\mathbb{R}^{r\times r}$. Then the eigen-decomposition of ${\mathbf{X}}$ can be written as
\begin{align}\label{eq:eigenFastUpdate}
{\mathbf{X}}&=(1-\varepsilon_{-}){\mathbf{V}}_{-}^r\Lambda_{-}^r{{\mathbf{V}}_{-}^r}^{\top}+\frac{\varepsilon_{-}}{n-r}\left({\mathbf{I}}-{\mathbf{V}}_{-}^r{{\mathbf{V}}_{-}^r}^{\top}\right).
\end{align}
Thus,
\begin{align*}
\log({\mathbf{X}})={\mathbf{V}}_{-}^r\log\left((1-\varepsilon_{-})\Lambda_{-}^r\right){{\mathbf{V}}_{-}^r}^{\top}+\log\left(\frac{\varepsilon_{-}}{n-r}\right)\left({\mathbf{I}}-{\mathbf{V}}_{-}^r{{\mathbf{V}}_{-}^r}^{\top}\right).
\end{align*}
In particular, when ${\mathbf{X}}$ admits such structure, with the eigen-decomposition of $\tilde{{\mathbf{X}}}$ explicitly given, $\log({\mathbf{X}})$ can be computed in $O(rn^2)$ time, and consequently, the update in \eqref{eq:LowRankMDvonNeumannUpdate} can be computed using a single rank-$r$ SVD and additional $O(rn^2)$ time.

Moreover, computing the matrices ${\mathbf{Y}}^r,{\mathbf{V}}^r$ via fast iterative methods is usually based on iteratively applying some orthonormal $n\times r$ matrix ${\mathbf{U}}$ to the matrix $\log({\mathbf{X}})-\eta\nabla{}f({\mathbf{X}})$ (see for instance \cite{Musco15}), which as discussed, under the above assumption on the structure of ${\mathbf{X}}$ takes the form:
\begin{align}\label{eq:QRiterations}
&\left({\log({\mathbf{X}})-\eta\nabla{}f({\mathbf{X}})}\right){\mathbf{U}} \nonumber= \\
& \Big({\mathbf{V}}_{-}^r\log\left((1-\varepsilon_{-})\Lambda_{-}^r\right){{\mathbf{V}}_{-}^r}^{\top}+\log\left(\frac{\varepsilon_{-}}{n-r}\right)\left({\mathbf{I}}-{\mathbf{V}}_{-}^r{{\mathbf{V}}_{-}^r}^{\top}\right)- \eta\nabla{}f({\mathbf{X}})\Big){\mathbf{U}} =\nonumber \\
& {\mathbf{V}}_{-}^r\log\left((1-\varepsilon_{-})\Lambda_{-}^r\right){{\mathbf{V}}_{-}^r}^{\top}{\mathbf{U}}+\log\left(\frac{\varepsilon_{-}}{n-r}\right)\left({\mathbf{U}}-{\mathbf{V}}_{-}^r{{\mathbf{V}}_{-}^r}^{\top}{\mathbf{U}}\right)- \eta\nabla{}f({\mathbf{X}}){\mathbf{U}}.
\end{align}
Examining the various matrix multiplications in the RHS of \eqref{eq:QRiterations}, we can see that the matrix $\left({\log({\mathbf{X}})-\eta\nabla{}f({\mathbf{X}})}\right){\mathbf{U}}$ can be computed in $O(r^2n + r\cdot \textrm{nnz}(\nabla{}f({\mathbf{X}}))$ time, which can be much more efficient than $O(rn^2)$ for problems in which the gradient is typically sparse.
Moreover, from the RHS of \eqref{eq:QRiterations} we can see that in order to apply the iterations of fast iterative SVD methods towards the computation of the new iterate ${\mathbf{X}}_+$, we only require to store the rank-$r$ eigen-decomposition of $\tilde{{\mathbf{X}}}$, which amounts to storing only a rank-$r$ matrix in memory (the matrix ${\mathbf{V}}^r_{-}$) and the $r$ eigenvalues (diagonal entries of $\Lambda^r_{-}$).

Note that while the first iterate ${\mathbf{X}}_1$ in Algorithm \ref{alg:LRMD} takes a slightly different form than Eq. \eqref{eq:eigenFastUpdate}, as long as the initialization matrix ${\mathbf{X}}_0$ is of rank $r$ and its eigen-decomposition is given, writing it in the form of \eqref{eq:eigenFastUpdate} is straightforward.

The complete description of our low-rank matrix exponentiated gradient algorithm is given below as Algorithm \ref{alg:LRMD}. To conclude this part, when using fast-iterative methods to compute the eigen-decomposition of ${\mathbf{Y}}_r$, aside from the computation of the gradient vector, our algorithm only requires to store matrices of rank at most $r$  (which requires storing only $r(n+1)$ scalars when the eigen-decomposition is explicitly maintained).

\begin{algorithm}
	\caption{Low Rank Matrix Exponentiated Gradient for Problem \eqref{eq:Model}}
	\label{alg:LRMD}
	\begin{algorithmic}
		\STATE \textbf{Input:} $\{\eta_t\}_{t\ge1}$ - sequence of (positive) step-sizes, $\lbrace\varepsilon_t\rbrace_{t\ge 0}\subset[0,1]$ - sequence of approximation parameters, $r\in\{1,2,\dots,n-1\}$ - SVD rank parameter, initialization matrix -  ${\mathbf{X}}_0\in\mathcal{S}_n$ such that $\textrm{rank}({\mathbf{X}}_0) = r$
		\STATE \textbf{Initialization:} ${\mathbf{X}}_1=(1-\varepsilon_0){\mathbf{X}}_0+\frac{\varepsilon_0}{n}{\mathbf{I}}$
		\FOR{$t = 1,2,...$}
            \STATE \ $ {\mathbf{Y}}_t=\exp(\log ({\mathbf{X}}_{t})-\eta_{t}\nabla f({\mathbf{X}}_{t}))$ (for deterministic optimization; not to be explicitly computed)

            \\ \OR
            \STATE \ $ {\mathbf{Y}}_t=\exp(\log ({\mathbf{X}}_{t})-\eta_{t}\widehat{\nabla}_{t})$ (for stochastic optimization;  not to be explicitly computed) \COMMENT{$\widehat{\nabla}_{t}$ is unbiased estimator for  $\nabla{}f({\mathbf{X}}_t)$}

            \STATE ${\mathbf{Y}}_t^r={\mathbf{V}}^r\Lambda^r{{\mathbf{V}}^r}^{\top}$ - rank-$r$ eigen-decomposition of ${\mathbf{Y}}_t$
			\STATE $a_{t}=\textrm{Tr}({\mathbf{Y}}_t^r)$
            \STATE ${\mathbf{X}}_{t+1}=(1-\varepsilon_t)\frac{{\mathbf{Y}}_t^r}{a_t}+\frac{\varepsilon_t}{n-r}({\mathbf{I}}-{\mathbf{V}}^r{{\mathbf{V}}^r}^{\top})$
        \ENDFOR
	\end{algorithmic}
\end{algorithm}

\subsection{Convergence of approximated sequences and computing certificates}\label{section:localConvergence}
As a starting point for our convergence analysis of Algorithm \ref{alg:LRMD}, we first state and prove the most general convergence results which do not rely on any additional assumption (such as strict complementarity or ``warm-start'' initialization). These results, which will also serve as the basis for all of our following theoretical derivations, naturally depend on how well the low-rank updates applied in Algorithm \ref{alg:LRMD} approximate their exact counterparts in Algorithm \ref{alg:exactMEG}. While without further assumptions it is not possible to a-priori guarantee much for these approximations, we conclude this section by providing computable certificates that can ensure during runtime that these approximations are indeed properly bounded, and hence, can be used to certify the correct convergence of the method in practice.

We note that the following two convergence theorems are quite generic and apply to any inexact update sequence $\lbrace {\mathbf{X}}_t\rbrace_{t\ge 1}$ which approximates the steps of the exact MEG method (Algorithm \ref{alg:exactMEG}), and do not directly rely on the updates applied in Algorithm \ref{alg:LRMD}.

\begin{theorem}\label{thm:converganceSmooth}[Convergence with deterministic updates]
Let $\lbrace {\mathbf{X}}_t\rbrace_{t\geq 1}$ be a sequence of points in $\mathcal{S}_n\cap\mathbb{S}^n_{++}$. Let $\{{\mathbf{Z}}_t\}_{t\geq 1}$ be a sequence such that ${\mathbf{Z}}_1={\mathbf{X}}_1$, and for all $t\geq 1$, ${\mathbf{Z}}_{t+1}$ is the accurate MEG update to ${\mathbf{X}}_t$ as defined in \eqref{eq:MDvonNeumannUpdate} with step-size $\eta_t=\eta\le\frac{1}{\beta}$. Then, for any $T\geq 1$ we have that
\begin{align*}
\min_{1\le t\le {T}}f({\mathbf{X}}_t) -f({\mathbf{X}}^*) & \le \frac{1}{\eta T}B({\mathbf{X}}^*,{\mathbf{X}}_1)+\frac{1}{\eta T}\sum_{t=1}^{T-1}\left(B({\mathbf{X}}^*,{\mathbf{X}}_{t+1})-B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})\right)
\\ & \ \ \ \ +\frac{2}{T}\sum_{t=1}^{T-1}\Vert\nabla f({\mathbf{X}}_{t+1})\Vert_{2}\sqrt{B({\mathbf{Z}}_{t+1},{\mathbf{X}}_{t+1})}.
\end{align*}
\end{theorem}

\begin{proof}
By the definition of ${\mathbf{Z}}_{t+1}$ in \eqref{eq:MDGeneralUpdate},
\[ {\mathbf{Z}}_{t+1}=\argmin_{{\mathbf{X}}\in\mathcal{S}_n}\{\langle\eta_t\nabla f({\mathbf{X}}_t)-\nabla\omega({\mathbf{X}}_t),{\mathbf{X}}\rangle+\omega({\mathbf{X}})\}. \]

Therefore, by the optimality condition for ${\mathbf{Z}}_{t+1}$, $\forall {\mathbf{X}}\in\mathcal{S}_n$,
\[ \langle \eta_t\nabla f({\mathbf{X}}_t)-\nabla\omega({\mathbf{X}}_t)+\nabla\omega({\mathbf{Z}}_{t+1}),{\mathbf{X}}-{\mathbf{Z}}_{t+1}\rangle\ge0. \]

Rearranging and using the three point lemma in \eqref{lemma:threePointLemma}, we get that $\forall {\mathbf{X}}\in\mathcal{S}_n$,
\begin{align} \label{ineq:smoothOptZ(t+1)}
\eta_t\langle\nabla f({\mathbf{X}}_t),{\mathbf{Z}}_{t+1}-{\mathbf{X}}\rangle & \le \langle\nabla\omega({\mathbf{X}}_t)-\nabla\omega({\mathbf{Z}}_{t+1})),{\mathbf{Z}}_{t+1}-{\mathbf{X}}\rangle \nonumber
\\ & = B({\mathbf{X}},{\mathbf{X}}_t)-B({\mathbf{X}},{\mathbf{Z}}_{t+1})-B({\mathbf{Z}}_{t+1},{\mathbf{X}}_t).
\end{align}

Using the $\beta$-smoothness of $f$, the gradient inequality and the inequality in \eqref{ineq:strongConvexityOfBregmanDistance}, we obtain
\begin{align*}
f({\mathbf{Z}}_{t+1}) & \le f({\mathbf{X}}_{t})+\langle\nabla f({\mathbf{X}}_t),{\mathbf{Z}}_{t+1}-{\mathbf{X}}_t\rangle+\frac{\beta}{2}\Vert{\mathbf{Z}}_{t+1}-{\mathbf{X}}_t\Vert_*^2
\\ & \le f({\mathbf{X}}^*)+\langle\nabla f({\mathbf{X}}_t),{\mathbf{Z}}_{t+1}-{\mathbf{X}}^*\rangle+\frac{\beta}{2}\Vert{\mathbf{Z}}_{t+1}-{\mathbf{X}}_t\Vert_*^2
\\ & \le f({\mathbf{X}}^*)+\langle\nabla f({\mathbf{X}}_t),{\mathbf{Z}}_{t+1}-{\mathbf{X}}^*\rangle+\betaB({\mathbf{Z}}_{t+1},{\mathbf{X}}_t).
\end{align*}

In addition, from the gradient inequality
\begin{align*}
f({\mathbf{Z}}_{t+1}) & \ge f({\mathbf{X}}_{t+1})-\langle\nabla f({\mathbf{X}}_{t+1}),{\mathbf{X}}_{t+1}-{\mathbf{Z}}_{t+1}\rangle.
\end{align*}

Combining the last two inequalities gives us
\begin{align*}
f({\mathbf{X}}_{t+1}) & \le f({\mathbf{X}}^*)+\langle\nabla f({\mathbf{X}}_t),{\mathbf{Z}}_{t+1}-{\mathbf{X}}^*\rangle+\betaB({\mathbf{Z}}_{t+1},{\mathbf{X}}_t)+\langle\nabla f({\mathbf{X}}_{t+1}),{\mathbf{X}}_{t+1}-{\mathbf{Z}}_{t+1}\rangle.
\end{align*}

Plugging in \eqref{ineq:smoothOptZ(t+1)} with ${\mathbf{X}}={\mathbf{X}}^*$ we get
\begin{align*}
f({\mathbf{X}}_{t+1}) & \le f({\mathbf{X}}^*)+\frac{1}{\eta_t}\left(B({\mathbf{X}}^*,{\mathbf{X}}_t)-B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})-B({\mathbf{Z}}_{t+1},{\mathbf{X}}_t)\right)
\\ & \ \ \ \ +\betaB({\mathbf{Z}}_{t+1},{\mathbf{X}}_t)+\langle\nabla f({\mathbf{X}}_{t+1}),{\mathbf{X}}_{t+1}-{\mathbf{Z}}_{t+1}\rangle
\\ & \le f({\mathbf{X}}^*)+\frac{1}{\eta_t}\left(B({\mathbf{X}}^*,{\mathbf{X}}_t)-B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})\right)+\left(\beta-\frac{1}{\eta_t}\right)B({\mathbf{Z}}_{t+1},{\mathbf{X}}_t)
\\ & \ \ \ \ +\Vert\nabla f({\mathbf{X}}_{t+1})\Vert_{2}\cdot\Vert{\mathbf{X}}_{t+1}-{\mathbf{Z}}_{t+1}\Vert_*
\\ & = f({\mathbf{X}}^*)+\frac{1}{\eta_t}\left(B({\mathbf{X}}^*,{\mathbf{X}}_t)-B({\mathbf{X}}^*,{\mathbf{X}}_{t+1})\right)+\left(\beta-\frac{1}{\eta_t}\right)B({\mathbf{Z}}_{t+1},{\mathbf{X}}_t)
\\ & \ \ \ \ +\frac{1}{\eta_t}\left(B({\mathbf{X}}^*,{\mathbf{X}}_{t+1})-B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})\right)+\Vert\nabla f({\mathbf{X}}_{t+1})\Vert_{2}\cdot\Vert{\mathbf{X}}_{t+1}-{\mathbf{Z}}_{t+1}\Vert_*
\\ & \le f({\mathbf{X}}^*)+\frac{1}{\eta_t}\left(B({\mathbf{X}}^*,{\mathbf{X}}_t)-B({\mathbf{X}}^*,{\mathbf{X}}_{t+1})\right)+\left(\beta-\frac{1}{\eta_t}\right)B({\mathbf{Z}}_{t+1},{\mathbf{X}}_t)
\\ & \ \ \ \ +\frac{1}{\eta_t}\left(B({\mathbf{X}}^*,{\mathbf{X}}_{t+1})-B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})\right)+2\Vert\nabla f({\mathbf{X}}_{t+1})\Vert_{2}\sqrt{B({\mathbf{Z}}_{t+1},{\mathbf{X}}_{t+1})},
\end{align*}
where the last inequality follows from \eqref{ineq:strongConvexityOfBregmanDistance}.

Averaging over $t=0,\ldots,T-1$ and taking $\eta_t=\eta\le\frac{1}{\beta}$ we reach,
\begin{align*}
\min_{1\le t\le {T}}f({\mathbf{X}}_t) -f({\mathbf{X}}^*) & \le  \frac{1}{T}\sum_{t=0}^{T-1} f({\mathbf{X}}_{t+1}) -f({\mathbf{X}}^*)
\\ & \le \frac{1}{\eta T}B({\mathbf{X}}^*,{\mathbf{X}}_1)+\frac{1}{\eta T}\sum_{t=0}^{T-1}\left(B({\mathbf{X}}^*,{\mathbf{X}}_{t+1})-B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})\right)
\\ & \ \ \ \ +\frac{2}{T\alpha}\sum_{t=0}^{T-1}\Vert\nabla f({\mathbf{X}}_{t+1})\Vert_{2}\sqrt{B({\mathbf{Z}}_{t+1},{\mathbf{X}}_{t+1})}
\\ & = \frac{1}{\eta T}B({\mathbf{X}}^*,{\mathbf{X}}_1)+\frac{1}{\eta T}\sum_{t=1}^{T-1}\left(B({\mathbf{X}}^*,{\mathbf{X}}_{t+1})-B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})\right)
\\ & \ \ \ \ +\frac{2}{T\alpha}\sum_{t=1}^{T-1}\Vert\nabla f({\mathbf{X}}_{t+1})\Vert_{2}\sqrt{B({\mathbf{Z}}_{t+1},{\mathbf{X}}_{t+1})},
\end{align*}
where the last equality follows since ${\mathbf{Z}}_1={\mathbf{X}}_1$.

\end{proof}

\begin{remark}Note that with the standard choice of step-size $\eta = 1/\beta$, whenever the sum $\sum_{t=1}^{T-1}\max\{B({\mathbf{X}}^*,{\mathbf{X}}_{t+1})-B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1}), \sqrt{B({\mathbf{Z}}_{t+1},{\mathbf{X}}_{t+1})}\}$ grows sublinearly in $T$, Theorem \ref{thm:converganceSmooth} indeed implies the convergence of Algorithm \ref{alg:LRMD} with deterministic gradients. In particular, when this sum is bounded by  a constant the rate is $O(1/T)$.
\end{remark}

\begin{theorem} \label{thm:StochasticConvergence}[Convergence with stochastic updates]
Let $\lbrace {\mathbf{X}}_t\rbrace_{t\geq 1}$ be a sequence of points in $\mathcal{S}_n\cap\mathbb{S}^n_{++}$.
Let $\{{\mathbf{Z}}_t\}_{t\geq 1}$ be a sequence such that ${\mathbf{Z}}_1 = {\mathbf{X}}_1$, and for all $t\geq 1$, ${\mathbf{Z}}_{t+1}$ is the accurate stochastic MEG update to ${\mathbf{X}}_t$ as defined in \eqref{eq:MDvonNeumannStochasticUpdate} with fixed step-size $\eta>0$. Then, after $T$ iterations, letting $\bar{{\mathbf{X}}} \sim \textrm{Uni}\{{\mathbf{X}}_1,\dots,{\mathbf{X}}_T\}$, it holds that
\begin{align*}
\mathbb{E}\left[f(\bar{{\mathbf{X}}})\right]-f({\mathbf{X}}^*) & \le \frac{B({\mathbf{X}}^*,{\mathbf{X}}_1)+\frac{G^2}{2}T\eta^2+\sum_{t=1}^T \mathbb{E}\left[B({\mathbf{X}}^*,{\mathbf{X}}_{t+1})-B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})\right]}{T\eta},
\end{align*}
where $G \geq \max_{t}\Vert{\widehat{\nabla}_t}\Vert_2$.
\end{theorem}

\begin{proof}

By the definition of ${\mathbf{Z}}_{t+1}$ in \eqref{eq:MDGeneralUpdate}, where the exact gradient is replaced with a stochastic gradient $\widehat{\nabla}_t$,
\[ {\mathbf{Z}}_{t+1}=\argmin_{{\mathbf{X}}\in\mathcal{S}_n}\{\langle\eta_t\widehat{\nabla}_t-\nabla\omega({\mathbf{X}}_t),{\mathbf{X}}\rangle+\omega({\mathbf{X}})\}. \]

Therefore, by the optimality condition for ${\mathbf{Z}}_{t+1}$, $\forall {\mathbf{X}}\in\mathcal{S}_n$,
\[ \langle \eta_t\widehat{\nabla}_t-\nabla\omega({\mathbf{X}}_t)+\nabla\omega({\mathbf{Z}}_{t+1}),{\mathbf{X}}-{\mathbf{Z}}_{t+1}\rangle\ge0. \]

Rearranging and using the three point lemma \eqref{lemma:threePointLemma}, we get that $\forall {\mathbf{X}}\in\mathcal{S}_n$,
\begin{align} \label{ineq:StochasticOptZ(t+1)}
\eta_t\langle\widehat{\nabla}_t,{\mathbf{Z}}_{t+1}-{\mathbf{X}}\rangle & \le \langle\nabla\omega({\mathbf{X}}_t)-\nabla\omega({\mathbf{Z}}_{t+1})),{\mathbf{Z}}_{t+1}-{\mathbf{X}}\rangle \nonumber
\\ & = B({\mathbf{X}},{\mathbf{X}}_t)-B({\mathbf{X}},{\mathbf{Z}}_{t+1})-B({\mathbf{Z}}_{t+1},{\mathbf{X}}_t).
\end{align}

Adding $\eta_t\langle \widehat{\nabla}_t,{\mathbf{X}}_t-{\mathbf{Z}}_{t+1}\rangle$ to both sides of \eqref{ineq:StochasticOptZ(t+1)} with ${\mathbf{X}}={\mathbf{X}}^*$, we obtain
\begin{align} \label{ineq:StochasticOptZ(t+1)_2}
\eta_t\langle\widehat{\nabla}_t,{\mathbf{X}}_{t}-{\mathbf{X}}^*\rangle & \le B({\mathbf{X}}^*,{\mathbf{X}}_t)-B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})-B({\mathbf{Z}}_{t+1},{\mathbf{X}}_t)+\eta_t\langle \widehat{\nabla}_t,{\mathbf{X}}_t-{\mathbf{Z}}_{t+1}\rangle.
\end{align}

Using H\"{o}lder's inequality, and \eqref{ineq:strongConvexityOfBregmanDistance} since $B(\cdot,\cdot)$ is $1$-strongly convex, we have that
\begin{align*}
-B({\mathbf{Z}}_{t+1},{\mathbf{X}}_t)+\eta_t\langle \widehat{\nabla}_t,{\mathbf{X}}_t-{\mathbf{Z}}_{t+1}\rangle & \le -\frac{1}{2}\Vert {\mathbf{Z}}_{t+1}-{\mathbf{X}}_t\Vert_*^2 + \eta_t\Vert\widehat{\nabla}_t\Vert_{2}\cdot\Vert{\mathbf{Z}}_{t+1}-{\mathbf{X}}_t\Vert_*
\\ & \le \max_{a\in\mathbb{R}}\left\{-\frac{1}{2}a^2+a\eta_t\Vert\widehat{\nabla}_t\Vert_{2}\right\}=\frac{1}{2}\eta_t^2\Vert\widehat{\nabla}_t\Vert_{2}^2.
\end{align*}

Plugging this into \eqref{ineq:StochasticOptZ(t+1)_2} we get
\begin{align*}
\eta_t\langle\widehat{\nabla}_t,{\mathbf{X}}_{t}-{\mathbf{X}}^*\rangle & \le B({\mathbf{X}}^*,{\mathbf{X}}_t)-B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})+\frac{1}{2}\eta_t^2\Vert\widehat{\nabla}_t\Vert_{2}^2
\\ & = B({\mathbf{X}}^*,{\mathbf{X}}_t)-B({\mathbf{X}}^*,{\mathbf{X}}_{t+1})+\frac{1}{2}\eta_t^2\Vert\widehat{\nabla}_t\Vert_{2}^2
\\ & \ \ \ \ +B({\mathbf{X}}^*,{\mathbf{X}}_{t+1})-B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1}).
\end{align*}

Summing over $t=1,...,T$ and dividing by $T$ we get
\begin{align*}
\frac{1}{T}\sum_{t=1}^T \eta_t \langle\widehat{\nabla}_t,{\mathbf{X}}_{t}-{\mathbf{X}}^*\rangle & \le \frac{1}{T}[B({\mathbf{X}}^*,{\mathbf{X}}_1)-B({\mathbf{X}}^*,{\mathbf{X}}_{T+1})]+\frac{1}{2T}\sum_{t=1}^T \eta_t^2\Vert\widehat{\nabla}_t\Vert_{2}^2
\\ & \ \ \ \ +\frac{1}{T}\sum_{t=1}^T \left[B({\mathbf{X}}^*,{\mathbf{X}}_{t+1})-B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})\right]
\\ & \le \frac{1}{T}B({\mathbf{X}}^*,{\mathbf{X}}_1)+\frac{1}{2T}\sum_{t=1}^T \eta_t^2\Vert\widehat{\nabla}_t\Vert_{2}^2\\ & \ \ \ \ +\frac{1}{T}\sum_{t=1}^T \left[B({\mathbf{X}}^*,{\mathbf{X}}_{t+1})-B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})\right].
\end{align*}

Taking expectation on both sides, we obtain
\begin{align*}
& \mathbb{E}\left[\frac{1}{T}\sum_{t=1}^T \eta_t \langle\widehat{\nabla}_t,{\mathbf{X}}_{t}-{\mathbf{X}}^*\rangle\right] \\ & \le \frac{1}{T}B({\mathbf{X}}^*,{\mathbf{X}}_1)+\frac{G^2}{2T}\sum_{t=1}^T \eta_t^2+\frac{1}{T}\sum_{t=1}^T \mathbb{E}\left[B({\mathbf{X}}^*,{\mathbf{X}}_{t+1})-B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})\right].
\end{align*}

It addition, using the law of total expectation, the gradient inequality, and taking $\eta_t=\eta$,
\begin{align*}
& \mathbb{E}\left[\frac{1}{T}\sum_{t=1}^T \eta \langle\widehat{\nabla}_t,{\mathbf{X}}_{t}-{\mathbf{X}}^*\rangle\right]  \\ & = \mathbb{E}\left[\frac{1}{T}\sum_{t=1}^T \eta \langle\mathbb{E}_{t}[\widehat{\nabla}_t\vert{\mathbf{X}}_t],{\mathbf{X}}_{t}-{\mathbf{X}}^*\rangle\right] = \mathbb{E}\left[\frac{1}{T}\sum_{t=1}^T \eta \langle\nabla f({\mathbf{X}}_t),{\mathbf{X}}_{t}-{\mathbf{X}}^*\rangle\right]
\\ & \ge \mathbb{E}\left[\frac{1}{T}\sum_{t=1}^T \eta (f({\mathbf{X}}_t)-f({\mathbf{X}}^*))\right] = \eta \left(\mathbb{E}\left[\frac{1}{T}\sum_{t=1}^T f({\mathbf{X}}_t)\right]-f({\mathbf{X}}^*)\right).
\end{align*}

Noticing that $\mathbb{E}\left[\frac{1}{T}\sum_{t=1}^T f({\mathbf{X}}_t)\right]=\mathbb{E}\left[f(\bar{{\mathbf{X}}})\right]$ and combining the last two inequalities we obtain
\begin{align*}
\mathbb{E}\left[f(\bar{{\mathbf{X}}})\right]-f({\mathbf{X}}^*) & \le \frac{B({\mathbf{X}}^*,{\mathbf{X}}_1)+\frac{G^2}{2}T\eta^2+\sum_{t=1}^T \mathbb{E}\left[B({\mathbf{X}}^*,{\mathbf{X}}_{t+1})-B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})\right]}{T\eta}.
\end{align*}

\end{proof}

\begin{remark}
Note that with a standard step-size of the form $\eta = C/\sqrt{T}$, for a suitable constant $C>0$, whenever the sum $\sum_{t=1}^{T}\mathbb{E}[B({\mathbf{X}}^*,{\mathbf{X}}_{t+1})-B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})]$ grows sublinearly in $\sqrt{T}$, Theorem \ref{thm:StochasticConvergence} indeed implies the convergence of Algorithm \ref{alg:LRMD} with stochastic gradients. In particular, when this sum is bounded by a constant the rate is $O(1/\sqrt{T})$.
\end{remark}

\subsubsection{Computing certificates for convergence}\label{sec:certificate}

As can be seen in Theorems \ref{thm:converganceSmooth}, \ref{thm:StochasticConvergence}, the convergence of our Algorithm \ref{alg:LRMD} naturally depends on the approximation errors $\left(B({\mathbf{X}}^*,{\mathbf{X}}_{t+1})-B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})\right)$ and $B({\mathbf{Z}}_{t+1},{\mathbf{X}}_{t+1})$. Thus, naturally, a significant portion of our analysis will be devoted to bounding these errors. Our first step is the following lemma.

\begin{lemma} \label{lemma:sufficientConvegenceCondition}
Let ${\mathbf{X}}\in\mathcal{S}_n$, ${\mathbf{X}}\succ 0$. Let ${\mathbf{Z}}_{+}$ be the accurate MEG update to ${\mathbf{X}}$ as defined in \eqref{eq:MDvonNeumannUpdate} or \eqref{eq:MDvonNeumannStochasticUpdate}, and ${\mathbf{X}}_{+}$ be the corresponding low-rank MEG update as defined in \eqref{eq:LowRankMDvonNeumannUpdate}.
Then, for any $\varepsilon\in(0,3/4]$,
\begin{align*}
\max\lbraceB({\mathbf{X}}^*,{\mathbf{X}}_{+})-B({\mathbf{X}}^*,{\mathbf{Z}}_{+}),B({\mathbf{Z}}_{+},{\mathbf{X}}_{+})\rbrace \le \max\left\lbrace 2\varepsilon, \log\left(\frac{(n-r)\lambda_{r+1}({\mathbf{Y}})}{\varepsilon \sum_{i=1}^n \lambda_i({\mathbf{Y}})}\right)\right\rbrace,
\end{align*}
where ${\mathbf{Y}} =\exp(\log ({\mathbf{X}})-\eta\nabla{}f({\mathbf{X}}))$ for deterministic gradients or ${\mathbf{Y}} =\exp(\log ({\mathbf{X}})-\eta\widehat{\nabla})$ for stochastic gradients with $\widehat{\nabla}$ being an unbiased estimator for $\nabla{}f({\mathbf{X}})$.
\end{lemma}

\begin{proof}

The von Neumann inequality claims that for any ${\mathbf{X}},{\mathbf{Y}}\in\mathbb{S}^n$, it holds that $\langle {\mathbf{X}},{\mathbf{Y}}\rangle \le \sum_{i=1}^n \lambda_i({\mathbf{X}})\lambda_i({\mathbf{Y}})$. Therefore,
\begin{align} \label{bound1ErrorTerm}
& B({\mathbf{X}}^*,{\mathbf{X}}_{+})-B({\mathbf{X}}^*,{\mathbf{Z}}_{+}) \nonumber \\ & = \textrm{Tr}({\mathbf{X}}^*\log({\mathbf{X}}^*)-{\mathbf{X}}^*\log({\mathbf{X}}_{+}))-\textrm{Tr}({\mathbf{X}}^*\log({\mathbf{X}}^*)-{\mathbf{X}}^*\log({\mathbf{Z}}_{+})) \nonumber
\\ & = \textrm{Tr}\left({\mathbf{X}}^*\left[\log({\mathbf{Z}}_{+})-\log({\mathbf{X}}_{+})\right]\right) \nonumber
\\ & \le \sum_{i=1}^n \lambda_i({\mathbf{X}}^*)\lambda_i(\log({\mathbf{Z}}_{+})-\log({\mathbf{X}}_{+}) \nonumber
\\ & \le \lambda_1(\log({\mathbf{Z}}_{+})-\log({\mathbf{X}}_{+})).
\end{align}

Using similar arguments, we get the bound
\begin{align} \label{bound2ErrorTerm}
B({\mathbf{Z}}_{+},{\mathbf{X}}_{+}) & = \textrm{Tr}({\mathbf{Z}}_{+}\log({\mathbf{Z}}_{+})-{\mathbf{Z}}_{+}\log({\mathbf{X}}_{+})) \nonumber
\\ & \le \sum_{i=1}^n \lambda_i({\mathbf{Z}}_{+})\lambda_i(\log({\mathbf{Z}}_{+})-\log({\mathbf{X}}_{+})
 \le \lambda_1(\log({\mathbf{Z}}_{+})-\log({\mathbf{X}}_{+})),
\end{align}
where the last inequality holds since ${\mathbf{Z}}_{t+1}\in\mathcal{S}_n$.

Combining \eqref{bound1ErrorTerm} and \eqref{bound2ErrorTerm}, we obtain the following bound:
\begin{align*}
\max\lbraceB({\mathbf{X}}^*,{\mathbf{X}}_{+})-B({\mathbf{X}}^*,{\mathbf{Z}}_{+}),B({\mathbf{Z}}_{+},{\mathbf{X}}_{+})\rbrace
 \le \lambda_1(\log({\mathbf{Z}}_{+})-\log({\mathbf{X}}_{+})).
\end{align*}

Let ${\mathbf{Y}}={\mathbf{V}}\Lambda{\mathbf{V}}^{\top}$ denote the eigen-decomposition of ${\mathbf{Y}}$. Therefore, its rank-$r$ approximation is ${\mathbf{Y}}^r={\mathbf{V}}^r\Lambda^r{{\mathbf{V}}^r}^{\top}$. ${\mathbf{Z}}_{+}$ and ${\mathbf{X}}_{+}$ have the same eigenvectors, and therefore $\log({\mathbf{Z}}_{+})-\log({\mathbf{X}}_{+})={\mathbf{V}}\left[\log\left(\frac{\Lambda}{b}\right)-\log\left(\begin{array}{cc}(1-\varepsilon)\frac{\Lambda^r}{a} & 0 \\ 0 & \frac{\varepsilon}{n-r}{\mathbf{I}}\end{array}\right)\right]{\mathbf{V}}^{\top}$.
Denote ${\mathbf{D}}=\log\left(\frac{\Lambda}{b}\right)-\log\left(\begin{array}{cc}(1-\varepsilon)\frac{\Lambda^r}{a} & 0 \\ 0 & \frac{\varepsilon}{n-r}{\mathbf{I}}\end{array}\right)$.
It holds that
\begin{align}
\forall j\in[r] &:\ \lambda_j({\mathbf{D}}) = \log\left(\frac{\frac{\lambda_j({\mathbf{Y}})}{b}}{(1-\varepsilon)\frac{\lambda_j({\mathbf{Y}})}{a}}\right); \label{StochasticBigEig}
\\ \forall j>r & :\ \lambda_j({\mathbf{D}}) = \log\left(\frac{\frac{\lambda_j({\mathbf{Y}})}{b}}{\frac{\varepsilon}{n-r}}\right). \label{StochasticLittleEig}
\end{align}
It is important to note that $\lambda_1({\mathbf{D}}),\ldots,\lambda_n({\mathbf{D}})$ are not necessarily ordered in a non-increasing order.

From \eqref{StochasticBigEig}, for all $j\le r$ it holds that
\begin{align*}
\lambda_j({\mathbf{D}}) & = \log\left(\frac{\lambda_j({\mathbf{Y}})\sum_{i=1}^{r}\lambda_i({\mathbf{Y}})}{(1-\varepsilon)\lambda_j({\mathbf{Y}})\sum_{i=1}^n\lambda_i({\mathbf{Y}})}\right)
\\ & = \log\left(\frac{\sum_{i=1}^{r}\lambda_i({\mathbf{Y}})}{\sum_{i=1}^n\lambda_i({\mathbf{Y}})}\right)+ \log\left(\frac{1}{1-\varepsilon}\right)
\le \log\left(\frac{1}{1-\varepsilon}\right) \underset{(a)}{\le} 2\varepsilon,
\end{align*}
where (a) holds for any $\varepsilon\le\frac{3}{4}$.

For $j\ge r+1$, note that using the definition of the eigenvalues of ${\mathbf{D}}$ in \eqref{StochasticLittleEig}, it holds that
\begin{align*}
\lambda_j({\mathbf{D}}) & = \log\left(\frac{(n-r)\lambda_j({\mathbf{Y}})}{\varepsilon \sum_{i=1}^n \lambda_i({\mathbf{Y}})}\right)
\le \log\left(\frac{(n-r)\lambda_{r+1}({\mathbf{Y}})}{\varepsilon \sum_{i=1}^n \lambda_i({\mathbf{Y}})}\right) = \lambda_{r+1}({\mathbf{D}}).
\end{align*}

The last two inequalities put together give us the desired result.

\end{proof}

Lemma \ref{lemma:sufficientConvegenceCondition} gives us a sufficient condition so that the errors due to the low rank MEG updates in Algorithm \ref{alg:LRMD} do not become too significant.
As seen, a sufficient condition to guarantee that $\max\lbraceB({\mathbf{X}}^*,{\mathbf{X}}_{+})-B({\mathbf{X}}^*,{\mathbf{Z}}_{+}),B({\mathbf{Z}}_{+},{\mathbf{X}}_{+})\rbrace\le 2\varepsilon$ is that
\begin{equation} \label{convergenceCondition}
\log\left(\frac{(n-r)\lambda_{r+1}({\mathbf{Y}})}{\varepsilon \sum_{i=1}^n \lambda_i({\mathbf{Y}})}\right)\le 2\varepsilon.
\end{equation}

Unfortunately, the condition in \eqref{convergenceCondition} is useless as it involves computing the full matrix of ${\mathbf{Y}}$ in order to calculate the scalar $b=\sum_{i=1}^n {\lambda_i({\mathbf{Y}})}$, which is exactly what we are striving to avoid. Instead, $b$ can be approximated with some $b_k:=\sum_{i=1}^k {\lambda_i({\mathbf{Y}})}$ for some $k\le n$.
Thus, instead of the full eigen-decomposition of ${\mathbf{Y}}$, in order to check if  our Algorithm \ref{alg:LRMD} is converging with a tolerable error, we can check the weaker condition of $\log\left(\frac{(n-r)\lambda_{r+1}({\mathbf{Y}})}{\varepsilon b_k}\right)\le 2\varepsilon$, which is likewise a sufficient condition since $b_k \leq b$, and can be computed much more efficiently for $k<<n$. In particular, computing $b_{r+1}$, where $r$ is the SVD rank parameter in Algorithm \ref{alg:LRMD} requires to increase the rank of the SVD computations in Algorithm \ref{alg:LRMD} by only one\footnote{Since $\lambda_{r+1}({\mathbf{Y}})$ is also required in \eqref{convergenceCondition} there is no point in computing $b_k$ for $k< r+1$ which will result in a worse approximation of the scalar $b$.} . In Section \ref{section:experiments} we present numerical evidence which demonstrate that in practice the accurate certificate computed using $b=b_n$ does not give any additional benefit over the use of $b_{r+1}$. Moreover, in the experiments we demonstrate that this computationally-cheap certificate indeed holds true from the very early stages of the run (already from the first iteration in most cases), and hence, may indeed be appealing for practical uses.

\section{Provable Bounds on Low-Rank Approximations under Strict Complementarity and Warm-Start Initialization}\label{sec:main_tools}
In this section we lay down two of the technical foundations towards obtaining provable convergence results for our Algorithm \ref{alg:LRMD}  from a ``warm-start" initialization under the strict complementarity condition (Assumption \ref{ass:strictcomp}).

In the following lemma we present our main technical result: for any point sufficiently close to an optimal solution ${\mathbf{X}}^*$ for which Assumption \ref{ass:strictcomp} holds, the error in using low rank MEG updates is indeed guaranteed to be properly bounded.

\begin{lemma} \label{lemma:StochasticRadius}
Let ${\mathbf{X}}^*\in\mathcal{X}^*$ be an optimal solution for which Assumption \ref{ass:strictcomp} holds with $\delta >0$. Let $\textrm{rank}({\mathbf{X}}^*):=r^*$ and let $n\not=r\ge r^*$ be the SVD rank parameter. Let ${\mathbf{X}}\in\mathcal{S}_n,{\mathbf{X}}\succ 0$ be a matrix such that $\lambda_{r^*+1}({\mathbf{X}})\le\frac{\varepsilon_{-}}{n-r}$ for some $\varepsilon_{-}\in[0,1]$ and let $\widetilde{\nabla}\in\mathbb{S}^n$ be such that $\Vert\widetilde{\nabla}-\nabla f({\mathbf{X}})\Vert\le\xi$ for some $\xi\ge0$ (for deterministic updates $\xi=0$). Denote $G=\sup_{{\mathbf{X}}\in\mathcal{S}_n}\Vert \nabla f({\mathbf{X}})\Vert_2$. \ If
\begin{equation} \label{eq:StochasticRadiusOfconvergence}
\sqrt{B({\mathbf{X}}^*,{\mathbf{X}})} \le \frac{1}{\sqrt{2}}\left[2\beta+\left(1+\frac{2\sqrt{2r^*}}{\lambda_{r^*}({\mathbf{X}}^*)}\right)G\right]^{-1}\left(\delta-\frac{1}{\eta}\log\left(\frac{\varepsilon_{-}}{\varepsilon}\right)+\frac{2\varepsilon}{\eta}-2\xi\right),
\end{equation}
then, for any $\eta>0$ the matrix ${\mathbf{Y}} = \exp\left({\log({\mathbf{X}})-\eta\widetilde{\nabla}}\right)$ and scalar $b=\sum_{i=1}^n\lambda_i({\mathbf{Y}})$ satisfy:
\begin{align}
\log\left(\frac{(n-r)\lambda_{r+1}({\mathbf{Y}})}{\varepsilon b}\right) \le 2\varepsilon.
\end{align}
As a result, for ${\mathbf{Z}}_+ := {\mathbf{Y}}/b$ and ${\mathbf{X}}_+$ as defined in Eq. \eqref{eq:LowRankMDvonNeumannUpdate}, it holds that
\begin{align*}
\max\lbraceB({\mathbf{X}}^*,{\mathbf{X}}_{+})-B({\mathbf{X}}^*,{\mathbf{Z}}_{+}),B({\mathbf{Z}}_{+},{\mathbf{X}}_{+})\rbrace \le 2\varepsilon.
\end{align*}
\end{lemma}

\begin{proof}

It holds that
\begin{align} \label{ineq:StochasticLittleEigBound}
\log\left(\frac{(n-r)\lambda_{r^*+1}({\mathbf{Y}})}{\varepsilon \sum_{i=1}^n \lambda_i({\mathbf{Y}})}\right) & = \log\left(\frac{n-r}{\varepsilon}\right) + \log(\lambda_{r+1}({\mathbf{Y}}))-\log\left( \sum_{i=1}^n \lambda_i({\mathbf{Y}})\right) \nonumber
\\ & \le \log\left(\frac{n-r}{\varepsilon}\right) + \log(\lambda_{r^*+1}({\mathbf{Y}}))-\log\left( \sum_{i=1}^n \lambda_i({\mathbf{Y}})\right).
\end{align}
We will now separately bound the last two terms.

\begin{align}  \label{ineq:StochasticInProofLambdaRPlus1}
\log(\lambda_{r^*+1}({\mathbf{Y}})) & = \lambda_{r^*+1}(\log({\mathbf{Y}})) \underset{(a)}{\le} \lambda_{r^*+1}(\log({\mathbf{X}})-\eta\nabla f({\mathbf{X}}^*)) + \eta\lambda_1(\nabla f({\mathbf{X}}^*)-\widetilde{\nabla}) \nonumber
\\ & \le \sum_{i=1}^{r^*+1}\lambda_{i}(\log({\mathbf{X}})-\eta\nabla f({\mathbf{X}}^*))-\sum_{i=1}^{r^*}\lambda_{i}(\log({\mathbf{X}})-\eta\nabla f({\mathbf{X}}^*))\nonumber \\ & \ \ \  +\eta\Vert\widetilde{\nabla}-\nabla f({\mathbf{X}}^*)\Vert_2 \nonumber
\\ & \underset{(b)}{\le} \sum_{i=1}^{r^*+1}\lambda_{i}(\log({\mathbf{X}}))-\eta\sum_{i=1}^{r^*+1}\lambda_{n-i+1}(\nabla f({\mathbf{X}}^*))-\sum_{i=1}^{r^*}\lambda_{i}(\log({\mathbf{X}})-\eta_t\nabla f({\mathbf{X}}^*)) \nonumber \\ & \ \ \ +\eta\Vert\widetilde{\nabla}-\nabla f({\mathbf{X}}^*)\Vert_2 \nonumber
\\ & \underset{(c)}{=} \sum_{i=1}^{r^*}\lambda_{i}(\log({\mathbf{X}}))+\log\left(\lambda_{r^*+1}({\mathbf{X}})\right)-\eta r^*\lambda_n(\nabla f({\mathbf{X}}^*))-\eta \lambda_{n-r^*}(\nabla f({\mathbf{X}}^*))\nonumber \\ & \ \ \ -\sum_{i=1}^{r^*}\lambda_{i}(\log({\mathbf{X}})-\eta\nabla f({\mathbf{X}}^*))+\eta\Vert\widetilde{\nabla}-\nabla f({\mathbf{X}}^*)\Vert_2,
\end{align}
where (a) follows from Weyl's inequality, (b) follows from Ky Fan's inequality for eigenvalues and (c) follows from Lemma \ref{lemma:opt-grad}.

Since $f$ is $\beta$-smooth,
\begin{align} \label{ineq:StochasticNorm-nablat-gradopt}
\Vert\widetilde{\nabla}-\nabla f({\mathbf{X}}^*)\Vert_2 \le \Vert\widetilde{\nabla}-\nabla f({\mathbf{X}})\Vert_2+\Vert\nabla f({\mathbf{X}})-\nabla f({\mathbf{X}}^*)\Vert_2 \le \xi+\beta\Vert{\mathbf{X}}-{\mathbf{X}}^*\Vert_*.
\end{align}

Denote $\log({\mathbf{X}})={\mathbf{W}}{\mathbf{D}}{\mathbf{W}}^{\top}$ and ${\mathbf{W}}_{r^*}$ to be the matrix with the $r^*$ first columns of ${\mathbf{W}}$.
\begin{align} \label{ineq:inProofSumReigs}
-\sum_{i=1}^{r^*}\lambda_{i}(\log({\mathbf{X}})-\eta\nabla f({\mathbf{X}}^*)) & \le -{\mathbf{W}}_{r^*}{\mathbf{W}}_{r^*}^{\top}\bullet(\log({\mathbf{X}})-\eta\nabla f({\mathbf{X}}^*)) \nonumber
\\ & = -\sum_{i=1}^{r^*}\lambda_i(\log({\mathbf{X}}))+\eta{\mathbf{W}}_{r^*}{\mathbf{W}}_{r^*}^{\top}\bullet\nabla f({\mathbf{X}}^*).
\end{align}

Let ${\mathbf{X}}^*={\mathbf{V}}^*\Lambda^*{{\mathbf{V}}^*}^{\top}$ denote the eigen-decomposition of ${\mathbf{X}}^*$. From Lemma \ref{lemma:opt-grad} this implies that ${\mathbf{V}}^*{{\mathbf{V}}^*}^{\top}\bullet\nabla f({\mathbf{X}}^*)=r^*\lambda_{n}(\nabla f({\mathbf{X}}^*))$. Therefore,
\begin{align} \label{ineq:inProofInnerProd}
{\mathbf{W}}_{r^*}{\mathbf{W}}_{r^*}^{\top}\bullet\nabla f({\mathbf{X}}^*) & = {\mathbf{V}}^*{{\mathbf{V}}^*}^{\top}\bullet\nabla f({\mathbf{X}}^*)+{\mathbf{W}}_{r^*}{\mathbf{W}}_{r^*}^{\top}\bullet\nabla f({\mathbf{X}}^*)-{\mathbf{V}}^*{{\mathbf{V}}^*}^{\top}\bullet\nabla f({\mathbf{X}}^*) \nonumber
\\ & \le {\mathbf{V}}^*{{\mathbf{V}}^*}^{\top}\bullet\nabla f({\mathbf{X}}^*)+\Vert{\mathbf{W}}_{r^*}{\mathbf{W}}_{r^*}^{\top}-{\mathbf{V}}^*{{\mathbf{V}}^*}^{\top}\Vert_{*}\Vert\nabla f({\mathbf{X}}^*)\Vert_{2} \nonumber
\\ & = r^*\lambda_{n}(\nabla f({\mathbf{X}}^*))+\Vert{\mathbf{W}}_{r^*}{\mathbf{W}}_{r^*}^{\top}-{\mathbf{V}}^*{{\mathbf{V}}^*}^{\top}\Vert_{*}\Vert\nabla f({\mathbf{X}}^*)\Vert_{2} \nonumber
\\ & \le r^*\lambda_{n}(\nabla f({\mathbf{X}}^*))+\sqrt{2r^*}\Vert{\mathbf{W}}_{r^*}{\mathbf{W}}_{r^*}^{\top}-{\mathbf{V}}^*{{\mathbf{V}}^*}^{\top}\Vert_{F}\Vert\nabla f({\mathbf{X}}^*)\Vert_{2} \nonumber
\\ & \underset{(a)}{\le} r^*\lambda_{n}(\nabla f({\mathbf{X}}^*))+\frac{2\sqrt{2r^*}\Vert{\mathbf{X}}-{\mathbf{X}}^*\Vert_{F}}{\lambda_{r^*}({\mathbf{X}}^*)}\Vert\nabla f({\mathbf{X}}^*)\Vert_{2} \nonumber
\\ & \le r^*\lambda_{n}(\nabla f({\mathbf{X}}^*))+\frac{2\sqrt{2r^*}\Vert{\mathbf{X}}-{\mathbf{X}}^*\Vert_{*}}{\lambda_{r^*}({\mathbf{X}}^*)}\Vert\nabla f({\mathbf{X}}^*)\Vert_{2}.
\end{align}
Here, (a) follows using the Davis-Kahan $\sin(\theta)$ theorem (see for instance Theorem 2 in \cite{Davis-Kahan}).

Plugging \eqref{ineq:StochasticNorm-nablat-gradopt} and the bounds \eqref{ineq:inProofSumReigs} and \eqref{ineq:inProofInnerProd} into \eqref{ineq:StochasticInProofLambdaRPlus1}, we get
\begin{align}
\log(\lambda_{r^*+1}({\mathbf{Y}})) & \le \log\left(\lambda_{r^*+1}({\mathbf{X}})\right)-\eta \lambda_{n-r^*}(\nabla f({\mathbf{X}}^*))+\eta\frac{2\sqrt{2r^*}\Vert{\mathbf{X}}-{\mathbf{X}}^*\Vert_{*}}{\lambda_{r^*}({\mathbf{X}}^*)}\Vert\nabla f({\mathbf{X}}^*)\Vert_{2} \nonumber
\\ & \ \ \ +\eta\xi+\eta\beta\Vert{\mathbf{X}}-{\mathbf{X}}^*\Vert_* \label{ineq:StochasticInProofSecondTermWithoutAssumption}
\\ & \le \log\left(\frac{\varepsilon_{-}}{n-r}\right)-\eta \lambda_{n-r^*}(\nabla f({\mathbf{X}}^*))+\eta\frac{2\sqrt{2r^*}\Vert{\mathbf{X}}-{\mathbf{X}}^*\Vert_{*}}{\lambda_{r^*}({\mathbf{X}}^*)}\Vert\nabla f({\mathbf{X}}^*)\Vert_{2} \nonumber
\\ & \ \ \ +\eta\xi+\eta\beta\Vert{\mathbf{X}}-{\mathbf{X}}^*\Vert_*. \label{ineq:StochasticInProofSecondTerm}
\end{align}
where the second inequality holds from the assumption that $\lambda_{r^*+1}({\mathbf{X}})\le\frac{\varepsilon_{-}}{n-r}$.

For the last term of \eqref{ineq:StochasticLittleEigBound},
\begin{align} \label{ineq:StochasticInProofThirdTerm}
-\log\left(\sum_{i=1}^n\lambda_i({\mathbf{Y}})\right) & = -\log\Big(\sum_{i=1}^n\lambda_i\left(\exp\left(\log({\mathbf{X}})-\eta_t\widetilde{\nabla}\right)\right)\Big) \nonumber
\\ & \underset{(a)}{\le} -\sum_{i=1}^n\lambda_i({\mathbf{X}})\cdot\log\left(\frac{\lambda_i\left(\exp\left(\log({\mathbf{X}})-\eta\widetilde{\nabla}\right)\right)}{\lambda_i({\mathbf{X}})}\right) \nonumber
\\ & = -\sum_{i=1}^n\lambda_i({\mathbf{X}})\cdot\left[\log\left(\lambda_i\left(\exp\left(\log({\mathbf{X}})-\eta\widetilde{\nabla}\right)\right)\right)-\log(\lambda_i({\mathbf{X}})\right] \nonumber
\\ & = -\sum_{i=1}^n\lambda_i({\mathbf{X}})\cdot\left[\lambda_i\left(\log({\mathbf{X}})-\eta\widetilde{\nabla}\right)-\log(\lambda_i({\mathbf{X}})\right] \nonumber
\\ & \underset{(b)}{\le} -\langle{\mathbf{X}},\log({\mathbf{X}})-\eta\widetilde{\nabla}\rangle +\sum_{i=1}^n\lambda_i({\mathbf{X}})\lambda_i(\log({\mathbf{X}})) \nonumber
\\ & \underset{(c)}{=} -\langle{\mathbf{X}},\log({\mathbf{X}})-\eta\widetilde{\nabla}\rangle +\sum_{i=1}^n\lambda_i({\mathbf{X}}\log({\mathbf{X}})) \nonumber
\\ & = -\textrm{Tr}({\mathbf{X}}\log({\mathbf{X}}))+\textrm{Tr}({\mathbf{X}}\log({\mathbf{X}}))+\eta\langle{\mathbf{X}},\widetilde{\nabla}\rangle \nonumber
\\ & = \eta\langle{\mathbf{X}}^*,\nabla f({\mathbf{X}}^*)\rangle + \eta\left(\langle{\mathbf{X}}-{\mathbf{X}}^*,\nabla f({\mathbf{X}}^*)\rangle+\langle{\mathbf{X}},\widetilde{\nabla}-\nabla f({\mathbf{X}}^*)\rangle\right) \nonumber
\\ & \underset{(d)}{=} \eta\lambda_{n}(\nabla f({\mathbf{X}}^*)) + \eta\left(\langle{\mathbf{X}}-{\mathbf{X}}^*,\nabla f({\mathbf{X}}^*)\rangle+\langle{\mathbf{X}},\widetilde{\nabla}-\nabla f({\mathbf{X}}^*)\rangle\right) \nonumber
\\ & \le \eta\lambda_{n}(\nabla f({\mathbf{X}}^*)) + \eta\Vert{\mathbf{X}}-{\mathbf{X}}^*\Vert_*\Vert\nabla f({\mathbf{X}}^*)\Vert_2+\eta\Vert{\mathbf{X}}\Vert_*\Vert\widetilde{\nabla}-\nabla f({\mathbf{X}}^*)\Vert_2 \nonumber
\\ & \underset{(e)}{\le} \eta\lambda_{n}(\nabla f({\mathbf{X}}^*)) + \eta\Vert{\mathbf{X}}-{\mathbf{X}}^*\Vert_*\Vert\nabla f({\mathbf{X}}^*)\Vert_2+\eta(\xi+\beta\Vert{\mathbf{X}}-{\mathbf{X}}^*\Vert_*)
\end{align}
where (a) follows since $-\log(\cdot)$ is convex, (b) follows from the von Neumann inequality, (c) follows since ${\mathbf{X}}$ and $\log({\mathbf{X}})$ have the same eigen-vectors, (d) follows from Lemma \ref{lemma:opt-grad}, and (e) follows from \eqref{ineq:StochasticNorm-nablat-gradopt} and since $\Vert{\mathbf{X}}\Vert_*=1$.

Plugging \eqref{ineq:StochasticInProofSecondTerm} and \eqref{ineq:StochasticInProofThirdTerm} into \eqref{ineq:StochasticLittleEigBound}, we get
\begin{align*}
& \log\left(\frac{(n-r)\lambda_{r+1}({\mathbf{Y}})}{\varepsilon \sum_{i=1}^n \lambda_i({\mathbf{Y}})}\right)
\\ & \le -\eta\delta +\log\left(\frac{\varepsilon_{-}}{\varepsilon}\right)+ 2\eta\xi + \eta\left(2\beta+\left(1+\frac{2\sqrt{2r^*}}{\lambda_{r^*}({\mathbf{X}}^*)}\right)\Vert\nabla f({\mathbf{X}}^*)\Vert_{2}\right)\Vert{\mathbf{X}}-{\mathbf{X}}^*\Vert_*
\\ & \le -\eta\delta +\log\left(\frac{\varepsilon_{-}}{\varepsilon}\right)+ 2\eta\xi + \sqrt{2}\eta\left(2\beta+\left(1+\frac{2\sqrt{2r^*}}{\lambda_{r^*}({\mathbf{X}}^*)}\right)G\right)\sqrt{B({\mathbf{X}}^*,{\mathbf{X}})},
\end{align*}
where the last inequality follows from \eqref{ineq:strongConvexityOfBregmanDistance} due to the $1$-strong convexity of $B(\cdot,\cdot)$ and the definition of $G$.

Therefore, if ${\mathbf{X}}$ is a matrix such that
\begin{equation*}
\sqrt{B({\mathbf{X}}^*,{\mathbf{X}})} \le \frac{1}{\sqrt{2}}\left[2\beta+\left(1+\frac{2\sqrt{2r^*}}{\lambda_{r^*}({\mathbf{X}}^*)}\right)G\right]^{-1}\left(\delta-\frac{1}{\eta}\log\left(\frac{\varepsilon_{-}}{\varepsilon}\right)+\frac{2\varepsilon}{\eta}-2\xi\right),
\end{equation*}
then
\begin{align*}
\log\left(\frac{(n-r)\lambda_{r+1}({\mathbf{Y}})}{\varepsilon \sum_{i=1}^n \lambda_i({\mathbf{Y}})}\right)\le 2\varepsilon.
\end{align*}

\end{proof}

Since our provable convergence results only apply from a ``warm-start" initialization, we bring the following lemma which will be used to give precise conditions on the initialization parameters $({\mathbf{X}}_0,\varepsilon_0)$ in Algorithm \ref{alg:LRMD}, which are used to produce the first iterate ${\mathbf{X}}_1$, so that it is close enough to an optimal solution of interest.

\begin{lemma}[Warm-start Initialization]\label{lem:warmstart}
Let ${\mathbf{X}}^*\in\mathcal{X}^*$ be an optimal solution to Problem \eqref{eq:Model}. Let $\textrm{rank}({\mathbf{X}}^*):=r^*$, and let ${\mathbf{X}}\in\mathcal{S}_n$ be such that $\textrm{rank}({\mathbf{X}}) \geq \textrm{rank}({\mathbf{X}}^*)$. Let $\varepsilon\in(0,3/4]$ and fix some $R>0$. Suppose
\begin{align}\label{eq:init:highRank}
\textrm{Tr}({\mathbf{X}}^*\log{\mathbf{X}}^* - {\mathbf{X}}^*{\mathbf{V}}_{{\mathbf{X}}}\log(\Lambda_{{\mathbf{X}}}){\mathbf{V}}_{{\mathbf{X}}}^{\top}) +  \lambda_1({\mathbf{X}}^*)\log\left({\frac{n}{\varepsilon}}\right)(r^*-\Vert{{\mathbf{V}}_{{\mathbf{X}}^*}^{\top}{\mathbf{V}}_{{\mathbf{X}}}}\Vert_F^2)  + 4\varepsilon \leq R^2
\end{align}

holds, where  ${\mathbf{V}}_{{\mathbf{X}}}\Lambda_{{\mathbf{X}}}{\mathbf{V}}_{{\mathbf{X}}}^{\top}$ denotes the compact-form eigen-decomposition of ${\mathbf{X}}$ (i.e., $\Lambda_{{\mathbf{X}}}\succ\mathbf{0}$), and similarly ${\mathbf{V}}_{{\mathbf{X}}^*}\Lambda_{{\mathbf{X}}^*}{\mathbf{V}}_{{\mathbf{X}}^*}^{\top}$ denotes the compact-form eigen-decomposition of ${\mathbf{X}}^*$. Then, the matrix ${\mathbf{W}} = (1-\varepsilon){\mathbf{X}} + \frac{\varepsilon}{n }{\mathbf{I}}$ satisfies $B({\mathbf{X}}^*,{\mathbf{W}}) \leq R^2$.

Moreover, if $\textrm{rank}({\mathbf{X}}) = \textrm{rank}({\mathbf{X}}^*)$ Then, the condition in \eqref{eq:init:highRank} can be replaced with
\begin{align}\label{eq:init:equalRank}
\textrm{Tr}({\mathbf{X}}^*\log{\mathbf{X}}^* - {\mathbf{X}}^*{\mathbf{V}}_{{\mathbf{X}}}\log(\Lambda_{{\mathbf{X}}}){\mathbf{V}}_{{\mathbf{X}}}^{\top}) +  \frac{2\lambda_1({\mathbf{X}}^*)}{\lambda_{r^*}({\mathbf{X}}^*)^2}\log\left({\frac{n}{\varepsilon}}\right)\Vert{{\mathbf{X}}^*-{\mathbf{X}}}\Vert_F^2  + 4\varepsilon \leq R^2.
\end{align}
\end{lemma}
\begin{proof}
From simple calculations we have
\begin{align*}
-\log({\mathbf{W}}) &= -\log\left({(1-\varepsilon){\mathbf{X}} + \frac{\varepsilon}{n}{\mathbf{I}}}\right)  \\
&= -{\mathbf{V}}_{{\mathbf{X}}}\log\left({(1-\varepsilon)\Lambda_{{\mathbf{X}}}+\frac{\varepsilon}{n}{\mathbf{I}}}\right){\mathbf{V}}_{{\mathbf{X}}}^{\top} - \log\left({\frac{\varepsilon}{n}}\right)({\mathbf{I}}-{\mathbf{V}}_{{\mathbf{X}}}{\mathbf{V}}_{{\mathbf{X}}}^{\top}) \\
&\preceq -{\mathbf{V}}_{{\mathbf{X}}}\log\left({(1-\varepsilon)\Lambda_{{\mathbf{X}}}}\right){\mathbf{V}}_{{\mathbf{X}}}^{\top} - \log\left({\frac{\varepsilon}{n}}\right)({\mathbf{I}}-{\mathbf{V}}_{{\mathbf{X}}}{\mathbf{V}}_{{\mathbf{X}}}^{\top})\\
&=-{\mathbf{V}}_{{\mathbf{X}}}\left({\log(1-\varepsilon){\mathbf{I}}+\log(\Lambda_{{\mathbf{X}}})}\right){\mathbf{V}}_{{\mathbf{X}}}^{\top} + \log\left({\frac{n}{\varepsilon}}\right)({\mathbf{I}}-{\mathbf{V}}_{{\mathbf{X}}}{\mathbf{V}}_{{\mathbf{X}}}^{\top}) \\
&\preceq -\log(1-\varepsilon){\mathbf{I}} - {\mathbf{V}}_{{\mathbf{X}}}\log(\Lambda_{{\mathbf{X}}}){\mathbf{V}}_{{\mathbf{X}}}^{\top} + \log\left({\frac{n}{\varepsilon}}\right)({\mathbf{I}}-{\mathbf{V}}_{{\mathbf{X}}}{\mathbf{V}}_{{\mathbf{X}}}^{\top}).
\end{align*}

Thus,
\begin{align*}
B({\mathbf{X}}^*,{\mathbf{W}}) &= \textrm{Tr}({\mathbf{X}}^*\log({\mathbf{X}}^*) - {\mathbf{X}}^*\log({\mathbf{W}})) \\
&\leq \textrm{Tr}\left({{\mathbf{X}}^*\log({\mathbf{X}}^*)  -{\mathbf{X}}^*{\mathbf{V}}_{{\mathbf{X}}}\log(\Lambda_{{\mathbf{X}}}){\mathbf{V}}_{{\mathbf{X}}}^{\top}}\right) \\
&~~~ - \textrm{Tr}\left({{\mathbf{X}}^*\left({\log(1-\varepsilon){\mathbf{I}}-\log\left({\frac{n}{\varepsilon}}\right)({\mathbf{I}}-{\mathbf{V}}_{{\mathbf{X}}}{\mathbf{V}}_{{\mathbf{X}}}^{\top})}\right)}\right)\\
&=\textrm{Tr}({\mathbf{X}}^*\log{\mathbf{X}}^* - {\mathbf{X}}^*{\mathbf{V}}_{{\mathbf{X}}}\log(\Lambda_{{\mathbf{X}}}){\mathbf{V}}_{{\mathbf{X}}}^{\top})\\
&~~~+\log\left({\frac{n}{\varepsilon}}\right)\textrm{Tr}\left({{\mathbf{X}}^*({\mathbf{I}}-{\mathbf{V}}_{{\mathbf{X}}}{\mathbf{V}}_{{\mathbf{X}}}^{\top})}\right)-\log(1-\varepsilon).
\end{align*}

We continue to bound the last two terms.
\begin{align*}
\textrm{Tr}\left({{\mathbf{X}}^*({\mathbf{I}}-{\mathbf{V}}_{{\mathbf{X}}}{\mathbf{V}}_{{\mathbf{X}}}^{\top})}\right) &= \textrm{Tr}\left({{\mathbf{V}}_{{\mathbf{X}}^*}\Lambda_{{\mathbf{X}}^*}{\mathbf{V}}_{{\mathbf{X}}^*}^{\top}({\mathbf{I}}-{\mathbf{V}}_{{\mathbf{X}}}{\mathbf{V}}_{{\mathbf{X}}}^{\top})}\right)\\
&= \textrm{Tr}\left({\Lambda_{{\mathbf{X}}^*}({\mathbf{I}}-{\mathbf{V}}_{{\mathbf{X}}^*}^{\top}{\mathbf{V}}_{{\mathbf{X}}}{\mathbf{V}}_{{\mathbf{X}}}^{\top}{\mathbf{V}}_{{\mathbf{X}}^*})}\right)\\
&\leq \lambda_1({\mathbf{X}}^*)\textrm{Tr}\left({{\mathbf{I}}-{\mathbf{V}}_{{\mathbf{X}}^*}^{\top}{\mathbf{V}}_{{\mathbf{X}}}{\mathbf{V}}_{{\mathbf{X}}}^{\top}{\mathbf{V}}_{{\mathbf{X}}^*}}\right)\\
&=\lambda_1({\mathbf{X}}^*)\left({r^*-\textrm{Tr}\left({{\mathbf{V}}_{{\mathbf{X}}^*}^{\top}{\mathbf{V}}_{{\mathbf{X}}}{\mathbf{V}}_{{\mathbf{X}}}^{\top}{\mathbf{V}}_{{\mathbf{X}}^*}}\right)}\right)\\
&= \lambda_1({\mathbf{X}}^*)\left({r^*-\Vert{{\mathbf{V}}_{{\mathbf{X}}^*}^{\top}{\mathbf{V}}_{{\mathbf{X}}}}\Vert_F^2}\right).
\end{align*}

Also, for all $\varepsilon\in(0,3/4]$ we have $- \log(1-\varepsilon) \leq \frac{\varepsilon}{1-\varepsilon} \leq 4\varepsilon$.

Thus, we conclude that
\begin{align*}
B({\mathbf{X}}^*,{\mathbf{W}}) &\leq\textrm{Tr}({\mathbf{X}}^*\log{\mathbf{X}}^* - {\mathbf{X}}^*{\mathbf{V}}_{{\mathbf{X}}}\log(\Lambda_{{\mathbf{X}}}){\mathbf{V}}_{{\mathbf{X}}}^{\top}) \\
&~~~+   \log\left({\frac{n}{\varepsilon}}\right)\lambda_1({\mathbf{X}}^*)\left({r^*-\Vert{{\mathbf{V}}_{{\mathbf{X}}^*}^{\top}{\mathbf{V}}_{{\mathbf{X}}}}\Vert_F^2}\right)  + 4\varepsilon.
\end{align*}

In case $\textrm{rank}({\mathbf{X}})=\textrm{rank}({\mathbf{X}}^*)$ we can replace the upper-bound on $\textrm{Tr}\left({{\mathbf{X}}^*({\mathbf{I}}-{\mathbf{V}}_{{\mathbf{X}}}{\mathbf{V}}_{{\mathbf{X}}}^{\top})}\right)$ with the following:
\begin{align*}
\textrm{Tr}\left({{\mathbf{X}}^*({\mathbf{I}}-{\mathbf{V}}_{{\mathbf{X}}}{\mathbf{V}}_{{\mathbf{X}}}^{\top})}\right) &= \textrm{Tr}\left({{\mathbf{V}}_{{\mathbf{X}}^*}\Lambda_{{\mathbf{X}}^*}{\mathbf{V}}_{{\mathbf{X}}^*}^{\top}({\mathbf{I}}-{\mathbf{V}}_{{\mathbf{X}}}{\mathbf{V}}_{{\mathbf{X}}}^{\top})}\right)\\
&= \textrm{Tr}\left({\Lambda_{{\mathbf{X}}^*}({\mathbf{I}}-{\mathbf{V}}_{{\mathbf{X}}^*}^{\top}{\mathbf{V}}_{{\mathbf{X}}}{\mathbf{V}}_{{\mathbf{X}}}^{\top}{\mathbf{V}}_{{\mathbf{X}}^*})}\right)\\
&\leq \lambda_1({\mathbf{X}}^*)\textrm{Tr}\left({{\mathbf{I}}-{\mathbf{V}}_{{\mathbf{X}}^*}^{\top}{\mathbf{V}}_{{\mathbf{X}}}{\mathbf{V}}_{{\mathbf{X}}}^{\top}{\mathbf{V}}_{{\mathbf{X}}^*}}\right)\\
&= \lambda_1({\mathbf{X}}^*)\left({r^*-\textrm{Tr}\left({{\mathbf{V}}_{{\mathbf{X}}^*}{\mathbf{V}}_{{\mathbf{X}}^*}^{\top}{\mathbf{V}}_{{\mathbf{X}}}{\mathbf{V}}_{{\mathbf{X}}}^{\top}}\right)}\right)\\
&= \frac{\lambda_1({\mathbf{X}}^*)}{2}\Vert{{\mathbf{V}}_{{\mathbf{X}}^*}{\mathbf{V}}_{{\mathbf{X}}^*}^{\top}-{\mathbf{V}}_{{\mathbf{X}}}{\mathbf{V}}_{{\mathbf{X}}}^{\top}}\Vert_F^2\\
&\leq \frac{2\lambda_1({\mathbf{X}}^*)}{\lambda_{r^*}({\mathbf{X}}^*)^2}\Vert{{\mathbf{X}}^*-{\mathbf{X}}}\Vert_F^2,
\end{align*}
where the last inequality follows from the Davis-Kahan $\sin\theta$ theorem (see for instance Theorem 2 in \cite{Davis-Kahan}).

\end{proof}

\section{Provable Local Convergence with Deterministic Updates under Strict Complementarity} \label{section:smoothCase}

In this section we prove a convergence rate of $O(1/t)$ for our  Algorithm \ref{alg:LRMD} from ``warm-start" initialization when exact gradients are used and the strict complementarity condition holds.

In the following lemma we show that there is a decrease in the Bregman distance between the iterates of the algorithm and the optimal solution up to a small error. This step is important since as discussed in Section \ref{section:localConvergence} we can only provably control the errors due to the low rank updates in a certain ball around an optimal solution, and hence, we must show that when the algorithm is initialized inside this ball, the iterates do not leave it.

\begin{lemma} \label{lemma:bregmanDecrease}
Let ${\mathbf{X}}^*\in\mathcal{X}^*$ be an optimal solution for which Assumption \ref{ass:strictcomp} holds with $\delta >0$. Let $\{{\mathbf{X}}_t\}_{t\ge1}$ be the sequence generated by the low rank MEG algorithm, Algorithm \ref{alg:LRMD}, on Problem \eqref{eq:Model} with $n\not=r\ge r^*=\textrm{rank}({\mathbf{X}}^*)$, $\eta_t\le\frac{1}{\beta}$ and $\varepsilon_t\in(0,3/4]$ such that $\lambda_{r^*+1}({\mathbf{X}}_t)\le\frac{\varepsilon_{t}}{n-r}$ for all $t\ge1$. Let ${\mathbf{Z}}_{t+1}$ be the standard MEG step update to ${\mathbf{X}}_t$, as defined in \eqref{eq:MDvonNeumannUpdate}. Then for any $t\ge1$, if ${\mathbf{X}}_t$ satisfies the condition in \eqref{eq:StochasticRadiusOfconvergence} with $\xi=0$, it holds that
\[ B({\mathbf{X}}^*,{\mathbf{X}}_{t+1}) \le B({\mathbf{X}}^*,{\mathbf{X}}_{t})+2\varepsilon_t. \]
\end{lemma}

\begin{proof}
From the optimality of ${\mathbf{X}}^*$, $\beta$-smoothness of $f$ and plugging in \eqref{ineq:strongConvexityOfBregmanDistance} due to the $1$-strong convexity of $B(\cdot,\cdot)$, it holds that
\begin{align*}
f({\mathbf{X}}^*) & \le f({\mathbf{Z}}_{t+1}) \le f({\mathbf{X}}_t)+\langle \nabla f({\mathbf{X}}_t),{\mathbf{Z}}_{t+1}-{\mathbf{X}}_t\rangle+\frac{\beta}{2}\Vert{\mathbf{Z}}_{t+1}-{\mathbf{X}}_t\Vert_F^2
\\ & \le f({\mathbf{X}}_t)+\langle \nabla f({\mathbf{X}}_t),{\mathbf{Z}}_{t+1}-{\mathbf{X}}_t\rangle+\frac{\beta}{2}\Vert{\mathbf{Z}}_{t+1}-{\mathbf{X}}_t\Vert_*^2
\\ & \le f({\mathbf{X}}_t)+\langle \nabla f({\mathbf{X}}_t),{\mathbf{Z}}_{t+1}-{\mathbf{X}}_t\rangle+\betaB({\mathbf{Z}}_{t+1},{\mathbf{X}}_{t}).
\end{align*}

Therefore,
\begin{align} \label{ineq:partOfFejerProof}
\langle \nabla f({\mathbf{X}}_t),{\mathbf{Z}}_{t+1}-{\mathbf{X}}_t\rangle \ge f({\mathbf{X}}^*)-f({\mathbf{X}}_t)-\betaB({\mathbf{Z}}_{t+1},{\mathbf{X}}_{t}).
\end{align}

By the definition of ${\mathbf{Z}}_{t+1}$ in \eqref{eq:MDGeneralUpdate},
\[ {\mathbf{Z}}_{t+1}=\argmin_{{\mathbf{Z}}\in\mathcal{S}_n}\{\langle\eta_t\nabla f({\mathbf{X}}_t)-\nabla\omega({\mathbf{X}}_t),{\mathbf{Z}}\rangle+\omega({\mathbf{Z}})\}. \]

Therefore, by the optimality condition for ${\mathbf{Z}}_{t+1}$ it holds that
\begin{align*}
0 & \ge \langle -\eta_t\nabla f({\mathbf{X}}_t)+\nabla\omega({\mathbf{X}}_t)-\nabla\omega({\mathbf{Z}}_{t+1}),{\mathbf{X}}^*-{\mathbf{Z}}_{t+1}\rangle
\\ & = \langle \nabla\omega({\mathbf{X}}_t)-\nabla\omega({\mathbf{Z}}_{t+1}),{\mathbf{X}}^*-{\mathbf{Z}}_{t+1}\rangle + \langle \eta_t\nabla f({\mathbf{X}}_t),{\mathbf{Z}}_{t+1}-{\mathbf{X}}_t\rangle + \langle \eta_t\nabla f({\mathbf{X}}_t),{\mathbf{X}}_t-{\mathbf{X}}^*\rangle
\\ & \underset{(a)}{=} B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})+B({\mathbf{Z}}_{t+1},{\mathbf{X}}_{t})-B({\mathbf{X}}^*,{\mathbf{X}}_{t}) + \langle \eta_t\nabla f({\mathbf{X}}_t),{\mathbf{Z}}_{t+1}-{\mathbf{X}}_t\rangle\\ & \ \ \  + \langle \eta_t\nabla f({\mathbf{X}}_t),{\mathbf{X}}_t-{\mathbf{X}}^*\rangle
\\ & \underset{(b)}{\ge} B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})+B({\mathbf{Z}}_{t+1},{\mathbf{X}}_{t})-B({\mathbf{X}}^*,{\mathbf{X}}_{t}) + \eta_t\left[f({\mathbf{X}}^*)-f({\mathbf{X}}_t)-\betaB({\mathbf{Z}}_{t+1},{\mathbf{X}}_{t})\right]
\\ & \ \ \ + \langle \eta_t\nabla f({\mathbf{X}}_t),{\mathbf{X}}_t-{\mathbf{X}}^*\rangle
\\ & \underset{(c)}{\ge} B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})-B({\mathbf{X}}^*,{\mathbf{X}}_{t}) + \eta_t\left(f({\mathbf{X}}^*)-f({\mathbf{X}}_t) + \langle\nabla f({\mathbf{X}}_t),{\mathbf{X}}_t-{\mathbf{X}}^*\rangle\right).
\end{align*}
Here (a) follows from the three point lemma in \eqref{lemma:threePointLemma}, (b) follows from plugging in \eqref{ineq:partOfFejerProof}, and (c) holds since $\eta_t\le\frac{1}{\beta}$.

From the gradient inequality $f({\mathbf{X}}^*)-f({\mathbf{X}}_t) + \langle\nabla f({\mathbf{X}}_t),{\mathbf{X}}_t-{\mathbf{X}}^*\rangle\ge0$, and therefore,
\[ B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1}) \le B({\mathbf{X}}^*,{\mathbf{X}}_{t}). \]

If ${\mathbf{X}}_t$ satisfies the condition in \eqref{eq:StochasticRadiusOfconvergence} with $\xi=0$ and an exact gradient, then from Lemma \ref{lemma:StochasticRadius}
\begin{align*}
B({\mathbf{X}}^*,{\mathbf{X}}_{t+1}) = B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})+B({\mathbf{X}}^*,{\mathbf{X}}_{t+1})-B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1}) \le B({\mathbf{X}}^*,{\mathbf{X}}_{t})+2\varepsilon_t.
\end{align*}

\end{proof}

We can now finally derive a concrete convergence rate for Algorithm \ref{alg:LRMD} with exact gradients. For simplicity and ease of presentation, in the following theorem we consider the case in which the SVD rank parameter in Algorithm \ref{alg:LRMD} is set exactly to that of a corresponding optimal solution, i.e., $r=\textrm{rank}({\mathbf{X}}^*)$. In the sequel (see Theorem \ref{thm:smoothCaseGen}) we consider the more general case in which $r\geq \textrm{rank}({\mathbf{X}}^*)$.

\begin{theorem}\label{thm:smoothCaseDeminish}
Fix an optimal solution ${\mathbf{X}}^*$ to Problem \eqref{eq:Model} for which Assumption \ref{ass:strictcomp} holds with parameter $\delta>0$ and denote $r^*=\textrm{rank}({\mathbf{X}}^*)$.
Let $\{{\mathbf{X}}_t\}_{t\geq 1}$ be the sequence of iterates generated by Algorithm \ref{alg:LRMD} with deterministic updates and with a fixed step-size $\eta_t=\eta=1/\beta$, and with SVD rank parameter $r=r^*$. Suppose that for all $t\geq 0$: $\varepsilon_t=\frac{\tilde{\varepsilon}_0}{2\max\{G^2,1\}}\frac{1}{(t+1+c)^3}$ for some $\tilde{\varepsilon}_0\leq R_0^2$,
where
\begin{align*}
R_0 := \frac{1}{4}\left[2\beta+\left(1+\frac{2\sqrt{2r^*}}{\lambda_{r^*}({\mathbf{X}}^*)}\right)G\right]^{-1}\delta,
\end{align*}
$c\ge 6\beta/\delta$, and $G\geq\sup_{{\mathbf{X}}\in\mathcal{S}_n}\Vert \nabla f({\mathbf{X}})\Vert_2$. Finally, assume the initialization matrix ${\mathbf{X}}_0$ satisfies $\textrm{rank}({\mathbf{X}}_0) = r^*$ and the parameters ${\mathbf{X}}_0,{\mathbf{X}}^*,\tilde{\varepsilon}_0,R_0$ satisfy the initialization condition in Eq. \eqref{eq:init:equalRank} (substituting ${\mathbf{X}}={\mathbf{X}}_0$, $\varepsilon=\tilde{\varepsilon}_0$, $R=R_0$).
Then, for any $T\geq 1$ it holds that
\begin{align*}
\min_{1\le t\le {T}}f({\mathbf{X}}_t) -f({\mathbf{X}}^*) & \le \left(\left(1+\frac{1}{2\max\{G^2,1\}}\right)\beta{}R_0^2+4R_0\right)\frac{1}{T}.
\end{align*}
\end{theorem}

\begin{proof}
We first observe that by our choice of the sequence $\{\varepsilon_t\}_{t\geq 0}$ we have that the sequence $\left\{-\frac{1}{\eta}\log\left(\frac{\varepsilon_{t-1}}{\varepsilon_{t}}\right)\right\}_{t\geq 1}$ is monotone non-decreasing and thus, for all $t\geq 1$:
\begin{align*}
 -\frac{1}{\eta}\log\left(\frac{\varepsilon_{t-1}}{\varepsilon_{t}}\right) &\ge -\frac{1}{\eta}\log\left(\frac{\varepsilon_{0}}{\varepsilon_{1}}\right)=-3\beta\log\left(\frac{c+2}{c+1}\right)\\
 &\geq -3\beta\log\left(1 + \frac{1}{c}\right) \geq -\frac{3\beta}{c} \geq -\frac{\delta}{2},
\end{align*}
where the last inequality follows from plugging our choice for $c$.

Thus, we have that for all $t\geq 1$:
\begin{align}\label{eq:thm:smoothConv:1}
&\frac{1}{2\sqrt{2}}\left[2\beta+\left(1+\frac{2\sqrt{2r^*}}{\lambda_{r^*}({\mathbf{X}}^*)}\right)G\right]^{-1}\delta \nonumber
\\ & \le \frac{1}{\sqrt{2}}\left[2\beta+\left(1+\frac{2\sqrt{2r^*}}{\lambda_{r^*}({\mathbf{X}}^*)}\right)G\right]^{-1}\left(\delta-\frac{1}{\eta}\log\left(\frac{\varepsilon_{t-1}}{\varepsilon_{t}}\right)\right) .
\end{align}

Thus, in order to invoke Lemma \ref{lemma:StochasticRadius} for all $t\geq 1$ it suffices to prove that for all $t\geq 1$: i. $\lambda_{r^*+1}({\mathbf{X}}_t) \leq \frac{\varepsilon_{t-1}}{n-r}$ and ii. $\sqrt{B({\mathbf{X}}^*,{\mathbf{X}}_t)} \leq \textrm{LHS of } \eqref{eq:thm:smoothConv:1}$.

The requirement $\lambda_{r^*+1}({\mathbf{X}}_t) \leq \frac{\varepsilon_{t-1}}{n-r}$ holds trivially by the design of the algorithm and since the SVD parameter satisfies $r=r^*$.

We now prove by induction that indeed for all $t\geq 1$, $\sqrt{B({\mathbf{X}}^*,{\mathbf{X}}_t)} \leq \textrm{LHS of } \eqref{eq:thm:smoothConv:1}$. The base case $t=1$ clearly holds due to the choice of initialization. Now, if the assumption holds for all $i\in\{1,...t\}$, then invoking Lemma \ref{lemma:bregmanDecrease} for all $i=1,...,t$ we have that
\begin{align} \label{ineq:recursionOfRadious}
B({\mathbf{X}}^*,{\mathbf{X}}_{t+1}) & \le B({\mathbf{X}}^*,{\mathbf{X}}_{t})+2\varepsilon_{t} \le \ldots \le B({\mathbf{X}}^*,{\mathbf{X}}_{1})+2\sum_{i=1}^{t}\varepsilon_i \nonumber
\\ &  \leq B({\mathbf{X}}^*,{\mathbf{X}}_{1})+\frac{\tilde{\varepsilon}_0}{\max\{G^2,1\}}\sum_{i=1}^{\infty}\frac{1}{(i+1)^3}  \le B({\mathbf{X}}^*,{\mathbf{X}}_{1})+\tilde{\varepsilon}_0 \leq 2R_0^2,
\end{align}
where the last inequality follows from our initialization assumption and Lemma \ref{lem:warmstart}.

Thus, we obtain that
\begin{align*}
\sqrt{B({\mathbf{X}}^*,{\mathbf{X}}_{t+1})} & \le \sqrt{2}R_0  = \frac{1}{2\sqrt{2}}\left[2\beta+\left(1+\frac{2\sqrt{2r}}{\lambda_r({\mathbf{X}}^*)}\right)G\right]^{-1}\delta,
\end{align*}
and the induction holds.

Invoking Lemma \ref{lemma:StochasticRadius} for all $t\geq 1$ guarantees that for all $t\ge 1$, $B({\mathbf{Z}}_{t+1},{\mathbf{X}}_{t+1})\le 2\varepsilon_t$ and $B({\mathbf{X}}^*,{\mathbf{X}}_{t+1})-B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})\le 2\varepsilon_t$. Plugging-in $\eta=1/\beta$ and our initialization assumption, it holds from Theorem \ref{thm:converganceSmooth} that for all $T\geq 1$,
\begin{align*}
\min_{1\le t\le {T}}f({\mathbf{X}}_t) -f({\mathbf{X}}^*) & \le \frac{\beta{}R_0^2}{T}+\frac{2\beta}{T}\sum_{t=1}^{T-1}\varepsilon_t+\frac{2\sqrt{2}G}{T}\sum_{t=1}^{T-1}\sqrt{\varepsilon_t}.
\end{align*}

In order to bound the RHS of this inequality, we note that the following inequalities hold:
\begin{align*}
\sum_{t=1}^{T-1}\sqrt{\varepsilon_t} &\leq \sum_{t=1}^{T-1}\frac{R_0}{\sqrt{2}\max\{G,1\}}\frac{1}{(t+1+c)^{3/2}} \\& < \frac{R_0}{\sqrt{2}\max\{G,1\}}\int_0^{\infty}\frac{1}{(t+1)^{3/2}}dt = \frac{\sqrt{2}R_0}{\max\{G,1\}}.\\
\sum_{t=1}^{T-1}\varepsilon_t &\leq \sum_{t=1}^{T-1}\frac{R_0^2}{2\max\{G^2,1\}}\frac{1}{(t+1+c)^{3}} < \frac{R_0^2}{2\max\{G^2,1\}}\int_0^{\infty}\frac{dt}{(t+1)^{3}} = \frac{R_0^2}{4\max\{G^2,1\}}.
\end{align*}

Plugging-in these  bounds we obtain the rate in the theorem.

\end{proof}

One caveat of Theorem \ref{thm:smoothCaseDeminish} is that it requires exact knowledge of the rank of the optimal solution since it requires to set $r=\textrm{rank}({\mathbf{X}}^*)$ in Algorithm \ref{alg:LRMD}. We will now present an analysis that removes this limitation.

Towards this end, we observe that from Lemma \ref{lemma:StochasticRadius}, in order for any $n\not=r>r^*$ to be suitable, it needs to satisfy the condition $\lambda_{r^*+1}({\mathbf{X}}_t)\le\frac{\varepsilon_{t-1}}{n-r}$ (in the case of $r=r^*$ we have shown that it holds trivially). Thus,  we begin by upper-bounding $\lambda_{r^*+1}({\mathbf{X}}_{t+1})$.
\begin{lemma}
Let ${\mathbf{X}}^*\in\mathcal{X}^*$ be an optimal solution which satisfies Assumption \ref{ass:strictcomp} with some $\delta>0$ and denote $r^*=\textrm{rank}({\mathbf{X}}^*)$. Suppose the rank parameter in Algorithm \ref{alg:LRMD} satisfies $n\not=r>r^*$. Denote $G=\sup_{{\mathbf{X}}\in\mathcal{S}_n}\Vert \nabla f({\mathbf{X}})\Vert_2$. Then for any $t\ge1$ it holds that
\begin{align} \label{ineq:boundOnLambda_rstarplus1}
 & \lambda_{r^*+1}({\mathbf{X}}_{t+1}) \nonumber \\ & \le \max\left\lbrace(1-\varepsilon_t)\lambda_{r^*+1}({\mathbf{X}}_t)\exp\left(-\eta_t\delta+\eta_t c \sqrt{B({\mathbf{X}}^*,{\mathbf{X}}_t)}+B({\mathbf{X}}^*,{\mathbf{X}}_t)\right),\frac{\varepsilon_t}{n-r}\right\rbrace,
\end{align}
where $c=\sqrt{2}\left(2\beta+\frac{2\sqrt{2{r^*}}}{\lambda_{r^*}({\mathbf{X}}^*)}G\right)$.

\end{lemma}

\begin{proof}

We have seen in \eqref{ineq:StochasticInProofSecondTermWithoutAssumption} with $\xi=0$ that
\begin{align*}
& \log(\lambda_{r^*+1}({\mathbf{Y}}_t)) \\ & \le \log\left(\lambda_{r^*+1}({\mathbf{X}}_t)\right)-\eta_t \lambda_{n-r^*}(\nabla f({\mathbf{X}}^*))+\eta_t\frac{2\sqrt{2r^*}G\Vert{\mathbf{X}}_t-{\mathbf{X}}^*\Vert_{*}}{\lambda_{r^*}({\mathbf{X}}^*)} \nonumber
+\eta_t\beta\Vert{\mathbf{X}}_t-{\mathbf{X}}^*\Vert_*.
\end{align*}

In addition,
\begin{align*}
-\log\left(\sum_{i=1}^r\lambda_i({\mathbf{Y}}_t)\right) & \le -\log\left(\sum_{i=1}^{r^*}\lambda_i({\mathbf{Y}}_t)\right) = -\log\left(\sum_{i=1}^{r^*}\lambda_i\left(\exp\left(\log({\mathbf{X}}_t)-\eta_t\nabla f({\mathbf{X}}_t)\right)\right)\right)
\\ & \underset{(a)}{\le} -\sum_{i=1}^{r^*}\lambda_i({\mathbf{X}}^*)\cdot\log\left(\frac{\lambda_i\left(\exp\left(\log({\mathbf{X}}_t)-\eta_t\nabla f({\mathbf{X}}_t)\right)\right)}{\lambda_i({\mathbf{X}}^*)}\right)
\\ & = -\sum_{i=1}^{r^*}\lambda_i({\mathbf{X}}^*)\cdot\left[\log\left(\lambda_i\left(\exp\left(\log({\mathbf{X}}_t)-\eta_t\nabla f({\mathbf{X}}_t)\right)\right)-\log(\lambda_i({\mathbf{X}}^*)\right)\right]
\\ & = -\sum_{i=1}^{r^*}\lambda_i({\mathbf{X}}^*)\lambda_i\left(\log({\mathbf{X}}_t)-\eta_t\nabla f({\mathbf{X}}_t)\right)+\sum_{i=1}^{r^*}\lambda_i({\mathbf{X}}^*)\log(\lambda_i({\mathbf{X}}^*))
\\ & = -\sum_{i=1}^{n}\lambda_i({\mathbf{X}}^*)\lambda_i\left(\log({\mathbf{X}}_t)-\eta_t\nabla f({\mathbf{X}}_t)\right)+\sum_{i=1}^{r^*}\lambda_i({\mathbf{X}}^*)\log(\lambda_i({\mathbf{X}}^*))
\\ & \underset{(b)}{\le} -\langle{\mathbf{X}}^*,\log({\mathbf{X}}_t)-\eta_t\nabla f({\mathbf{X}}_t)\rangle +\sum_{i=1}^{r^*}\lambda_i({\mathbf{X}}^*)\lambda_i(\log({\mathbf{X}}^*))
\\ & \underset{(c)}{=} -\langle{\mathbf{X}}^*,\log({\mathbf{X}}_t)-\eta_t\nabla f({\mathbf{X}}_t)\rangle +\sum_{i=1}^{r^*}\lambda_i({\mathbf{X}}^*\log({\mathbf{X}}^*))
\\ & = -\textrm{Tr}({\mathbf{X}}^*\log({\mathbf{X}}_t))+\textrm{Tr}({\mathbf{X}}^*\log({\mathbf{X}}^*))+\eta_t\langle{\mathbf{X}}^*,\nabla f({\mathbf{X}}_t)\rangle
\\ & = B({\mathbf{X}}^*,{\mathbf{X}}_t)+\eta_t\langle{\mathbf{X}}^*,\nabla f({\mathbf{X}}^*)\rangle + \eta_t\langle{\mathbf{X}}^*,\nabla f({\mathbf{X}}_t)-\nabla f({\mathbf{X}}^*)\rangle
\\ & = B({\mathbf{X}}^*,{\mathbf{X}}_t)+\eta_t\lambda_{n}(\nabla f({\mathbf{X}}^*)) + \eta_t\langle{\mathbf{X}}^*,\nabla f({\mathbf{X}}_t)-\nabla f({\mathbf{X}}^*)\rangle
\\ & \le B({\mathbf{X}}^*,{\mathbf{X}}_t)+\eta_t\lambda_{n}(\nabla f({\mathbf{X}}^*)) +\eta_t\Vert{\mathbf{X}}^*\Vert_*\cdot\Vert\nabla f({\mathbf{X}}_t)-\nabla f({\mathbf{X}}^*)\Vert_2
\\ & \underset{(d)}{\le} B({\mathbf{X}}^*,{\mathbf{X}}_t)+\eta_t\lambda_{n}(\nabla f({\mathbf{X}}^*)) +\eta_t\Vert\nabla f({\mathbf{X}}_t)-\nabla f({\mathbf{X}}^*)\Vert_2
\\ & \le B({\mathbf{X}}^*,{\mathbf{X}}_t)+\eta_t\lambda_{n}(\nabla f({\mathbf{X}}^*))+\eta_t\beta\Vert{\mathbf{X}}_t-{\mathbf{X}}^*\Vert_*
\end{align*}
where (a) follows since $-\log(\cdot)$ is convex, (b) follows from the von Neumann inequality, (c) follows since ${\mathbf{X}}_t$ and $\log({\mathbf{X}}_t)$ have the same eigen-vectors, and (d) follows since $\Vert{\mathbf{X}}^*\Vert_*=1$.

Then
\begin{align*}
\log\left(\frac{\lambda_{r^*+1}({\mathbf{Y}}_t)}{a_t}\right) & = \log(\lambda_{r^*+1}({\mathbf{Y}}_t))-\log\left(\sum_{i=1}^{r^*}\lambda_i({\mathbf{Y}}_t)\right)
\\ & \le \log\left(\lambda_{r^*+1}({\mathbf{X}})\right)-\eta_t \delta +2\eta_t\left(\beta+\frac{\sqrt{2r^*}G}{\lambda_{r^*}({\mathbf{X}}^*)}\right)\Vert{\mathbf{X}}-{\mathbf{X}}^*\Vert_*+B({\mathbf{X}}^*,{\mathbf{X}})
\\ & \le \log\left(\lambda_{r^*+1}({\mathbf{X}}_t)\right)-\eta_t\delta+\eta_t c \sqrt{B({\mathbf{X}}^*,{\mathbf{X}}_t)}+B({\mathbf{X}}^*,{\mathbf{X}}_t),
\end{align*}
where the last inequality holds from \eqref{ineq:strongConvexityOfBregmanDistance} and $c=\sqrt{2}\left(2\beta+\frac{2\sqrt{2{r^*}}}{\lambda_{r^*}({\mathbf{X}}^*)}G\right)$.

Therefore,
\begin{align*}
\frac{\lambda_{r^*+1}({\mathbf{Y}}_t)}{a_t} & \le \exp\left(\log\left(\lambda_{r^*+1}({\mathbf{X}}_t)\right)-\eta_t\delta+\eta_t c \sqrt{B({\mathbf{X}}^*,{\mathbf{X}}_t)}+B({\mathbf{X}}^*,{\mathbf{X}}_t)\right)
\\ & = \lambda_{r^*+1}({\mathbf{X}}_t)\exp\left(-\eta_t\delta+\eta_t c \sqrt{B({\mathbf{X}}^*,{\mathbf{X}}_t)}+B({\mathbf{X}}^*,{\mathbf{X}}_t)\right).
\end{align*}

Using the definition of the update of ${\mathbf{X}}_{t+1}$,
\begin{align*}
& \lambda_{r^*+1}({\mathbf{X}}_{t+1})  = \max\left\lbrace(1-\varepsilon_t)\frac{\lambda_{r^*+1}({\mathbf{Y}}_t)}{a_t},\frac{\varepsilon_t}{n-r}\right\rbrace
\\ & \le \max\left\lbrace(1-\varepsilon_t)\lambda_{r^*+1}({\mathbf{X}}_t)\exp\left(-\eta_t\delta+\eta_t c \sqrt{B({\mathbf{X}}^*,{\mathbf{X}}_t)}+B({\mathbf{X}}^*,{\mathbf{X}}_t)\right),\frac{\varepsilon_t}{n-r}\right\rbrace.
\end{align*}

\end{proof}

\begin{theorem}\label{thm:smoothCaseGen}Fix an optimal solution ${\mathbf{X}}^*$ to Problem \eqref{eq:Model} for which Assumption \ref{ass:strictcomp} holds with parameter $\delta>0$ and denote $r^*=\textrm{rank}({\mathbf{X}}^*)$. Let $\{{\mathbf{X}}_t\}_{t\geq 1}$ be the sequence of iterates generated by Algorithm \ref{alg:LRMD} with deterministic updates and with a fixed step-size $\eta_t=\eta=1/\beta$, and with SVD rank parameter $n\not=r\geq{}r^*$. Suppose that for all $t\geq 0$: $\varepsilon_t=\frac{3\tilde{\varepsilon}_0}{2\max\{G^2,1\}}\frac{1}{(t+c+1)^3}$ for some $\tilde{\varepsilon}_0\leq R_0^2$,
where
\begin{align*}
R_0 &:=\min\left\lbrace\frac{1}{4}\left[2\beta+\left(1+\frac{2\sqrt{2r^*}}{\lambda_{r^*}({\mathbf{X}}^*)}\right)G\right]^{-1}\delta,~\frac{G}{\beta}\right\rbrace,
\end{align*}
$c\ge\frac{6\max\{\beta,1\}}{\delta}$, and $G\geq\sup_{{\mathbf{X}}\in\mathcal{S}_n}\Vert \nabla f({\mathbf{X}})\Vert_2$. Finally, assume the initialization matrix ${\mathbf{X}}_0$ satisfies $\lambda_{r^*+1}({\mathbf{X}}_0) \leq r\varepsilon_0/(n(n-r))$ and the parameters ${\mathbf{X}}_0,{\mathbf{X}}^*,\tilde{\varepsilon}_0,R_0$ satisfy the initialization condition in Eq. \eqref{eq:init:highRank} (substituting ${\mathbf{X}}={\mathbf{X}}_0$, $\varepsilon=\tilde{\varepsilon}_0$, $R=R_0$).
Then, for any $T\geq 1$ it holds that
\begin{align*}
\min_{1\le t\le {T}}f({\mathbf{X}}_t) -f({\mathbf{X}}^*) & \le \left(\left(1+\frac{1}{\max\{G^2,1\}}\right)\beta R_0^2+4\sqrt{3}R_0\right)\frac{1}{T}.
\end{align*}
\end{theorem}

\begin{proof}

We first observe that by our choice of the sequence $\lbrace \varepsilon_t\rbrace_{t\ge0}$ and since $\eta_t$ is fixed for all $t$, we have that the sequence $\left\lbrace-\max\left\lbrace\frac{1}{\eta_{t}},1\right\rbrace\log\left(\frac{\varepsilon_{t-1}}{\varepsilon_{t}}\right)\right\rbrace_{t\ge1}$ is monotone non-decreasing, and thus, for all $t\ge1$:
\begin{align*}
-\max\left\lbrace\frac{1}{\eta_{t}},1\right\rbrace\log\left(\frac{\varepsilon_{t-1}}{\varepsilon_{t}}\right) & \ge -\max\left\lbrace\frac{1}{\eta_{1}},1\right\rbrace\log\left(\frac{\varepsilon_{0}}{\varepsilon_{1}}\right)
= -3\max\lbrace\beta,1\rbrace\log\left(\frac{c+2}{c+1}\right)
\\ & \ge -3\max\lbrace\beta,1\rbrace\log\left(1+\frac{1}{c}\right)
\ge -\frac{3\max\lbrace\beta,1\rbrace}{c} \ge -\frac{\delta}{2},
\end{align*}
where the last inequality follows from plugging in our choice for $c$.

Thus, we have that for all $t\geq 1$:
\begin{align}\label{eq:thm:smoothConv:4}
& \min\left\lbrace\frac{1}{\sqrt{2}}\left[2\beta+\left(1+\frac{2\sqrt{2r^*}}{\lambda_{r^*}({\mathbf{X}}^*)}\right)G\right]^{-1}\left(\delta-\max\left\lbrace\frac{1}{\eta_{t}},1\right\rbrace\log\left(\frac{\varepsilon_{t-1}}{\varepsilon_{t}}\right)\right),\sqrt{2}\eta G\right\rbrace \nonumber
\\ & \ge \min\left\lbrace\frac{1}{\sqrt{2}}\left[2\beta+\left(1+\frac{2\sqrt{2r^*}}{\lambda_{r^*}({\mathbf{X}}^*)}\right)G\right]^{-1}\frac{\delta}{2},\sqrt{2}\eta G\right\rbrace.
\end{align}

Thus, in order to invoke Lemma \ref{lemma:StochasticRadius} for all $t\geq 1$ it suffices to prove that for all $t\geq 1$: i. $\lambda_{r^*+1}({\mathbf{X}}_t) \leq \frac{\varepsilon_{t-1}}{n-r}$ and ii. $\sqrt{B({\mathbf{X}}^*,{\mathbf{X}}_t)} \leq \textrm{RHS of } \eqref{eq:thm:smoothConv:4}$.

We now prove by induction that indeed for all $t\geq 1$, $\sqrt{B({\mathbf{X}}^*,{\mathbf{X}}_t)} \leq \textrm{RHS of } \eqref{eq:thm:smoothConv:4}$ and $\lambda_{r^*+1}({\mathbf{X}}_{t})\le\frac{\varepsilon_t}{n-r}$. The base case holds due to our initialization, and by noticing that from the definition of ${\mathbf{X}}_1$ the condition $\lambda_{r^*+1}({\mathbf{X}}_0) \leq r\varepsilon_0/(n(n-r))$ implies that $\lambda_{r^*+1}({\mathbf{X}}_1)\le \frac{r\varepsilon_0}{n(n-r)}+\frac{\varepsilon_0}{n}=\frac{\varepsilon_0}{n-r}$. Now, if the assumptions holds for all $i\in\lbrace 1,\ldots,t-1\rbrace$, then from \eqref{ineq:boundOnLambda_rstarplus1}
\begin{align*}
 & \lambda_{r^*+1}({\mathbf{X}}_{t})
\\ & \le \max\left\lbrace(1-\varepsilon_t)\lambda_{r^*+1}({\mathbf{X}}_{t-1})\exp\left(-\eta_t\delta+\eta_t c \sqrt{B({\mathbf{X}}^*,{\mathbf{X}}_{t-1})}+B({\mathbf{X}}^*,{\mathbf{X}}_{t-1})\right),\frac{\varepsilon_t}{n-r}\right\rbrace
\\ & \le \max\left\lbrace\frac{\varepsilon_{t-1}}{n-r}\exp\left(-\eta_t\delta+\eta_t c \sqrt{B({\mathbf{X}}^*,{\mathbf{X}}_{t-1})}+B({\mathbf{X}}^*,{\mathbf{X}}_{t-1})\right),\frac{\varepsilon_t}{n-r}\right\rbrace
\\ & \underset{(a)}{\le} \max\left\lbrace\frac{\varepsilon_{t-1}}{n-r}\exp\left(-\eta_t\delta+\eta_t c \sqrt{B({\mathbf{X}}^*,{\mathbf{X}}_{t-1})}+\sqrt{2}\eta_t G\sqrt{B({\mathbf{X}}^*,{\mathbf{X}}_{t-1})}\right),\frac{\varepsilon_t}{n-r}\right\rbrace
\\ & \underset{(b)}{\le} \max\left\lbrace\frac{\varepsilon_{t-1}}{n-r}\exp\left(-\max\left\lbrace\frac{1}{\eta_{t-1}},1\right\rbrace\log\left(\frac{\varepsilon_{t-2}}{\varepsilon_{t-1}}\right)\right),\frac{\varepsilon_{t}}{n-r}\right\rbrace
\\ & \le \max\left\lbrace\frac{\varepsilon_{t-1}}{n-r}\exp\left(-\log\left(\frac{\varepsilon_{t-2}}{\varepsilon_{t-1}}\right)\right),\frac{\varepsilon_{t}}{n-r}\right\rbrace
\\ & \le \max\left\lbrace\frac{\varepsilon_{t-1}}{n-r}\exp\left(-\log\left(\frac{\varepsilon_{t-1}}{\varepsilon_{t}}\right)\right),\frac{\varepsilon_{t}}{n-r}\right\rbrace = \frac{\varepsilon_t}{n-r},
\end{align*}
where (a) and (b) hold due to the bound in the induction hypothesis of $\textrm{the RHS of } \eqref{eq:thm:smoothConv:4}$.
Invoking Lemma \ref{lemma:StochasticRadius}, it holds that $B({\mathbf{X}}^*,{\mathbf{X}}_{i}) \le B({\mathbf{X}}^*,{\mathbf{Z}}_{i})+2\varepsilon_i$ for all $i \in \lbrace1,\ldots,t\rbrace$.
Using recursion and invoking Lemma \ref{lemma:bregmanDecrease}, we obtain
\begin{align*}
B({\mathbf{X}}^*,{\mathbf{X}}_{t}) & \le B({\mathbf{X}}^*,{\mathbf{X}}_{t-1})+2\varepsilon_{t-1} \le \ldots \le B({\mathbf{X}}^*,{\mathbf{X}}_{1})+2\sum_{i=1}^{t-1}\varepsilon_i \nonumber
\\ & \le B({\mathbf{X}}^*,{\mathbf{X}}_{1})+\frac{3\tilde{\varepsilon}_0}{\max\{G^2,1\}}\sum_{i=1}^{T}\frac{1}{(t+c+1)^3}
\\ & \le B({\mathbf{X}}^*,{\mathbf{X}}_{1})+3\tilde{\varepsilon}_0\sum_{i=1}^{\infty}\frac{1}{(t+1)^3} \le B({\mathbf{X}}^*,{\mathbf{X}}_{1})+\tilde{\varepsilon}_0\le 2R_0^2,
\end{align*}
where the last inequality holds from our initialization assumption and Lemma \ref{lem:warmstart}.

Thus, we obtain that for all $t\ge1$:
\begin{align*}
& \sqrt{B({\mathbf{X}}^*,{\mathbf{X}}_t)}\le \sqrt{2}R_0 = \min\left\lbrace\frac{1}{\sqrt{2}}\left[2\beta+\left(1+\frac{2\sqrt{2r^*}}{\lambda_{r^*}({\mathbf{X}}^*)}\right)G\right]^{-1}\frac{\delta}{2},\sqrt{2}\eta G\right\rbrace,
\end{align*}
and the induction holds.

Invoking Lemma \ref{lemma:StochasticRadius} for all $t\ge1$ guarantees that for all $t\ge1$,
$B({\mathbf{X}}^*,{\mathbf{X}}_{t+1}) - B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1}) \le 2\varepsilon_t$ and $B({\mathbf{Z}}_{t+1},{\mathbf{X}}_{t+1})\le 2\varepsilon_t$.

In addition, we note that the following inequalities hold for any $T\ge 1$:
\[ \sum_{t=1}^{T-1}\sqrt{\varepsilon_t}=\frac{\sqrt{3\tilde{\varepsilon}_0}}{\sqrt{2}\max\{G,1\}}\sum_{t=1}^{T-1} \frac{1}{(t+c+1)^{1.5}} \le \frac{\sqrt{3\tilde{\varepsilon}_0}}{\sqrt{2}\max\{G,1\}}\sum_{t=1}^{\infty} \frac{1}{(t+1)^{1.5}} < \frac{\sqrt{6}R_0}{\max\{G,1\}}, \]
\[ \sum_{t=1}^{T-1}\varepsilon_t=\frac{3\tilde{\varepsilon}_0}{2\max\{G^2,1\}}\sum_{t=1}^{T-1} \frac{1}{(t+c+1)^3}\le \frac{3\tilde{\varepsilon}_0}{2\max\{G^2,1\}}\sum_{t=1}^{\infty} \frac{1}{(t+1)^3} < \frac{R_0^2}{2\max\{G^2,1\}}.\]

Plugging in the last two bounds and our choice $\eta=\frac{1}{\beta}$ into Theorem \ref{thm:converganceSmooth}, we obtain that for all $T\ge1$,
\begin{align*}
\min_{1\le t\le {T}}f({\mathbf{X}}_t) -f({\mathbf{X}}^*) & \le \frac{\betaB({\mathbf{X}}^*,{\mathbf{X}}_{1})}{T}+\frac{2\beta}{T}\sum_{t=1}^{T-1}\varepsilon_t+\frac{2\sqrt{2}G}{T}\sum_{t=1}^{T-1}\sqrt{\varepsilon_t}
\\ & \le \frac{\beta R_0^2}{T}+\frac{\beta R_0^2}{\max\{G^2,1\}T}+\frac{4\sqrt{3}R_0}{T}.
\end{align*}

\end{proof}

\section{Provable Local Convergence with Stochastic Updates under Strict Complementarity} \label{section:stochasticCase}

In this section we turn to consider the stochastic setting in which $f(\cdot)$ is given by a first-order stochastic oracle that for any ${\mathbf{X}}\in\mathbb{S}^n$ it returns a random matrix $\widehat{\nabla}\in\mathbb{S}^n$ such that $\mathbb{E}[\widehat{\nabla}\vert{\mathbf{X}}]=\nabla f({\mathbf{X}})$, $\Vert\widehat{\nabla}\Vert\le G$ and $\textbf{Var}[\widehat{\nabla}\vert{\mathbf{X}}] = \mathbb{E}[\Vert{\widehat{\nabla}-\nabla{}f({\mathbf{X}})}\Vert^2|{\mathbf{X}}]\le\sigma^2$, for some $G>0$ and $\sigma>0$.

Our main result for this section is the proof that when initialized with a ``warm-start" point and suitable choice of parameters and assuming strict complementarity holds, Algorithm \ref{alg:LRMD} with mini-batches of stochastic gradients converges in expectation  to an optimal solution with rate $O(1/\sqrt{t})$.

Using Algorithm \ref{alg:LRMD} with mini-batches is crucial to the application of Lemma \ref{lemma:StochasticRadius}, in order to bound the gradient error parameter $\xi$ (at least with some positive probability). Thus, throughout this section we consider Algorithm \ref{alg:LRMD} with stochastic updates and with a fixed mini-batch size $L$, where the stochastic gradient on iteration $t$ is given by $\widehat{\nabla}_t = \frac{1}{L}\sum_{i=1}^L\widehat{\nabla}_t^{(i)}$, such that $\widehat{\nabla}_t^{(1)},\dots,\widehat{\nabla}_t^{(L)}$ denote $L$ i.i.d. queries to the stochastic oracle with the current iterate ${\mathbf{X}}_t$. Naturally, this is important since the variance of the mini-batched gradient satisfies $\textbf{Var}[\widehat{\nabla}_t|{\mathbf{X}}_t] = \sigma^2/L$.

In the following lemma we bound the decrease in the Bregman distance between an iterate and the expectation of its standard stochastic MEG update step computed upon it.

\begin{lemma} \label{lemma:StochasticBregmanDecrease}
Let $\{{\mathbf{X}}_t\}_{t\ge1}$ be the sequence generated by Algorithm \ref{alg:LRMD} with stochastic updates, and let ${\mathbf{X}}^*$ be an optimal solution. For all $t\geq 1$ let ${\mathbf{Z}}_{t+1}$ be the standard stochastic MEG ${\mathbf{X}}_t$, as defined in \eqref{eq:MDvonNeumannStochasticUpdate}. Then for any $t\ge1$,
it holds that
\[ \mathbb{E}[B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})\vert{\mathbf{X}}_t] \le B({\mathbf{X}}^*,{\mathbf{X}}_{t})+\frac{1}{2}\eta_t^2G^2. \]
\end{lemma}

\begin{proof}

By the definition of ${\mathbf{Z}}_{t+1}$ in \eqref{eq:MDGeneralUpdate}, where the exact gradient is replaced with a stochastic gradient $\widehat{\nabla}_t$,
\[ {\mathbf{Z}}_{t+1}=\argmin_{{\mathbf{Z}}\in\mathcal{S}_n}\{\langle\eta_t\widehat{\nabla}_t-\nabla\omega({\mathbf{X}}_t),{\mathbf{Z}}\rangle+\omega({\mathbf{Z}})\}. \]

Therefore, by the optimality condition for ${\mathbf{Z}}_{t+1}$, it holds that
\begin{align*}
0 & \ge \langle -\eta_t\widehat{\nabla}_t+\nabla\omega({\mathbf{X}}_t)-\nabla\omega({\mathbf{Z}}_{t+1}),{\mathbf{X}}^*-{\mathbf{Z}}_{t+1}\rangle
\\ & = \langle \nabla\omega({\mathbf{X}}_t)-\nabla\omega({\mathbf{Z}}_{t+1}),{\mathbf{X}}^*-{\mathbf{Z}}_{t+1}\rangle + \langle \eta_t\widehat{\nabla}_t,{\mathbf{Z}}_{t+1}-{\mathbf{X}}_t\rangle + \langle \eta_t\widehat{\nabla}_t,{\mathbf{X}}_t-{\mathbf{X}}^*\rangle
\\ & \underset{(a)}{=} B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})+B({\mathbf{Z}}_{t+1},{\mathbf{X}}_{t})-B({\mathbf{X}}^*,{\mathbf{X}}_{t}) + \langle \eta_t\widehat{\nabla}_t,{\mathbf{Z}}_{t+1}-{\mathbf{X}}_t\rangle + \langle \eta_t\widehat{\nabla}_t,{\mathbf{X}}_t-{\mathbf{X}}^*\rangle
\\ & \underset{(b)}{\ge} B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})-B({\mathbf{X}}^*,{\mathbf{X}}_{t}) +\langle \eta_t\widehat{\nabla}_t,{\mathbf{X}}_t-{\mathbf{X}}^*\rangle \\ & \ \ \ +\frac{1}{2}\Vert{\mathbf{Z}}_{t+1}-{\mathbf{X}}_t\Vert^2_* - \eta_t\Vert\widehat{\nabla}_t\Vert_2\Vert{\mathbf{Z}}_{t+1}-{\mathbf{X}}_t\Vert_*
\\ & \ge B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})-B({\mathbf{X}}^*,{\mathbf{X}}_{t}) +\langle \eta_t\widehat{\nabla}_t,{\mathbf{X}}_t-{\mathbf{X}}^*\rangle + \min_{a\in\mathbb{R}}\left\{\frac{1}{2}a^2-a\eta_t\Vert\widehat{\nabla}_t\Vert_{2}\right\}
\\ & = B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})-B({\mathbf{X}}^*,{\mathbf{X}}_{t}) +\langle \eta_t\widehat{\nabla}_t,{\mathbf{X}}_t-{\mathbf{X}}^*\rangle -\frac{1}{2}\eta_t^2\Vert\widehat{\nabla}_t\Vert_{2}^2
\\ & \underset{(c)}{\ge} B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})-B({\mathbf{X}}^*,{\mathbf{X}}_{t}) + \eta_t\left(f({\mathbf{X}}^*)-f({\mathbf{X}}_t) + \langle\widehat{\nabla}_t,{\mathbf{X}}_t-{\mathbf{X}}^*\rangle\right)-\frac{1}{2}\eta_t^2\Vert\widehat{\nabla}_t\Vert_{2}^2.
\end{align*}
where (a) holds from the three point lemma in \eqref{lemma:threePointLemma}, (b) follows from $1$-strong convexity of the Bregman distance as in \eqref{ineq:strongConvexityOfBregmanDistance} and H\"{o}lder's inequality, and (c) follows from the optimality of ${\mathbf{X}}^*$.

Taking expectation with respect to ${\mathbf{X}}_t$, we get
\begin{align*}
& \mathbb{E}[B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})\vert{\mathbf{X}}_t] \\ & \le B({\mathbf{X}}^*,{\mathbf{X}}_{t})-\eta_t\left(f({\mathbf{X}}^*)-f({\mathbf{X}}_t) + \langle\nabla f({\mathbf{X}}_t),{\mathbf{X}}_t-{\mathbf{X}}^*\rangle\right)+\frac{1}{2}\eta_t^2\mathbb{E}[\Vert\widehat{\nabla}_t\Vert_{2}^2\vert{\mathbf{X}}_t].
\end{align*}

Using the gradient inequality, this implies that
\[ \mathbb{E}[B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})\vert{\mathbf{X}}_t] \le B({\mathbf{X}}^*,{\mathbf{X}}_{t})+\frac{1}{2}\eta_t^2G^2. \]

\end{proof}

In our analysis of the deterministic case we have used Lemma \ref{lemma:bregmanDecrease} to show that the Bregman distance between the iterates and the optimal solution does not increase over time, up to some small easily controlled error resulting from the inexact low-rank computations. This property does not hold anymore  when using stochastic updates. Therefore, towards obtaining a convergence rate for the stochastic setting, we  introduce a martingale argument to prove that with high probability the iterates of Algorithm \ref{alg:LRMD} remain inside the ball around the optimal solution inside-which the convergence of the algorithm could be guaranteed.

\begin{lemma} \label{lemma:concentrationBound}
Fix $p\in(0,1)$. Let $\lbrace {\mathbf{X}}_t\rbrace_{t\in [T]}$ be the sequence generated by $T$ iterations of Algorithm \ref{alg:LRMD} with stochastic updates and with fixed mini-batch size $L>0$, and suppose that $\lbrace\eta_t\rbrace_{t\in [T-1]}$ is a non-increasing sequence of non-negative scalars. For any $t\ge1$ let ${\mathbf{Z}}_{t+1}$ be the standard stochastic MEG update to ${\mathbf{X}}_t$ as defined in \eqref{eq:MDvonNeumannStochasticUpdate} and ${\mathbf{Z}}_1={\mathbf{X}}_1$. Then for any ${\mathbf{X}}^*\in\mathcal{X}^*$, if the inequality
\begin{align} \label{ineq:conditionOfConcentrationBoundLemma}
\sqrt{\sum_{i=1}^{T-1}\eta_{i}^2\sigma^2/L}\ge \frac{1}{3}\eta_1G\sqrt{2\log\left(\frac{T}{p}\right)}
\end{align}
holds, then with probability at least $1-p$, it holds  for all $t\in[T]$ that
\begin{align*}
B({\mathbf{X}}^*,{\mathbf{X}}_t) & \le B({\mathbf{X}}^*,{\mathbf{X}}_1)+\frac{1}{2}G^2\sum_{i=1}^{t-1}\eta_i^2+\sum_{i=1}^{t}(B({\mathbf{X}}^*,{\mathbf{X}}_{i})-B({\mathbf{X}}^*,{\mathbf{Z}}_{i}))\\ & \ \ \ +\sqrt{16\sigma^2/L\sum_{i=1}^{T-1}\eta_{i}^2}\sqrt{\log\left(\frac{T}{p}\right)}.
\end{align*}
\end{lemma}

\begin{proof}

For all $t\in[T]$ define the random variable
\[ W_t:=B({\mathbf{X}}^*,{\mathbf{Z}}_t)-\frac{1}{2}G^2\sum_{i=1}^{t-1}\eta_i^2-\sum_{i=1}^{t-1}(B({\mathbf{X}}^*,{\mathbf{X}}_{i})-B({\mathbf{X}}^*,{\mathbf{Z}}_{i})). \]

$W_1,\ldots,W_T$ form a sub-martingale sequence with respect to the sequence ${\mathbf{X}}_1,\ldots,{\mathbf{X}}_{T-1}$, since for all $t\in[T-1]$ it holds that
\begin{align*}
& \mathbb{E}[W_{t+1}\vert {\mathbf{X}}_1,\ldots,{\mathbf{X}}_{t}] \\ & = \mathbb{E}[B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})\vert{\mathbf{X}}_{t}]-\frac{1}{2}G^2\sum_{i=1}^{t}\eta_i^2-\sum_{i=1}^{t}\mathbb{E}[B({\mathbf{X}}^*,{\mathbf{X}}_{i})-B({\mathbf{X}}^*,{\mathbf{Z}}_{i})\vert{\mathbf{X}}_i]
\\ & = \mathbb{E}[B({\mathbf{X}}^*,{\mathbf{Z}}_{t+1})\vert{\mathbf{X}}_{t}]-\frac{1}{2}G^2\sum_{i=1}^{t}\eta_i^2-\sum_{i=1}^{t}(B({\mathbf{X}}^*,{\mathbf{X}}_{i})-B({\mathbf{X}}^*,{\mathbf{Z}}_{i}))
\\ & \underset{(a)}{\le} B({\mathbf{X}}^*,{\mathbf{X}}_{t})+\frac{1}{2}\eta_t^2G^2-\frac{1}{2}G^2\sum_{i=1}^{t}\eta_i^2-\sum_{i=1}^{t}(B({\mathbf{X}}^*,{\mathbf{X}}_{i})-B({\mathbf{X}}^*,{\mathbf{Z}}_{i}))
\\ & = B({\mathbf{X}}^*,{\mathbf{Z}}_{t})-\frac{1}{2}G^2\sum_{i=1}^{t-1}\eta_i^2-\sum_{i=1}^{t-1}(B({\mathbf{X}}^*,{\mathbf{X}}_{i})-B({\mathbf{X}}^*,{\mathbf{Z}}_{i})) = W_t,
\end{align*}
where (a) follows from Lemma \ref{lemma:StochasticBregmanDecrease}.

We will show that the sub-martingale has bounded-differences to upper-bound its variance.
It holds for all $2\le t\le T$ that
\begin{align} \label{eq:boundedDifferences}
W_t-\mathbb{E}[W_{t}\vert W_1,\ldots,W_{t-1}] & = B({\mathbf{X}}^*,{\mathbf{Z}}_{t})-\mathbb{E}[B({\mathbf{X}}^*,{\mathbf{Z}}_{t})\vert{\mathbf{X}}_{t-1}].
\end{align}

First, note that for any $i\in[n]$ it holds that $\lambda_i\left(\log\left(\frac{{\mathbf{Y}}_{t-1}}{b_{t-1}}\right)\right)=\log\left(\lambda_i\left(\frac{{\mathbf{Y}}_{t-1}}{b_{t-1}}\right)\right)=\log\left(\frac{1}{b_{t-1}}\lambda_i\left({\mathbf{Y}}_{t-1}\right)\right)=\log\left(\lambda_i\left({\mathbf{Y}}_{t-1}\right)\right)-\log(b_{t-1})$.
Therefore, $\log\left(\frac{{\mathbf{Y}}_{t-1}}{b_{t-1}}\right)=\log\left({\mathbf{Y}}_{t-1}\right)-\log(b_{t-1}{\mathbf{I}})$. This implies that
\begin{align} \label{eq:traceYt/bt}
\textrm{Tr}\left({\mathbf{X}}^*\log\left({\mathbf{Z}}_t\right)\right) & =\textrm{Tr}\left({\mathbf{X}}^*\log\left(\frac{{\mathbf{Y}}_{t-1}}{b_{t-1}}\right)\right) = \textrm{Tr}\left({\mathbf{X}}^*\log\left({\mathbf{Y}}_{t-1}\right)\right)-\textrm{Tr}\left({\mathbf{X}}^*\log(b_{t-1}{\mathbf{I}})\right) \nonumber
\\ & = \textrm{Tr}\left({\mathbf{X}}^*\log\left({\mathbf{Y}}_{t-1}\right)\right)-\log(b_{t-1}).
\end{align}

Now, by the definition of ${\mathbf{Z}}_t$ in \eqref{eq:MDvonNeumannStochasticUpdate}, we obtain
\begin{align} \label{ineq:boundedDifferences1}
& B({\mathbf{X}}^*,{\mathbf{Z}}_{t})-\mathbb{E}[B({\mathbf{X}}^*,{\mathbf{Z}}_{t})\vert{\mathbf{X}}_{t-1}] \nonumber
\\ & = \mathbb{E}[\textrm{Tr}({\mathbf{X}}^*\log({\mathbf{Z}}_t))\vert{\mathbf{X}}_{t-1}]-\textrm{Tr}({\mathbf{X}}^*\log({\mathbf{Z}}_t)) \nonumber
\\ & \underset{(a)}{=} \mathbb{E}[\textrm{Tr}({\mathbf{X}}^*(\log({\mathbf{X}}_{t-1})-\eta_{t-1}\widehat{\nabla}_{t-1}))\vert{\mathbf{X}}_{t-1}]-\textrm{Tr}({\mathbf{X}}^*(\log({\mathbf{X}}_{t-1})-\eta_{t-1}\widehat{\nabla}_{t-1})) \nonumber \\ & \ \ \ -\mathbb{E}[\log(b_{t-1})\vert{\mathbf{X}}_{t-1}]+\log(b_{t-1}) \nonumber
\\ & = \eta_{t-1}\textrm{Tr}({\mathbf{X}}^*(\widehat{\nabla}_{t-1}-\nabla f({\mathbf{X}}_{t-1})))-\mathbb{E}[\log(b_{t-1})\vert{\mathbf{X}}_{t-1}]+\log(b_{t-1}) \nonumber
\\ & \underset{(b)}{\le} \eta_{t-1}\Vert{\mathbf{X}}^*\Vert_*\Vert\widehat{\nabla}_{t-1}-\nabla f({\mathbf{X}}_{t-1})\Vert_2-\mathbb{E}[\log(b_{t-1})\vert{\mathbf{X}}_{t-1}]+\log(b_{t-1}) \nonumber
\\ & = \eta_{t-1}\Vert\widehat{\nabla}_{t-1}-\nabla f({\mathbf{X}}_{t-1})\Vert_2-\mathbb{E}[\log(b_{t-1})\vert{\mathbf{X}}_{t-1}]+\log(b_{t-1}),
\end{align}
where (a) follows from \eqref{eq:traceYt/bt}, and (b) follows from H\"{o}lder's inequality.

In addition,
\begin{align} \label{ineq:boundedDifferences2}
\log(b_{t-1}) & = \log\left(\textrm{Tr}(\exp(\log({\mathbf{X}}_{t-1})-\eta_{t-1}\widehat{\nabla}_{t-1}))\right) \nonumber
\\ & \underset{(a)}{\le} \log\left(\textrm{Tr}(\exp(\log({\mathbf{X}}_{t-1}))\exp(-\eta_{t-1}\widehat{\nabla}_{t-1}))\right) \nonumber
\\ & =\log\left(\textrm{Tr}({\mathbf{X}}_{t-1}\exp(-\eta_{t-1}\widehat{\nabla}_{t-1}))\right) \nonumber
\\ & \underset{(b)}{\le} \log\left(\sum_{i=1}^{n}\lambda_i({\mathbf{X}}_{t-1})\lambda_i(\exp(-\eta_{t-1}\widehat{\nabla}_{t-1}))\right) \nonumber
\\ & = \log\left(\sum_{i=1}^{n}\lambda_i({\mathbf{X}}_{t-1})\exp(\lambda_i(-\eta_{t-1}\widehat{\nabla}_{t-1}))\right) \nonumber
\\ & \le \log\left(\sum_{i=1}^{n}\lambda_i({\mathbf{X}}_{t-1})\exp(\eta_{t-1}\Vert\widehat{\nabla}_{t-1}\Vert_2)\right)= \eta_{t-1}\Vert\widehat{\nabla}_{t-1}\Vert_2 \le \eta_{t-1}G,
\end{align}
where (a) follows from the Golden-Thompson inequality\footnote{The Golden-Thompson inequality states that for any two matrices ${\mathbf{A}},{\mathbf{B}}\in\mathbb{S}^n$ it holds that $\textrm{Tr}(\exp({\mathbf{A}}+{\mathbf{B}})) \leq \textrm{Tr}(\exp({\mathbf{A}})\exp({\mathbf{B}}))$.}, and (b) follows from the von Neumann inequality.

Similarly to the proof of \eqref{ineq:StochasticInProofThirdTerm} and using H\"{o}lder's inequality, we get
\begin{align} \label{ineq:boundedDifferences3}
\mathbb{E}[-\log(b_{t-1})\vert{\mathbf{X}}_{t-1}] & \le \mathbb{E}[\eta_{t-1}\langle {\mathbf{X}}_{t-1},\widehat{\nabla}_{t-1}\rangle\vert{\mathbf{X}}_{t-1}] \nonumber
\\ & \le \mathbb{E}[\eta_{t-1}\Vert{\mathbf{X}}_{t-1}\Vert_*\Vert\widehat{\nabla}_{t-1}\Vert_2\vert{\mathbf{X}}_{t-1}] \le \eta_{t-1}G.
\end{align}

Plugging \eqref{ineq:boundedDifferences1},\eqref{ineq:boundedDifferences2} and \eqref{ineq:boundedDifferences3} into \eqref{eq:boundedDifferences} we have
\begin{align*}
W_t-\mathbb{E}[W_{t}\vert W_1,\ldots,W_{t-1}] & \le \eta_{t-1}\Vert\widehat{\nabla}_{t-1}-\nabla f({\mathbf{X}}_{t-1})\Vert_2+2\eta_{t-1}G
\le 4\eta_{t-1}G.
\end{align*}

We will now upper-bound the conditional variance. For any $2\le t\le T$ it holds that
\begin{align} \label{ineq:varBound}
\textbf{Var}[W_{t}\vert W_1,\ldots,W_{t-1}] & = \textbf{Var}[B({\mathbf{X}}^*,{\mathbf{Z}}_{t})\vert {\mathbf{X}}_{t-1}] \nonumber
\\ & \underset{(a)}{=} \textbf{Var}[\textrm{Tr}({\mathbf{X}}^*(\log({\mathbf{X}}_{t-1})-\eta_{t-1}\widehat{\nabla}_{t-1}))-\log(b_{t-1})\vert {\mathbf{X}}_{t-1}] \nonumber
\\ & \underset{(b)}{\le} 2\textbf{Var}[\eta_{t-1}\textrm{Tr}({\mathbf{X}}^*\widehat{\nabla}_{t-1})\vert {\mathbf{X}}_{t-1}]+2\textbf{Var}[\log(b_{t-1})\vert {\mathbf{X}}_{t-1}],
\end{align}
where (a) follows from \eqref{eq:traceYt/bt}, and (b) follows since $\textbf{Var}(X+Y)\le2\textbf{Var}(X)+2\textbf{Var}(Y)$.

\begin{align} \label{ineq:varBound1}
\textbf{Var}[\eta_{t-1}\textrm{Tr}({\mathbf{X}}^*\widehat{\nabla}_{t-1})\vert {\mathbf{X}}_{t-1}] & = \textbf{Var}[\eta_{t-1}\textrm{Tr}({\mathbf{X}}^*\widehat{\nabla}_{t-1})-\eta_{t-1}\textrm{Tr}({\mathbf{X}}^*\nabla f({\mathbf{X}}_{t-1}))\vert {\mathbf{X}}_{t-1}] \nonumber
\\ & \le \eta_{t-1}^2\mathbb{E}\left[\left(\textrm{Tr}({\mathbf{X}}^*(\widehat{\nabla}_{t-1}-\nabla f({\mathbf{X}}_{t-1}))\right)^2\Big\vert {\mathbf{X}}_{t-1}\right] \nonumber
\\ & \underset{(a)}{\le} \eta_{t-1}^2\mathbb{E}\left[\left(\Vert{\mathbf{X}}^*\Vert_*\Vert\widehat{\nabla}_{t-1}-\nabla f({\mathbf{X}}_{t-1})\Vert_2\right)^2\Big\vert {\mathbf{X}}_{t-1}\right] \nonumber
\\ & \le \eta_{t-1}^2\sigma^2/L,
\end{align}
where (a) follows from H\"{o}lder's inequality.

\begin{align} \label{ineq:varBound2}
& \textbf{Var}[\log(b_{t-1})\vert{\mathbf{X}}_{t-1}] = \textbf{Var}\left[\log\left(\textrm{Tr}(\exp(\log({\mathbf{X}}_{t-1})-\eta_{t-1}\widehat{\nabla}_{t-1}))\right)\Big\vert{\mathbf{X}}_{t-1}\right] \nonumber
\\ & = \textbf{Var}\left[\log\left(\textrm{Tr}(\exp(\log({\mathbf{X}}_{t-1})-\eta_{t-1}\widehat{\nabla}_{t-1}))\right)+\log\left(\exp(\eta_{t-1}\lambda_n(\nabla f({\mathbf{X}}_{t-1})))\right)\Big\vert{\mathbf{X}}_{t-1}\right] \nonumber
\\ & = \textbf{Var}\left[\log\left(\textrm{Tr}(\exp(\log({\mathbf{X}}_{t-1})-\eta_{t-1}\widehat{\nabla}_{t-1}))\exp(\eta_{t-1}\lambda_n(\nabla f({\mathbf{X}}_{t-1})))\right)\Big\vert{\mathbf{X}}_{t-1}\right] \nonumber
\\ & \le \mathbb{E}\left[\left(\log\left(\textrm{Tr}(\exp(\log({\mathbf{X}}_{t-1})-\eta_{t-1}\widehat{\nabla}_{t-1}))\exp(\eta_{t-1}\lambda_n(\nabla f({\mathbf{X}}_{t-1})))\right)\right)^2\Big\vert{\mathbf{X}}_{t-1}\right] \nonumber
\\ & \underset{(a)}{\le} \mathbb{E}\left[\left(\log\left(\textrm{Tr}(\exp(\log({\mathbf{X}}_{t-1}))\exp(-\eta_{t-1}\widehat{\nabla}_{t-1}))\exp(\eta_{t-1}\lambda_n(\nabla f({\mathbf{X}}_{t-1})))\right)\right)^2\Big\vert{\mathbf{X}}_{t-1}\right] \nonumber
\\ & =  \mathbb{E}\left[\left(\log\left(\textrm{Tr}({\mathbf{X}}_{t-1}\exp(-\eta_{t-1}\widehat{\nabla}_{t-1}))\exp(\eta_{t-1}\lambda_n(\nabla f({\mathbf{X}}_{t-1})))\right)\right)^2\Big\vert{\mathbf{X}}_{t-1}\right] \nonumber
\\ & \underset{(b)}{\le} \mathbb{E}\left[\left(\log\left(\sum_{i=1}^{n}\lambda_i({\mathbf{X}}_{t-1})\lambda_i(\exp(-\eta_{t-1}\widehat{\nabla}_{t-1}))\exp(\eta_{t-1}\lambda_n(\nabla f({\mathbf{X}}_{t-1})))\right)\right)^2\Bigg\vert{\mathbf{X}}_{t-1}\right] \nonumber
\\ & \underset{(c)}{\le} \mathbb{E}\left[\left(\log\left(\sum_{i=1}^{n}\lambda_i({\mathbf{X}}_{t-1})\exp(\eta_{t-1}\lambda_i(\nabla f({\mathbf{X}}_{t-1})-\widehat{\nabla}_{t-1}))\right)\right)^2\Bigg\vert{\mathbf{X}}_{t-1}\right] \nonumber
\\ & \le \mathbb{E}\left[\left(\log\left(\sum_{i=1}^{n}\lambda_i({\mathbf{X}}_{t-1})\exp(\eta_{t-1}\Vert\widehat{\nabla}_{t-1}-\nabla f({\mathbf{X}}_{t-1})\Vert_2)\right)\right)^2\Bigg\vert{\mathbf{X}}_{t-1}\right] \nonumber
\\ & = \eta_{t-1}^2\mathbb{E}\left[\Vert\widehat{\nabla}_{t-1}-\nabla f({\mathbf{X}}_{t-1})\Vert_2^2\big\vert{\mathbf{X}}_{t-1}\right]
\le \eta_{t-1}^2\sigma^2/L,
\end{align}
where (a) follows from the Golden–Thompson inequality, (b) follows from the von Neumann inequality, and (c) follows from Weyl's inequality.

Plugging \eqref{ineq:varBound1} and \eqref{ineq:varBound2} into \eqref{ineq:varBound}, we get
\begin{align*}
\textbf{Var}[W_{t}\vert W_1,\ldots,W_{t-1}] & \le 4\eta_{t-1}^2\sigma^2/L.
\end{align*}

Now, using a standard concentration argument for sub-martingales (see Theorem 7.3 in \cite{concentrationBound}, which we apply with parameters $\sigma_i=4\eta_{i-1}^2\sigma^2/L$, $\phi_i=0$, $a_i=0$, and $M=\max_{t}\lbrace4\eta_{t-1}G\rbrace=4\eta_{1}G$), we have that for any $\Delta>0$ and $t\in[T]$,
\begin{align*}
\Pr(W_t\ge W_1+\Delta) & \le \exp\left(\frac{-\Delta^2}{\sum_{i=1}^{t-1}8\eta_{i}^2\sigma^2/L+8/3\eta_1G\Delta}\right)
\\ & \le \exp\left(\frac{-\Delta^2}{\sum_{i=1}^{T-1}8\eta_{i}^2\sigma^2/L+8/3\eta_1G\Delta}\right).
\end{align*}

Let $\Delta=\sqrt{2\sum_{i=1}^{T-1}8\eta_{i}^2\sigma^2/L}\sqrt{\log\left(\frac{1}{p'}\right)}$, and suppose that it satisfies $\sum_{i=1}^{T-1}8\eta_{i}^2\sigma^2/L \ge 8/3\eta_1G\Delta$, which in turn is equivalent to the condition
\begin{align*}
\sqrt{\sum_{i=1}^{T-1}\eta_{i}^2\sigma^2/L}\ge \frac{1}{3}\eta_1G\sqrt{2\log\left(\frac{1}{p'}\right)}.
\end{align*}

Then, it holds that $\Pr(W_t\ge W_1+\Delta)\le p'$, which equivalently implies with probability at least $1-p'$
\begin{align*}
B({\mathbf{X}}^*,{\mathbf{Z}}_t) \le B({\mathbf{X}}^*,{\mathbf{Z}}_1)+\frac{1}{2}G^2\sum_{i=1}^{t-1}\eta_i^2+\sum_{i=1}^{t-1}(B({\mathbf{X}}^*,{\mathbf{X}}_{i})-B({\mathbf{X}}^*,{\mathbf{Z}}_{i}))+\Delta.
\end{align*}

Plugging in $B({\mathbf{X}}^*,{\mathbf{Z}}_t)=B({\mathbf{X}}^*,{\mathbf{X}}_t)-(B({\mathbf{X}}^*,{\mathbf{X}}_t)-B({\mathbf{X}}^*,{\mathbf{Z}}_t))$ and recalling that ${\mathbf{Z}}_1={\mathbf{X}}_1$, we have that with probability at least $1-p'$ it holds that
\begin{align*}
B({\mathbf{X}}^*,{\mathbf{X}}_t) \le B({\mathbf{X}}^*,{\mathbf{X}}_1)+\frac{1}{2}G^2\sum_{i=1}^{t-1}\eta_i^2+\sum_{i=1}^{t}(B({\mathbf{X}}^*,{\mathbf{X}}_{i})-B({\mathbf{X}}^*,{\mathbf{Z}}_{i}))+\Delta.
\end{align*}

The Lemma now follows from setting $p' = p/T$ and using the union-bound.

\end{proof}

We can now finally state and prove our local convergence result for Algorithm \ref{alg:LRMD} with stochastic updates. We prove the local convergence under the assumption that the SVD rank parameter $r$ in Algorithm \ref{alg:LRMD} is set to exactly the rank of the corresponding optimal solution $r^*$ (similarly to our first theorem in the deterministic setting - Theorem \ref{thm:smoothCaseDeminish}). The possibility of extending this result to the more general case $r\geq r^*$ (as done in our second theorem for the deterministic setting - Theorem \ref{thm:smoothCaseGen}) remains an open question.

\begin{theorem} \label{thm:StochasticPutTogether}
Let ${\mathbf{X}}^*\in\mathcal{X}^*$ be an optimal solution such that $\textrm{rank}({\mathbf{X}}^*)=r^*$ and suppose ${\mathbf{X}}^*$ satisfies Assumption \ref{ass:strictcomp} with some parameter $\delta >0$. Consider running Algorithm \ref{alg:LRMD} for $T$ iterations with SVD rank parameter $r=r^*$, a fixed minibatch-size $L\ge\max\left\lbrace\frac{16\sigma^2}{G^2R_0^2}\log\left(8T\right),\frac{128G^2}{\delta^2}\log(nT)\right\rbrace$ and a fixed step-size  $\eta=\frac{R_0}{2G\sqrt{T}}$, where
\begin{align*}
R_0 :=  \frac{1}{8}\left[2\beta+\left(1+\frac{2\sqrt{2r^*}}{\lambda_{r^*}({\mathbf{X}}^*)}\right)G\right]^{-1}\delta.
\end{align*}
Suppose further that for all $t\geq 0$: $\varepsilon_t=\frac{9}{32}\frac{\tilde{\varepsilon}_0}{(t+c+1)^2}$ and $c\ge\frac{\delta\sqrt{\tilde{\varepsilon}_0}}{16G\sqrt{T}}$ for some  $\tilde{\varepsilon}_0 \leq R_0^2$. Finally, assume the initialization matrix ${\mathbf{X}}_0$ satisfies $\textrm{rank}({\mathbf{X}}_0) = r^*$, and the parameters ${\mathbf{X}}_0,{\mathbf{X}}^*,\tilde{\varepsilon}_0,R_0$ satisfy the initialization condition in Eq. \eqref{eq:init:equalRank} (substituting ${\mathbf{X}}={\mathbf{X}}_0$, $\varepsilon=\tilde{\varepsilon}_0$, $R=R_0$).
Then, for any $T$ sufficiently large and $\bar{{\mathbf{X}}}\sim\textrm{Uni}\{1,\dots,T\}$ it holds with probability at least $1/2$ that
\[f(\bar{{\mathbf{X}}})-f({\mathbf{X}}^*)=\frac{12GR_0}{\sqrt{T}}.\]
\end{theorem}

\begin{proof}

For all $t\in[T-1]$ it holds that
\[ \Vert  (\nabla f({\mathbf{X}}_t)-\widehat{\nabla}_t)^2\Vert \le 2\Vert  \nabla_t\Vert^2+2\Vert\widehat{\nabla}_t\Vert^2 \le 4G^2.\]
Thus, using a standard Hoeffding concentration argument (see for instance \cite{HoeffdingBound}), we have that with a batch-size of $L\ge\frac{128G^2}{\delta^2}\log(nT)$,
\begin{align*}
\Pr\left(\Vert\nabla f({\mathbf{X}}_t)-\widehat{\nabla}_t\Vert\ge\frac{\delta}{4}\right) & \le 2n\cdot\exp\left({-\frac{\delta^2}{16}}\Bigg/{\frac{8G^2}{\frac{128G^2}{\delta^2}\log(nT)}}\right)
\\ & = 2n\cdot\exp\left(-\log(nT)\right)
= \frac{2}{T}.
\end{align*}
Therefore, for a large enough $T$, with probability at least $9/10$, it holds for all $t\in[T-1]$ that
\begin{align} \label{xsiBound}
\Vert \widehat{\nabla}_t-\nabla f({\mathbf{X}}_t)\Vert_2\le\frac{\delta}{4}.
\end{align}

With a constant step-size $\eta=\frac{R_0}{2G\sqrt{T}}$, and constant batch-size $L\ge\frac{16\sigma^2}{G^2 R_0^2}\log\left(8T\right)$, condition \eqref{ineq:conditionOfConcentrationBoundLemma} holds with $p=1/8$ for any sufficiently large $T$. Thus, invoking Lemma \ref{lemma:concentrationBound}, we have with probability at least $7/8$ it holds that for all $t\in[T]$,
\begin{align} \label{ineq:StochasticRecursionBound}
& B({\mathbf{X}}^*,{\mathbf{X}}_t) \nonumber \\ & \le B({\mathbf{X}}^*,{\mathbf{X}}_1)+\frac{1}{2}G^2\sum_{i=1}^{t-1}\eta_i^2+\sum_{i=1}^{t}(B({\mathbf{X}}^*,{\mathbf{X}}_{i})-B({\mathbf{X}}^*,{\mathbf{Z}}_{i}))+\sqrt{16\sigma^2/L\sum_{i=1}^{T-1}\eta_{i}^2}\sqrt{\log\left(\frac{T}{p}\right)} \nonumber
\\ & \le \frac{13}{8}B({\mathbf{X}}^*,{\mathbf{X}}_1)+\sum_{i=1}^{t}(B({\mathbf{X}}^*,{\mathbf{X}}_{i})-B({\mathbf{X}}^*,{\mathbf{Z}}_{i})).
\end{align}

In addition, note that by our choice of the sequence $\lbrace \varepsilon_t\rbrace_{t\ge0}$, we have that the sequence $\left\lbrace-\frac{1}{\eta}\log\left(\frac{\varepsilon_{t-1}}{\varepsilon_{t}}\right)\right\rbrace_{t\ge1}$ is monotone non-decreasing, and thus, for all $t\ge1$:
\begin{align*}
-\frac{1}{\eta}\log\left(\frac{\varepsilon_{t-1}}{\varepsilon_{t}}\right) & \ge -\frac{1}{\eta}\log\left(\frac{\varepsilon_{0}}{\varepsilon_{1}}\right) = -\frac{4G\sqrt{T}}{R_0}\log\left(\frac{c+2}{c+1}\right)
\\ & \ge -\frac{4G\sqrt{T}}{R_0}\log\left(1+\frac{1}{c}\right) \ge -\frac{4G\sqrt{T}}{R_0}\frac{1}{c} \ge -\frac{\delta}{4},
\end{align*}
where the last inequality follows from plugging in our choice for $c$.

Thus, we have that for all $t\geq 1$:
\begin{align}\label{eq:thm:smoothConv:8}
&\frac{1}{\sqrt{2}}\left[2\beta+\left(1+\frac{2\sqrt{2r^*}}{\lambda_{r^*}({\mathbf{X}}^*)}\right)G\right]^{-1}\frac{\delta}{4} \nonumber
\\ & \le \frac{1}{\sqrt{2}}\left[2\beta+\left(1+\frac{2\sqrt{2r^*}}{\lambda_{r^*}({\mathbf{X}}^*)}\right)G\right]^{-1}\left(\frac{\delta}{2}-\frac{1}{\eta_t}\log\left(\frac{\varepsilon_{t-1}}{\varepsilon_{t}}\right)\right).
\end{align}

Thus, in order to invoke Lemma \ref{lemma:StochasticRadius} for all $t\geq 1$ it suffices to prove that for all $t\geq 1$: i. $\lambda_{r^*+1}({\mathbf{X}}_t) \leq \frac{\varepsilon_{t-1}}{n-r}$ and ii. $\sqrt{B({\mathbf{X}}^*,{\mathbf{X}}_t)} \leq \textrm{LHS of } \eqref{eq:thm:smoothConv:8}$.

The requirement $\lambda_{r^*+1}({\mathbf{X}}_t) \leq \frac{\varepsilon_{t-1}}{n-r}$ holds trivially by the design of the algorithm and since the SVD parameter satisfies $r=r^*$.

We now prove by induction that indeed for all $t\geq 1$, it holds with probability at least $7/8$ that $\sqrt{B({\mathbf{X}}^*,{\mathbf{X}}_t)} \leq \textrm{LHS of } \eqref{eq:thm:smoothConv:8}$.
The base case holds due to the initialization choice.
Now, if the assumption holds for all $i\in[t-1]$ then invoking Lemma \ref{lemma:StochasticRadius} it holds that $B({\mathbf{X}}^*,{\mathbf{X}}_{i})-B({\mathbf{X}}^*,{\mathbf{Z}}_{i})\le 2\varepsilon_i$ for all $i\in[t]$. From \eqref{ineq:StochasticRecursionBound}, this implies that with probability at least $7/8$,
\begin{align*}
B({\mathbf{X}}^*,{\mathbf{X}}_t) & \le \frac{13}{8}B({\mathbf{X}}^*,{\mathbf{X}}_1)+\sum_{i=1}^{t}2\varepsilon_i \nonumber
=\frac{13}{8}B({\mathbf{X}}^*,{\mathbf{X}}_1)+\frac{9\tilde{\varepsilon}_0}{16}\sum_{i=1}^{t}\frac{1}{(i+c+1)^2}
\\ & \le \frac{13}{8}B({\mathbf{X}}^*,{\mathbf{X}}_1)+\frac{9\tilde{\varepsilon}_0}{16}\sum_{i=1}^{\infty}\frac{1}{(i+1)^2} \le 2R_0^2.
\end{align*}

Thus, we obtain that with probability at least $7/8$ for all $t\ge1$:
\begin{align} \label{BtBound}
\sqrt{B({\mathbf{X}}^*,{\mathbf{X}}_t)} & \le \sqrt{2}R_0 = \frac{1}{\sqrt{2}}\left[2\beta+\left(1+\frac{2\sqrt{2r^*}}{\lambda_{r^*}({\mathbf{X}}^*)}\right)G\right]^{-1}\frac{\delta}{4},
\end{align}
and the induction holds.

By combining \eqref{xsiBound}, \eqref{BtBound}, and invoking Lemma \ref{lemma:StochasticRadius} we obtain with probability at least $3/4$ that for all $t\ge1$,
\begin{align} \label{ineq:boundOnStochasticEpsInProof}
B({\mathbf{X}}^*,{\mathbf{X}}_{t})-B({\mathbf{X}}^*,{\mathbf{Z}}_{t})\le2\varepsilon_t.
\end{align}

From Theorem \ref{thm:StochasticConvergence}, it holds that for any $\eta>0$, after $T$ iterations,
\begin{align*}
\mathbb{E}\left[f(\bar{{\mathbf{X}}})\right]-f({\mathbf{X}}^*) & \le \frac{B({\mathbf{X}}^*,{\mathbf{X}}_1)+\frac{G^2}{2}T\eta^2+\sum_{t=1}^T \mathbb{E}[B({\mathbf{X}}^*,{\mathbf{X}}_{t})-B({\mathbf{X}}^*,{\mathbf{Z}}_{t})]}{T\eta}.
\end{align*}
Then, from Markov's inequality, it holds with probability at least $3/4$ that
\begin{align} \label{ineq:stochasticBound}
f(\bar{{\mathbf{X}}})-f({\mathbf{X}}^*) & \le 4\frac{B({\mathbf{X}}^*,{\mathbf{X}}_1)+\frac{G^2}{2}T\eta^2+\sum_{t=1}^T \mathbb{E}[B({\mathbf{X}}^*,{\mathbf{X}}_{t})-B({\mathbf{X}}^*,{\mathbf{Z}}_{t})]}{T\eta}.
\end{align}

In addition, note that
\begin{align} \label{ineq:boundEpsStochastic}
\sum_{t=1}^T 2\varepsilon_{t} & = \frac{9\tilde{\varepsilon}_0}{16}\sum_{t=1}^T \frac{1}{(t+c+1)^2} \le \frac{9\tilde{\varepsilon}_0}{16}\sum_{t=1}^{\infty} \frac{1}{(t+1)^2} \le \frac{3}{8}R_0^2.
\end{align}

Thus, by plugging \eqref{ineq:boundOnStochasticEpsInProof}, \eqref{ineq:boundEpsStochastic}, and our choice for $\eta$ into \eqref{ineq:stochasticBound}, we have with probability at least $1/2$ that
\begin{align*}
f(\bar{{\mathbf{X}}})-f({\mathbf{X}}^*) & \le \frac{12GR_0}{\sqrt{T}}.
\end{align*}

\end{proof}

The following corollary states the overall sample complexity of our method for obtaining $\varepsilon$ approximation error from ``warm-start'' initialization, which, up to logarithmic factors, is optimal in $\varepsilon$.

\begin{corollary}
Under the assumptions of Theorem \ref{thm:StochasticPutTogether}, the overall sample complexity to achieve $f(\bar{{\mathbf{X}}})-f({\mathbf{X}}^*)\le\varepsilon$ with probability at least $1/2$, when initializing from a ``warm-start", is upper-bounded by
\[ \tilde{O}\left(\frac{1}{\varepsilon^2}\max\left\lbrace \sigma^2,G^2\lambda^2_{r^*}({\mathbf{X}}^*)\frac{1}{r^*}\right\rbrace\right). \]
\end{corollary}

\begin{proof}
We have from Theorem \ref{thm:StochasticPutTogether}, that with probability $1/2$ in order to achieve $f(\bar{{\mathbf{X}}})-f({\mathbf{X}}^*)\le\varepsilon$, we need to run
\begin{align*}
\tilde{O}\left(\frac{G^2R_0^2}{\varepsilon^2}\right)
\end{align*}
iterations.

The sample complexity is given by the mini-batch size times the number of iterations it takes to reach an $\varepsilon$ error.
Therefore, it is bounded by
\begin{align*}
& \tilde{O}\left(\frac{G^2R_0^2}{\varepsilon^2}\cdot
\max\left\lbrace\frac{\sigma^2}{G^2 R_0^2},\frac{G^2}{\delta^2}\right\rbrace\right)
= \tilde{O}\left(\frac{1}{\varepsilon^2}\max\left\lbrace \sigma^2,\frac{G^4R_0^2}{\delta^2}\right\rbrace\right)
\\ & = \tilde{O}\left(\frac{1}{\varepsilon^2}\max\left\lbrace \sigma^2,G^4\left(\frac{\lambda_{r^*}({\mathbf{X}}^*)}{(\beta+G)\lambda_{r^*}({\mathbf{X}}^*)+\sqrt{r^*}G}\right)^2\right\rbrace\right)
\\ & = \tilde{O}\left(\frac{1}{\varepsilon^2}\max\left\lbrace \sigma^2,G^4\left(\frac{\lambda_{r^*}({\mathbf{X}}^*)}{\sqrt{r^*}G}\right)^2\right\rbrace\right)
 = \tilde{O}\left(\frac{1}{\varepsilon^2}\max\left\lbrace \sigma^2,G^2\lambda^2_{r^*}({\mathbf{X}}^*)\frac{1}{r^*}\right\rbrace\right),
\end{align*}
where the second to last equation follows since $\lambda_{r^*}({\mathbf{X}}^*)\le\frac{1}{r^*}\le\sqrt{r^*}$.

\end{proof}

\section{Experiments} \label{section:experiments}
In this section we present preliminary empirical evidence in support of our theoretical findings.

We consider the problem of recovering a low-rank matrix from quadratic measurements. Concretely, we consider the following optimization problem:
\begin{align} \label{eq:recoveryProb}
\min_{{\mathbf{X}}\succeq0\ \textrm{Tr}({\mathbf{X}})=\tau}\left\lbrace f({\mathbf{X}}):= \frac{1}{2}\sum_{i=1}^m \left(\textrm{Tr}({\mathbf{a}}_i^{\top}{\mathbf{X}}{\mathbf{b}}_i)-{\mathbf{y}}_i\right)^2\right\rbrace.
\end{align}

Throughout this section we focus on the deterministic setting in which full gradients of the objective  function in \eqref{eq:recoveryProb} are available.

We let ${\mathbf{M}}=(\sqrt{n}{\mathbf{V}})(\sqrt{n}{\mathbf{V}})^{\top}$ be the ground-truth low-rank matrix,  where ${\mathbf{V}}\in\mathbb{R}^{n\times r}$ is generated by taking a random matrix with standard Gaussian entries and then normalizing it to have unit Frobenius norm, i.e., $\Vert{\mathbf{V}}\Vert_F=1$.
We generate $m$ pairs of random uniformly-distributed unit vectors $\lbrace({\mathbf{a}}_i,{\mathbf{b}}_i)\rbrace_{i=1}^m\subset\mathbb{R}^n\times\mathbb{R}^n$. ${\mathbf{y}}_0\in\mathbb{R}^m$, the vector of quadratic measurements of ${\mathbf{M}}$, is given by ${\mathbf{y}}_0(i)={\mathbf{a}}_i^{\top}{\mathbf{M}}{\mathbf{b}}_i$.
We then add noise to produce the noisy vector ${\mathbf{y}}={\mathbf{y}}_0+{\mathbf{n}}$, where ${\mathbf{n}}=\kappa\Vert{\mathbf{y}}_0\Vert_2{\mathbf{v}}$ for a random unit vector ${\mathbf{v}}\in\mathbb{R}^m$ and $\kappa\in\mathbb{R}$. In all of our experiments we set $\kappa=1/2$.

In all of our experiments, we scale our final estimate for the ground-truth matrix ${\mathbf{M}}$ to be $\frac{\textrm{Tr}({\mathbf{M}})}{\tau}{\mathbf{X}}^*$, where $\tau$ is the trace bound in \eqref{eq:recoveryProb} and ${\mathbf{X}}^*$ is obtained by solving \eqref{eq:recoveryProb},  in order for it to have the same trace as ${\mathbf{M}}$.
We measure the relative recovery error by $\left\Vert\frac{\textrm{Tr}({\mathbf{M}})}{\tau}{\mathbf{X}}^*-{\mathbf{M}}\right\Vert_F^2\Big/\left\Vert{\mathbf{M}}\right\Vert_F^2$. Similarly, the initial relative error of the matrix ${\mathbf{X}}_0$ used to initialize our algorithm is given by $\left\Vert\frac{\textrm{Tr}({\mathbf{M}})}{\tau}{\mathbf{X}}_0-{\mathbf{M}}\right\Vert_F^2\Big/\left\Vert{\mathbf{M}}\right\Vert_F^2$.
The signal-to-noise ratio is given by $\Vert{\mathbf{y}}_0\Vert^2 /\Vert{\mathbf{n}}\Vert^2=4$ (under our fixed choice of $\kappa=1/2$).

Since we cannot solve Problem \eqref{eq:recoveryProb} exactly, to verify the near-optimality of the found solution ${\mathbf{X}}^*$, we compute for it the dual-gap  which is given by
\begin{align*}
\max_{{\mathbf{Z}}\succeq0\ \textrm{Tr}({\mathbf{Z}})=\tau}{\langle{\mathbf{X}}^*-{\mathbf{Z}},\nabla f({\mathbf{X}}^*)\rangle} = \textrm{Tr}(({\mathbf{X}}^*-\tau{\mathbf{v}}_{n}{\mathbf{v}}_{n}^{\top})\nabla f({\mathbf{X}}^*)),
\end{align*}
where ${\mathbf{v}}_n$ is the eigen-vector corresponding to the smallest eigenvalue of $\nabla f({\mathbf{X}}^*)$. Note that due to convexity of the objective, the dual gap is always an upper-bound on the approximation error w.r.t. function value.

For all experiments we set the initialization matrix to
\begin{align}\label{eq:exp:init}
{\mathbf{X}}_0 = {\mathbf{V}}_r\textrm{diag}\left(\Pi_{\Delta_{\tau,r}}[\textrm{diag}(-\Lambda_r)]\right){\mathbf{V}}_r^{\top},
\end{align}
where $\Delta_{\tau,r}=\lbrace{\mathbf{z}}\in\mathbb{R}^r~|~{\mathbf{z}}\ge0,\ \sum_{i=1}^r {\mathbf{z}}_i=\tau\rbrace$ is the $\tau$-simplex in $\mathbb{R}^r$ and $\Pi_{\Delta_{\tau,r}}[\cdot]$ denotes the Euclidean projection over it, ${\mathbf{V}}_r\Lambda_r{\mathbf{V}}_r$ is the rank-$r$ eigen-decomposition of $-\nabla f(\tau{\mathbf{U}}\U^{\top})$, and ${\mathbf{U}}\in\mathbb{R}^{n\times r}$ is produced by taking a random matrix with standard Gaussian entries and normalizing it to have a unit Frobenius norm.

In all experiments we set the sequence $\{\varepsilon_t\}_{t\geq 0}$ in Algorithm \ref{alg:LRMD} to $\varepsilon_t=\frac{3}{5}\frac{1}{(t+c+1)^2}$ where $c=10$. We set the number of measurements in the objective \eqref{eq:recoveryProb} to $m=20nr$, the number of iterations to $T=200$, and the smoothness parameter to $\beta=0.4\sqrt{rn}$. When the ground-truth matrix is  rank-$1$ or rank-$5$ we set $\tau=0.5\cdot\textrm{Tr}({\mathbf{M}})$, and for rank-$20$ we set $\tau=0.65\cdot\textrm{Tr}({\mathbf{M}})$. We note that we choose the trace bound $\tau$ to be strictly smaller than $\textrm{Tr}({\mathbf{M}})$ to prevent the optimal solution from fitting some of the noise and resulting in a higher rank matrix. For every set of parameters ($r,n$) we take the averages of $20$ i.i.d runs.

\subsection{Rank of ground-truth matrix is known}

In our first line of experiments we assume that $\textrm{rank}({\mathbf{M}})$ is known and we set the SVD rank parameter in Algorithm \ref{alg:LRMD} to $r=\textrm{rank}({\mathbf{M}})$. This parameter is also used to set the initialization matrix ${\mathbf{X}}_0$.

We record several quantities of interest in Table \ref{table:results}. Importantly, Table \ref{table:results} indicates that i. for our recovery setting the strict complementarity condition (Assumption \ref{ass:strictcomp}) indeed holds and seems dimension-independent, and ii. in all cases, the convergence certification condition
\begin{align} \label{convergenceGuaranteeExperiment}
\log\left(\frac{(n-r)\lambda_{r+1}({\mathbf{Y}}_t)}{\varepsilon_t b_t}\right)\le 2\varepsilon_t, \quad b_t =\sum_{i=1}^n\lambda_i({\mathbf{Y}}_t)
\end{align}

(see Section \ref{sec:certificate}) holds from the early stages of the run (right from the first iteration in most cases) and throughout  all following iterations, implying the correct convergence of our low-rank MEG method, up to negligible error.

\begin{table}\renewcommand{\arraystretch}{1.3}
\begin{subtable}{\textwidth}
        {\small
        \begin{center}
 	    \begin{tabular}{| c | c | c | c | c | c |} \hline
   		dimension & avg. initialization & avg. 	recovery &  avg. gap in & avg. dual & avg. first iter. \\
    ($n$) &  error & error & $\nabla f({\mathbf{X}}^*)$ & gap & \eqref{convergenceGuaranteeExperiment} holds \\ \hline
   $100$ & $0.1593$ & $0.0439$ & $5.3665$ & $0.0065$ & $1$ \\ \hline
   $200$ & $0.1674$ & $0.0471$ & $5.2046$ & $0.0129$ & $1$ \\ \hline
   $400$ & $0.1732$ & $0.0487$ & $5.3049$ & $0.0264$ & $1$ \\ \hline
   $600$ & $0.1702$ & $0.0489$ & $5.2400$ & $0.0395$ & $1$ \\ \hline
  \end{tabular}
  \caption*{$\textrm{rank}({\mathbf{M}})=1$}\label{table:numericalResults1}
  \end{center}
  }
  \vskip -0.2in
\end{subtable}

\begin{subtable}{\textwidth}
        {\small
        \begin{center}
  \begin{tabular}{| c | c | c | c | c | c |} \hline
    dimension & avg. initialization & avg. recovery &  avg. gap in & avg. dual & avg. first iter. \\
    ($n$) &  error & error & $\nabla f({\mathbf{X}}^*)$ & gap & \eqref{convergenceGuaranteeExperiment} holds \\ \hline
   $100$ & $0.5763$ & $0.0780$ & $5.8528$ & $0.0079$ & $1$ \\ \hline
   $200$ & $0.2072$ & $0.0561$ & $5.6697$ & $0.0136$ & $1$ \\ \hline
   $400$ & $0.1307$ & $0.0451$ & $5.5942$ & $0.0271$ & $1$ \\ \hline
   $600$ & $0.1101$ & $0.0381$ & $5.5771$ & $0.0404$ & $1$ \\ \hline
  \end{tabular}
  \caption*{$\textrm{rank}({\mathbf{M}})=5$}\label{table:numericalResults5}
  \end{center}
}
  \vskip -0.2in
\end{subtable}

\begin{subtable}{\textwidth}
  {\small
  \begin{center}
  \begin{tabular}{| c | c | c | c | c | c |} \hline
    dimension & avg. initialization & avg. recovery &  avg. gap in  & avg. dual & avg. first iter. \\
    ($n$) &  error & error & $\nabla f({\mathbf{X}}^*)$ & gap & \eqref{convergenceGuaranteeExperiment} holds \\ \hline
   $100$ & $2.2782$ & $0.0638$ & $3.3870$ & $0.0167$ & $5.7$ \\ \hline
   $200$ & $1.2007$ & $0.0437$ & $3.3210$ & $0.0156$ & $8.6$ \\ \hline
   $300$ & $0.6729$ & $0.0365$ & $3.2839$ & $0.0265$ & $11.35$ \\ \hline
  \end{tabular}
  \caption*{$\textrm{rank}({\mathbf{M}})=20$}\label{table:numericalResults20}
  \end{center}
  }
  \vskip -0.2in
\end{subtable}
\caption{Results for the low-rank matrix recovery from quadratic measurements problem using the Low-Rank Matrix Exponentiated Gradient Method. The fourth column is the spectral gap $\lambda_{n-r}(\nabla{}f({\mathbf{X}}^*)) - \lambda_{n}(\nabla{}f({\mathbf{X}}^*))$, and the last column records the iteration of the algorithm at which  the convergence certificate \eqref{convergenceGuaranteeExperiment} begins to take effect.}
\label{table:results}
\end{table}\renewcommand{\arraystretch}{1}

We move on to compare the empirical convergence of the standard MEG algorithm, Algorithm \ref{alg:exactMEG}, and our low rank variant, Algorithm \ref{alg:LRMD}. For every run we initialize both algorithms with the same matrix (as given in Eq. \eqref{eq:exp:init}). In addition to measuring convergence in function value, we also measure the Bregman distance $B({\mathbf{W}}_t,{\mathbf{X}}_t)$ between the iterates $\lbrace {\mathbf{X}}_t\rbrace_{t\ge1}$ produced by our low-rank algorithm and the iterates $\lbrace {\mathbf{W}}_t\rbrace_{t\ge1}$ produced by the standard MEG algorithm. For all graphs we set $n=200$ and $\textrm{rank}({\mathbf{M}})=1,5,20$. The results are shown in Figure \ref{fig:comparison between methods}.

As can be seen, both the standard MEG method and our low rank variant converge very similarly in terms of the function value.  In addition, it can be seen that the Bregman distance $B({\mathbf{W}}_t,{\mathbf{X}}_t)$ decays very quickly, providing additional evidence for the correct convergence of our low rank variant.

\begin{figure}
\centering
  \begin{subfigure}[t]{0.49\textwidth}
    \includegraphics[width=\textwidth]{iteration_n200_r1_zoomIn}

    \label{fig:1}
  \end{subfigure}\hfil
  \begin{subfigure}[t]{0.49\textwidth}
    \includegraphics[width=\textwidth]{breg_n200_r1}

    \label{fig:2}
  \end{subfigure}\hfil
  \caption*{$\textrm{rank}({\mathbf{M}})=1$}
  \label{fig:distIteration200_1}
\medskip
  \begin{subfigure}[t]{0.49\textwidth}
    \includegraphics[width=\textwidth]{iteration_n200_r5_zoomIn}

    \label{fig:3}
  \end{subfigure}
  \begin{subfigure}[t]{0.49\textwidth}
    \includegraphics[width=\textwidth]{breg_n200_r5}

    \label{fig:4}
  \end{subfigure}
  \caption*{$\textrm{rank}({\mathbf{M}})=5$}
  \label{fig:distIteration200_5}
\medskip
  \begin{subfigure}[t]{0.49\textwidth}
    \includegraphics[width=\textwidth]{iteration_n200_r20_zoomIn}

    \label{fig:5}
  \end{subfigure}
  \begin{subfigure}[t]{0.49\textwidth}
    \includegraphics[width=\textwidth]{breg_n200_r20}

    \label{fig:6}
  \end{subfigure}
  \caption*{$\textrm{rank}({\mathbf{M}})=20$}
  \caption{Comparison between the low rank MEG algorithm and the standard MEG method for the low-rank matrix recovery from quadratic measurements problem with $n=200$. On the left we plot the convergence in function values for both variants and on the right we plot the Bregman distance between their iterates.} \label{fig:comparison between methods}
\end{figure}

Since the certificate \eqref{convergenceGuaranteeExperiment} requires computing all eigenvalues of the matrix ${\mathbf{Y}}_t$ (due to the scalar $b_t$) which makes no sense in practice when our goal is to use only low rank SVD computations, as discussed in Section \ref{sec:certificate}, we consider the harder-to-satisfy certificate
\begin{align} \label{convergenceGuaranteeExperiment2}
\log\left(\frac{(n-r)\lambda_{r+1}({\mathbf{Y}}_t)}{\varepsilon_t b^{r+1}_t}\right)\le 2\varepsilon_t, \quad b^{r+1}_t =\sum_{i=1}^{r+1}\lambda_i({\mathbf{Y}}_t).
\end{align}

Recall that when \eqref{convergenceGuaranteeExperiment2} holds, \eqref{convergenceGuaranteeExperiment} also holds, but not the other way around. Also, computing \eqref{convergenceGuaranteeExperiment2} only requires computing the top $r+1$ eigenvalues of ${\mathbf{Y}}_t$, i.e., increasing the SVD rank used in Algorithm \ref{alg:LRMD} by only 1. In Table \ref{table:at_bt} we provide evidence that \eqref{convergenceGuaranteeExperiment2} is indeed a reliable certificate: in all runs it  holds starting from the same iteration for which the tighter yet  computationally more expensive certificate \eqref{convergenceGuaranteeExperiment} holds.
\begin{table*}\renewcommand{\arraystretch}{1.3}
{\small
\begin{center}
  \begin{tabular}{| c | c | c | c |} \hline
      & avg. first iter.  & avg. first iter. \\
   $\textrm{rank}({\mathbf{M}})$  & $\log\left(\frac{(n-r)\lambda_{r+1}({\mathbf{Y}}_t)}{\varepsilon_tb_t}\right)\le 2\varepsilon_t$ holds & $\log\left(\frac{(n-r)\lambda_{r+1}({\mathbf{Y}}_t)}{\varepsilon_tb^{r+1}_t}\right)\le 2\varepsilon_t$ holds\\ \hline
   $1$ & $1$  & $1$\\ \hline
   $5$ & $1$ & $1$\\ \hline
   $20$ & $8.6$ & $8.6$\\ \hline
  \end{tabular}
  \caption{Comparison between the first iteration the convergence guarantee \eqref{convergenceGuaranteeExperiment} holds vs. the first iteration the weaker condition $\log\left(\frac{(n-r)\lambda_{r+1}({\mathbf{Y}}_t)}{\varepsilon_tb^{r+1}_t}\right)\le 2\varepsilon_t$ holds for $n=200$.
  }\label{table:at_bt}
\end{center}
}
\vskip -0.2in
\end{table*}\renewcommand{\arraystretch}{1}

\subsection{Rank of ground-truth matrix is overestimated}
In a second line of experiments we examine the case in which $\textrm{rank}({\mathbf{M}})$ - the rank of the ground-truth matrix is not precisely known, but only overestimated. That is, in all following experiments, both for the SVD computations and for the initialization in Algorithm \ref{alg:LRMD}, we use rank parameter $r > \textrm{rank}({\mathbf{M}})$.

We run all experiments with $n=200$. As can be seen in Table \ref{table:bigR}, the convergence guarantee of \eqref{convergenceGuaranteeExperiment} begins to hold starting from a very early stage.
In Figure \ref{fig:comparison between methods bigR} it can be seen that the results using $r>\textrm{rank}({\mathbf{M}})$ are very similar to the previously examined case $r = \textrm{rank}({\mathbf{M}})$, and that the standard MEG method and our low-rank variant perform very similarly also in this case.
\begin{table*}[!htb]\renewcommand{\arraystretch}{1.3}
{\small
\begin{center}
  \begin{tabular}{| c | c | c | c |} \hline
   $\textrm{rank}({\mathbf{M}})$  & SVD \textrm{rank} & avg. first iteration \eqref{convergenceGuaranteeExperiment} holds\\ \hline
   $1$ & $3$  & $2$\\ \hline
   $5$ & $10$ & $3.85$\\ \hline
   $20$ & $30$ & $9.25$\\ \hline
  \end{tabular}
  \caption{The first iteration the convergence certificate \eqref{convergenceGuaranteeExperiment} holds when using $r>\textrm{rank}({\mathbf{M}})$ for $n=200$.
  }\label{table:bigR}
\end{center}
}
\vskip -0.2in
\end{table*}\renewcommand{\arraystretch}{1}

\begin{figure}
\centering
  \begin{subfigure}[t]{0.49\textwidth}
    \includegraphics[width=\textwidth]{iteration_n200_r1_zoomIn_rsvd3}

    \label{fig:12}
  \end{subfigure}\hfil
  \begin{subfigure}[t]{0.49\textwidth}
    \includegraphics[width=\textwidth]{breg_n200_r1_rsvd3}

    \label{fig:22}
  \end{subfigure}\hfil
  \caption*{$\textrm{rank}({\mathbf{M}})=1$, SVD rank = $3$}
  \label{fig:bigRdistIteration200_1}
\medskip
  \begin{subfigure}[t]{0.49\textwidth}
    \includegraphics[width=\textwidth]{iteration_n200_r5_zoomIn_rsvd10}

    \label{fig:32}
  \end{subfigure}
  \begin{subfigure}[t]{0.49\textwidth}
    \includegraphics[width=\textwidth]{breg_n200_r5_rsvd10}

    \label{fig:42}
  \end{subfigure}
  \caption*{$\textrm{rank}({\mathbf{M}})=5$, SVD rank = $10$}
  \label{fig:bigRdistIteration200_5}
\medskip
  \begin{subfigure}[t]{0.49\textwidth}
    \includegraphics[width=\textwidth]{iteration_n200_r20_zoomIn_rsvd30}

    \label{fig:52}
  \end{subfigure}
  \begin{subfigure}[t]{0.49\textwidth}
    \includegraphics[width=\textwidth]{breg_n200_r20_rsvd30}

    \label{fig:62}
  \end{subfigure}
  \caption*{$\textrm{rank}({\mathbf{M}})=20$, SVD rank = $30$}
  \caption{Comparison between the low rank MEG and the standard MEG for the low-rank matrix recovery from quadratic measurements problem where $n=200$ and $r>\textrm{rank}({\mathbf{M}})$. On the left we plot the convergence in function values for both variants and on the right we plot the Bregman distance between their iterates.} \label{fig:comparison between methods bigR}
\end{figure}

\bibliographystyle{plain}

redpajama_set_name: RedPajamaArXiv