\section{Preliminaries}\label{Sec:Prelimnaries}
Consider the following optimization problem,
\begin{equation}\label{prob:BoundedMapped}
\tag{Constrained-OPT}
\max_{x \in \mathcal{S}}\ \ g(\mathcal{B}(x)) = \max_{v\in \mathcal{B}(\mathcal{S})}\ \ g(v)
\end{equation}
where $\mathcal{B}(\cdot):\mathbb{R}^{m}\rightarrow\mathbb{R}^{d}$ is a linear map, $g(\cdot)$ is a smooth, concave function, and $\mathcal{S}\subseteq \mathbb{R}^{m}$ is a compact, convex set. The variable $v = \mathcal{B}(x)$ is said to be a `projection' of the `lifted' decision variable $x$ of~\eqref{prob:BoundedMapped} and the map $\mathcal{B}(\cdot)$ can be interpreted as linear measurements of the decision variable. In this section, we briefly review the Frank-Wolfe algorithm~\cite{frank1956algorithm}, which is central to our algorithmic approach. We recall the modification done to the steps of Frank-Wolfe to adapt it to problems of type~\eqref{prob:BoundedMapped} by~\citet{sketchy} and give subsequent convergence results. Moreover, we also discuss~\eqref{prob:PrimalSDP}, which is a special case of~\eqref{prob:BoundedMapped}.

\begin{comment}
Any first-order method requires knowledge of the gradient of the objective function and for~\eqref{prob:BoundedMapped}, computing the gradient over the feasible domain $\mathcal{S}$ requires just the knowledge of projected iterate $v_t$. The use of linear map was inspired by~\citet{sketchy} which approximates the objective function value with a linear map and the projected variables keep track of the gradient of the objective function.
\end{comment}


\begin{comment}
\begin{remark}
A mapping $\mathcal{B}(X) = [[\langle B_i , X\rangle]_{i=1}^{d}]: \mathbb{S}^{n} \rightarrow \mathbb{R}^d$ of the variable preserves the following: A feasible point, $v^* = \mathcal{B}(X^*)$, is the optimal solution of 
\begin{equation}\label{prob:SDPBoundedMapped}
\begin{split}
\max\ &g(\mathcal{B}(X))\\
\textup{s.t.}\ &X \in \mathcal{S},
\end{split}
\end{equation}
if and only if $X^*$ is the optimal solution of~\eqref{prob:PrimalSDP}.
\end{remark}
\end{comment}

\subsection{An Approximate Frank-Wolfe Algorithm}
The Frank-Wolfe algorithm~\cite{frank1956algorithm} is an iterative algorithm for convex optimization over a compact, convex feasible region that solves a linear optimization problem at each iterate. The algorithm then, in each iteration, takes a fractional step towards a maximizer of the linear optimization subproblem which can be taken to be an extreme point of the feasible set. Frank-Wolfe is an example of a `projection-free' method since the fractional step ensures that the next iterate is feasible, eliminating the need to project onto the feasible region.

An approximate Frank-Wolfe algorithm (see, e.g.,~\cite{FW}) computes an $\epsilon$-optimal solution of~\eqref{prob:BoundedMapped} by solving the linear subproblem approximately (with some finite additive error) at each iteration.

\begin{definition}
A feasible point $\bar{x}\in \mathcal{S}$ is called $\epsilon$-optimal for the optimization problem $\max_{x\in \mathcal{S}} f(x)$ if $f(\bar{x}) \geq \max_{x\in \mathcal{S}} f(x) -\epsilon$.\end{definition}

For a concave objective function $g$ and a convex compact feasible domain $\mathcal{S}$, Algorithm~\ref{algo:FWapprox} details an approximate Frank-Wolfe algorithm to compute an $\epsilon$-optimal solution to~\eqref{prob:BoundedMapped}.

\begin{algorithm}[ht] 
\caption{Frank-Wolfe Algorithm with Approximate Solution to Subproblem~\eqref{prob:Subproblem}}
 \label{algo:FWapprox}
\SetAlgoLined
\SetKwInOut{Input}{Input}\SetKwInOut{Output}{Output}
\DontPrintSemicolon

\Input{Stopping criteria $\epsilon$, accuracy parameter $\eta$, probability $p$ for subproblem~\eqref{prob:Subproblem}, upper bound $C_g^u$ on the curvature constant}
\Output{$\epsilon$-optimal maximizer $x_t$ of $g(\mathcal{B}(x))$, $v_t = \mathcal{B}(x_t)$}
\BlankLine

 \SetKwFunction{CDD}{LMO}
 \SetKwFunction{UV}{UpdateVariable}
 \SetKwFunction{FWA}{FWApproxSubprob}
 \SetKwProg{Fn}{Function}{:}{\KwRet{$(x_t,v_t)$}}

 \Fn{\FWA}{
 Select initial point $x_0\in \mathcal{S}$ and set $v_0 = \mathcal{B}(x_0)$\;
 $t \leftarrow 0$\;
 $\gamma \leftarrow \frac{2}{t+2}$\;
 $(h_t, q_t) \leftarrow \CDD (\mathcal{B}^*(\nabla g(v_{t})), \frac{1}{2}\eta\gamma C_g^u,p)$\;
 
 \While{$\left\langle q_t - v_t, \nabla g(v_t)\right\rangle > \epsilon$}{
 
 $(x_{t+1},v_{t+1}) \leftarrow \UV (x_t, v_t, h_t, q_t,\gamma)$
 
 $t \leftarrow t+1$\;
 $\gamma \leftarrow \frac{2}{t+2}$\;
 
 $(h_t, q_t) \leftarrow \CDD (\mathcal{B}^*(\nabla g(v_{t})), \frac{1}{2}\eta\gamma C_g^u,p$)\;
 
 }
 }
 
 \SetKwProg{Fn}{Function}{:}{\KwRet{$(h,q)$}}
 \Fn{\CDD{$J, \delta$, $p$}}{
 
 Find $h\in \mathcal{S}$ such that with probability at least $1-\SubProb$, $\left\langle h, J\right\rangle \geq \max_{s\in \mathcal{S}} \left\langle s, J\right\rangle - \delta$ \;
 $q = \mathcal{B}(h)$\;
 
 }
 
 \SetKwProg{Fn}{Function}{:}{\KwRet{$(x,v)$}}
 \Fn{\UV{$x,v,h,q,\gamma$}}{
 $x \leftarrow (1-\gamma)x + \gamma h$\;
 $v \leftarrow (1-\gamma)v + \gamma q$\; 
 
 }
\end{algorithm}

In Algorithm~\ref{algo:FWapprox}, we assume that at each iteration, the linear optimization subproblem \texttt{LMO} is solved by a randomized method that succeeds with probability at least $1-p$. Moreover, we assume the randomness in this subproblem is independent across function calls.
In Algorithm~\ref{algo:FWapprox}, we not only keep track of the `lifted' iterate $x_t\in \mathbb{R}^n$ (decision variable) at each iteration but also the `projected' iterate $v_t\in\mathbb{R}^d$. The knowledge of these projected iterates $v_t$ is enough to compute the gradient of the objective function, and hence to compute the update direction. Thus, for a high dimensional decision variable, we can simply track the projected iterates $v_t$ to track the improvement in objective function value at each iterate. This is a key observation from the work of~\citet{sketchy} that is also crucial here.

\paragraph{Curvature constant.} The curvature constant $C_g$ is a measure of nonlinearity of the objective function $g$ over the feasible set $\mathcal{S}$. It is defined (see, e.g., \cite{FW}) as,
\begin{equation*}
C_g = \sup_{\substack{x,h\in \mathcal{S},\\ y = v + \gamma(h-v),\\ \gamma \in [0,1]}} -\frac{2}{\gamma^2} \left(g(y) - g(v) - \nabla g(v)^T (y-v)\right).
\end{equation*}
The value of the curvature constant provides insight into the deviation of the linearization of the function $g$ from the actual function value. If $g$ is twice-differentiable, then we often use the upper bound, $C_g \leq C_g^u = \lambda_{\textup{max}}(-\nabla^2 g) \textup{diam}(\mathcal{S})^2$, where $\textup{diam}(\mathcal{S})$ is the Euclidean diameter of the set $\mathcal{S}$.

\paragraph{Approximate solution to subproblem in \texttt{LMO}.} The subroutine \texttt{LMO} solves the linear maximization problem
\begin{equation}\label{prob:Subproblem}
\max_{h\in \mathcal{S}}\ \ \left\langle h, \mathcal{B}^*(\nabla g(v))\right\rangle
\end{equation}
approximately at each iteration $t$. The approximate maximizer $h_t$ is computed such that with probability at least $1-\SubProb$, the  additive error in the function value at $h_t$ is at most $\frac{1}{2}\eta\gamma C_g^u$, where $\eta \geq 0$ is a fixed accuracy parameter and $C_g^u$ is an upper bound on the curvature constant of $g$. Initially, this margin of error is high; but with each iteration we solve the subproblem~\eqref{prob:Subproblem} to a higher accuracy. When $\SubProb = 0$, Algorithm~\ref{algo:FWapprox} is equivalent to the Frank-Wolfe algorithm given by~\citet*{FW}. By setting $p>0$, we require the subroutine to adhere to the desired accuracy with some probability of failure.


\begin{comment}
For~\eqref{prob:PrimalSDP}, solving subproblem~\eqref{prob:Subproblem} is equivalent to computing the largest eigenvector of a matrix, i.e.,
\begin{equation*}
    H_t=
    \begin{cases}
      \alpha u_tu_t^T, & \text{if}\ \lambda_t\geq 0 \\
      0, & \textup{otherwise},
    \end{cases}
\end{equation*}
and $(\lambda_t,u_t)$ is the largest (eigenvalue, eigenvector) pair of matrix $\mathcal{B}^*(\nabla g(v_t))$ and $v_t = \mathcal{B}(X_t)$.

Defining $h_t = \mathcal{B}(H_t)$, the update step is written as,
\begin{equation}
v_{t+1} = (1-\gamma)v_{t} + \gamma h_t.
\end{equation}

\paragraph{Computing approximate eigenvectors (solving the subproblem~\eqref{prob:Subproblem}) in the Frank-Wolfe Algorithm.}
For large symmetric matrices,~\citet{money2005algorithm} has developed a technique which can be used to compute the required eigenvectors approximately. It projects the matrix onto a smaller subspace and computes the eigenvectors using Lanczos algorithm. The memory requirement is lesser since computations are done on a smaller matrix.
\end{comment}

\paragraph{Stopping criterion.}
Algorithm~\ref{algo:FWapprox} terminates when
\begin{equation*}
\left\langle q_t - v_t,\nabla g(v_t)\right\rangle \leq \epsilon
\end{equation*}
is satisfied. 
When this condition holds, it follows that $g(v_t) \geq g(v^{\star}) - \epsilon$, i.e., $v_t$ is an $\epsilon$-optimal solution of~\eqref{prob:BoundedMapped} \cite{FW}.

\paragraph{Convergence of the Frank-Wolfe algorithm.}
\begin{theorem} \label{thm:FWapproxConv}
Let $g:  \mathbb{R}^{d} \rightarrow \mathbb{R}$ be a concave and differentiable function and $x^{\star}$ an optimal solution of~\eqref{prob:BoundedMapped}. If $C_g^u$ is an upper bound on the curvature constant of $g$, and $\eta \geq 0$ is the accuracy parameter for subproblem~\eqref{prob:Subproblem}, then $x_t$, the $t$-th iterate of Algorithm~\ref{algo:FWapprox}, satisfies
\begin{equation}
-g(\mathcal{B}(x_{t})) + g(\mathcal{B}(x^{\star})) \leq \frac{2C_g^u(1+\eta)}{t+2}
\end{equation}
with probability $(1-\SubProb)^t \geq 1- t \SubProb$.
\end{theorem}
\begin{proof}
The results follows from \cite[Theorem~1]{FW} when the subproblem~\eqref{prob:Subproblem} is solved to an accuracy of $\frac{1}{2}\eta\gamma C_g^u$ with probability $1-p$.
\end{proof}

Thus, after $t =\frac{2C_g^u(1+\eta)}{\epsilon} - 2$ iterations, the solution $x_t$ satisfies
\begin{equation}\label{eqn:SubOptimalSolPenaltyProb1}
g(\mathcal{B}(x_{t})) \geq g(\mathcal{B}(x^\star)) - \epsilon
\end{equation}
with probability at least $1-tp$.
So, we can now compute an $\epsilon$-optimal solution to~\eqref{prob:BoundedMapped} in $t \sim \mathcal{O}(C_g^u/\epsilon)$ iterations with probability at least $1-tp$.

In Algorithm~\ref{algo:FWapprox}, we solve the subproblem~\eqref{prob:Subproblem} approximately to find an update direction that is an extreme point of the set $\mathcal{S}$. If the following  conditions on the constraint set $\mathcal{S}$ are satisfied, then we potentially generate the sampled representation of the solution to Problem~\eqref{prob:BoundedMapped} in low memory as illustrated in Section~\ref{sec:SDPBoundedTrace}.

\begin{description}
\item[\namedlabel{condition1}{LowMemoryComputations}] The image of the linear map $\mathcal{B}$ is low dimensional and the subproblem~\eqref{prob:Subproblem} can be solved in low memory.
\item[\namedlabel{condition2}{LowMemoryExtremePoints}] The extreme points of the feasible set $\mathcal{S}$ can be represented in low memory.
\end{description}

\subsection{SDP with Bounded Trace Constraint}\label{sec:SDPBoundedTrace}
For most of the paper, our focus is on semidefinite programming problems. As such, we briefly focus on~\eqref{prob:PrimalSDP}, a special SDP that has a trace constrained feasible set and a concave objective function $g$. The extreme points of the feasible set $\mathcal{S}= \{ X\succeq 0: \textup{Tr}(X)\leq \alpha\}$ are rank-1 PSD matrices and can be represented in low memory. Moreover, a solution of subproblem~\eqref{prob:Subproblem} at iterate $t$ is an eigenvector of the matrix $J = \mathcal{B}^*(\nabla g(\mathcal{B}(X_t)))$ corresponding to its largest eigenvalue. To solve this subproblem approximately, we find $w_t$ such that $\alpha w_t^TJw_t\geq \alpha \lambda_{\textup{max}}(J) - \frac{1}{2}\eta\gamma C_g^u$ and then select the update direction
\begin{equation}\label{eqn:HtforBoundedTrace}
    H_t=
    \begin{cases}
      \alpha w_tw_t^T, & \text{if}\ w_t^TJw_t\geq 0 \\
      0, & \textup{otherwise}.
    \end{cases}
\end{equation}

The computational complexity of each iteration depends on computing a rank-1 matrix that satisfies this inequality, i.e., solving an approximate eigenvalue problem.
We use power method with random start to compute $\lambda_t$ and the unit vector $w_t$.~\citet*{kuczynski1992estimating} provide error bounds for power method when the input matrix is PSD. In Lemma~\ref{lemma:PowerMethod}, we restate their result so it applies to the largest eigenvalue of a symmetric matrix $J$.

\begin{lemma}\label{lemma:PowerMethod}
Let $J\in \mathbb{S}^n$ and $\lambda = \max_i |\lambda_i(J)|$. For $\delta \geq 0$, $\SubProb \in [0,1)$ and $\alpha \geq 0$, the power method with random start computes an unit vector $w$ that satisfies
\begin{equation}
\alpha w^TJw \geq \alpha \lambda_{\textup{max}}(J) - \delta
\end{equation}
with probability at least $1-\SubProb$ after $k \geq \frac{\lambda\alpha}{\delta}\log\left(\frac{n}{p^2}\right)$ iterations. Each iteration of the power method consists of a matrix-vector multiplication with $J$ and the working memory is $n$ numbers.
\end{lemma}

\paragraph{} Note that, for~\eqref{prob:PrimalSDP}, it it sufficient to store the input parameters, the map $\mathcal{B}(\cdot)$ and the rank-1 updates. If we have an access to a black box performing matrix-vector multiplications with $J$, the Frank-Wolfe algorithm, when applied to~\eqref{prob:PrimalSDP}, has working memory bounded by $\mathcal{O}(n+d)$ at each step.


\section{Frank-Wolfe Algorithm with Gaussian Sampling}\label{Sec:Methodology}
We now explain how to modify Algorithm~\ref{algo:FWapprox}, when applied to~\eqref{prob:PrimalSDP}, to replace the matrix-valued iterates with Gaussian vectors such that their covariance is equal to the iterate value. We then show how to apply this approach to more general SDPs by incorporating the constraints into the objective with a penalty.

\subsection{Idea of Gaussian Sampling}\label{sec:GS}
Consider the Frank-Wolfe update for~\eqref{prob:PrimalSDP} at iterate $t$, i.e., $X_{t+1} = (1-\gamma_t)X_t + \gamma_t H_t$. Assume that, at iterate $t$, we have a zero-mean random vector $z_t$ with covariance $X_t$. The update direction $H_t = w_tw_t^T$ has rank one, so if $\zeta \sim \mathcal{N}(0,1)$, then $\zeta w_t \sim \mathcal{N}(0,H_t)$. Now, if we define $z_{t+1} = \sqrt{1-\gamma_t}z_t+\sqrt{\gamma_t}\zeta w_t$, then $\mathbb{E}\left[ z_{t+1}z_{t+1}^T\right] = (1-\gamma_t)X_t + \gamma_t H_t = X_{t+1}$.

\begin{comment}
\begin{remark}
Let $H_t = w_tw_t^T$ be a rank-$1$ matrix. If $\zeta\sim \mathcal{N}(0,1)$ is a random variable, then $\zeta w_t\sim \mathcal{N}(0,H_t)$. Moreover, let $z_t\sim \mathcal{N}(0,X_t)$ be a Gaussian random vector in $\mathbb{R}^n$ that is independent of $\zeta w_t$ and $\gamma \in [0,1]$. If $z_{t+1} = \sqrt{1 - \gamma}z_t + \sqrt{\gamma}\zeta w_t$, then $z_{t+1} \sim \mathcal{N}(0, (1-\gamma)X_t+\gamma H_t) \sim \mathcal{N}(0,X_{t+1})$.
\label{Proposition:H_decomp}
\end{remark}

\begin{remark}
Let $H = ww^T$ be a rank-$1$ matrix. If $\zeta\sim \mathcal{N}(0,1)$ is a random variable, then $\zeta w\sim \mathcal{N}(0,H)$.
\end{remark}\label{Proposition:H_decomp}

At iterate $t$, if we have two independent samples $z_t\sim \mathcal{N}(0,X_t)$ and $\zeta w_t \sim \mathcal{N}(0,H_t)$, it is possible to generate a sample $z_{t+1} \sim \mathcal{N}(0,X_{t+1})$ without explicitly creating the matrix $X_{t+1}$ using the update technique in the following remark.

\begin{remark} 
\label{Proposition:Gaussian}
Let $z_t\sim \mathcal{N}(0,X_t)$, $\zeta w_t\sim \mathcal{N}(0,H_t)$ be independent Gaussian random vectors in $\mathbb{R}^n$ and $\gamma \in [0,1]$. If $z_{t+1} = \sqrt{1 - \gamma}z_t + \sqrt{\gamma}\zeta w_t$, then $z_{t+1} \sim \mathcal{N}(0, (1-\gamma)X_t+\gamma H_t) \sim \mathcal{N}(0,X_{t+1})$.
\end{remark}
\end{comment}

A nonnegative weighted sum of the samples then gives rise to nonnegative weighted sum of their covariance matrices. Furthermore, we can generate a sample at the next iterate in $\mathcal{O}(n)$ memory when the update has rank at most one.

\subsubsection{Frank-Wolfe algorithm with Gaussian sampling}

Algorithm~\ref{Algo1} incorporates Gaussian sampling into Algorithm~\ref{algo:FWapprox} by replacing the matrix variables by the sampled representation. The main algorithm is similar to that in Algorithm~\ref{algo:FWapprox}. The difference lies in the functions \texttt{LMO} and \texttt{UpdateVariable}:
\begin{itemize}
\item \texttt{LMO}: In Algorithm~\ref{algo:FWapprox}, this function simply computes the update direction by computing the approximate maximizer of $\left\langle s, \mathcal{B}^*(\nabla g(v))\right\rangle$ over the feasible set with additive error at most $\frac{1}{2}\eta\gamma C_g^u$. In Algorithm~\ref{Algo1}, we replace this step by computing the vector $w$ such that $w_tw_t^T$ solves the subproblem~\eqref{prob:Subproblem} approximately and the update satisfies~\eqref{eqn:HtforBoundedTrace}.
\item \texttt{UpdateVariable}: In Algorithm~\ref{algo:FWapprox}, we update the decision variable $x$ and its projection $v = \mathcal{B}(x)$ at every iteration. In Algorithm~\ref{Algo1}, we only track the Gaussian sample $z\sim \mathcal{N}(0,X)$ and $v = \mathcal{B}(X)$ which requires only $\mathcal{O}(n+d)$ memory instead of $\mathcal{O}(n^2+d)$.
\end{itemize}

\begin{comment}
Combining this with Proposition~\ref{Proposition:H_decomp}, we define the sampling technique as follows:

\begin{enumerate}
\item At step $0$, generate a sample $z_{0} \sim \mathcal{N}(0,X_0)$, $\ z_{0} \in \mathbb{R}^n$. Set $t = 0$.
\item Compute the eigenvector $u_t$ corresponding to maximum eigenvalue of $\mathcal{B}^*(\nabla g(v_t))$ for the update direction and generate a sample $w_t = u_tx$, where $x\sim \mathcal{N}(0,1)$.
\item Recall the update, $X_{t+1} = (1-\gamma)X_t + \gamma H_t$. Using Proposition~\ref{Proposition:Gaussian}, update $z_{t+1} = (\sqrt{1-\gamma})z_t + \sqrt{\gamma}w_t$ in the sketch framework. Set $t = t+1$. Go to step 2.
\end{enumerate}
\end{comment}


Algorithm~\ref{Algo1} gives a detailed description of the sampled modification of Frank-Wolfe applied to~\eqref{prob:PrimalSDP}. The convergence rate of the Frank-Wolfe algorithm given in Theorem~\ref{thm:FWapproxConv} applies even with the incorporation of Gaussian sampling, i.e., Algorithm~\ref{Algo1} converges to $z_t$ such that $\mathbb{E}\left[ z_tz_t^T\right]$ is $\epsilon$-optimal for~\eqref{prob:PrimalSDP} after $t\sim \mathcal{O}\left(\frac{C_g^u}{\epsilon}\right)$ iterations. We summarize the result of Algorithm~\ref{Algo1} applied to~\eqref{prob:PrimalSDP} in Proposition~\ref{prop:MainAlgoResult}. 

\begin{algorithm}[ht] 
\caption{Frank-Wolfe Algorithm with Gaussian Sampling} \label{Algo1}

\SetAlgoLined
\SetKwInOut{Input}{Input}\SetKwInOut{Output}{Output}
\DontPrintSemicolon

\Input{Input data for~\eqref{prob:PrimalSDP}, stopping criteria $\epsilon$, accuracy parameter $\eta$, probability $p$ for subproblem~\eqref{prob:Subproblem}, upper bound $C_g^u$ on the curvature constant}

\Output{Sample $z \sim \mathcal{N}(0,\widehat{X}_{\epsilon})$ and $v = \mathcal{B}(\widehat{X}_{\epsilon})$ such that $\widehat{X}_{\epsilon}$ is an $\epsilon$-optimal solution of~\eqref{prob:PrimalSDP}}
\BlankLine

 \SetKwFunction{CDD}{LMO}
 \SetKwFunction{UV}{UpdateVariable}
 \SetKwProg{Fn2}{Function}{:}{\KwRet{$z_t$}}
 \SetKwFunction{FWGS}{FWGaussian}
 \SetKwProg{Fn}{Function}{:}{\KwRet{$(z_t,v_t)$}}
 \Fn{\FWGS}{
 Select initial point $X_0 \in \mathcal{S}$; $X_0 \leftarrow \frac{\alpha}{n}I \textup{(say)}$ and set $v_0 \leftarrow \mathcal{B}(X_0)$\;
 Sample $z_{0} \sim \mathcal{N}(0,X_0)$\;
 $t \leftarrow 0$\;
 $\gamma \leftarrow 2/(t+2)$\;
 $(w_t,q_t) \leftarrow \CDD (\mathcal{B}^*(\nabla g(v_t)), \frac{1}{2}\eta\gamma C_g^u,p)$\;
 
 
 \While{$\left\langle q_t - v_t, \nabla g(v_t)\right\rangle > \epsilon$}{
 $(z_{t+1}, v_{t+1}) \leftarrow \UV (z_{t},v_{t},w_t,q_t,\gamma)$\;
 
 $t \leftarrow t+1$\;
 $\gamma \leftarrow 2/(t+2)$\;
 $(w_t,q_t) \leftarrow \CDD (\mathcal{B}^*(\nabla g(v_t)), \frac{1}{2}\eta\gamma C_g^u,p)$\;
 
 
 }
 }
 
 \SetKwProg{Fn}{Function}{:}{\KwRet{($w,q$)}}
 \Fn{\CDD{$J$, $\delta$, $p$}}{
 Find a unit vector $w$ such that with probability at least $1-p$, $\alpha\langle ww^T, J\rangle \geq \max_{d\in \mathcal{S}} \alpha\langle d, J\rangle - \delta$\;
 $\lambda \leftarrow \langle ww^T, J\rangle$\;
 \uIf{$\lambda \geq 0$}{
 $q \leftarrow \mathcal{B}(\alpha ww^T)$\;
 }
 \Else{
 $q \leftarrow 0$\;
 $w \leftarrow 0$\;
 }
 }
 
 \SetKwProg{Fn}{Function}{:}{\KwRet{($z,v$)}}
 \Fn{\UV{$z,v,w,q,\gamma$}}{
 $z \leftarrow (\sqrt{1-\gamma})z + \sqrt{\gamma}w\zeta\;\;\textup{where $\zeta\sim \mathcal{N}(0,1)$}$\;
 $v \leftarrow (1-\gamma)v + \gamma q $\;
 
 }

\end{algorithm}


\begin{prop}\label{prop:MainAlgoResult}
The output of Algorithm~\ref{Algo1} is a zero-mean Gaussian random vector  $\widehat{z}_{\epsilon}$ with covariance $\widehat{X}_{\epsilon}$, where $\widehat{X}_{\epsilon}$ is an  $\epsilon$-optimal solution of~\eqref{prob:PrimalSDP}. The working memory of the algorithm is $\mathcal{O}(n+d)$, where the $d$ is the dimension of the image of $\mathcal{B}$.
\end{prop}

\subsection{SDP with Linear Equality Constraints}\label{sec:LintoBounded}
Consider an SDP with linear objective function and linear equality constraints written as,
\begin{equation}\label{prob:SDPlin}
\tag{SDP}
\max_X\ \ \langle C, X \rangle\ \ \ \textup{subject to}\ \ \ \begin{cases}&\mathcal{A}(X) = b\\
&X \succeq 0,\end{cases}
\end{equation}
where $\mathcal{A}(\cdot): \mathbb{S}^{n} \rightarrow \mathbb{R}^d$ is a linear map. We assume that the feasible region is bounded. In Section~\ref{sec:SDPBoundedTrace}, we saw that when Algorithm~\ref{algo:FWapprox} is applied to~\eqref{prob:PrimalSDP}, the rank of the update variable is at most one. And in Section~\ref{sec:GS}, we saw that in this case we can update and track the change in a Gaussian sample representing the matrix-valued decision variable without explicitly computing the matrix at intermediate steps. Because of the simplicity of solving subproblem~\eqref{prob:Subproblem} and updating the samples for trace constrained problems, we penalize the linear equality constraints in~\eqref{prob:SDPlin} with a smooth penalty function so that the feasible domain is reduced to the one in~\eqref{prob:PrimalSDP}. By penalizing the constraints, we can only generate a near-feasible point to~\eqref{prob:PrimalSDP}, however, the extreme points of the modified constraint set now have a concise representation. In the specific case of~\eqref{prob:relaxedMC}, it is possible to generate a feasible solution with relative error bounds on the objective function value if we know a near-feasible solution for the problem with infeasibility error bounded by $\|\mathcal{A}(X)-b\|_{\infty}$. This motivates us to use a penalty function~\eqref{eqn:penaltyPhiDef} that approximates $\|\mathcal{A}(X)-b\|_{\infty}$.


\paragraph{Penalty function.} For $M>0$, let $\phi_M(\cdot):\mathbb{R}^{d}\rightarrow \mathbb{R}$ be defined by
\begin{equation}\label{eqn:penaltyPhiDef}
\phi_M(v) = \frac{1}{M}\log\left(\sum_{i=1}^{d} e^{M(v_i)} + \sum_{i=1}^{d} e^{M(-v_i)}\right).
\end{equation}

This function, also known as \textsc{LogSumExp} (LSE), is a smoothed approximation of $\|v\|_{\infty}$ as the next well-known proposition shows.

\begin{prop}[Bound on penalty~\cite{grigoriadis1994fast}]\label{prop:BoundPenalty}
If $\phi_M(\cdot)$ is defined as in~\eqref{eqn:penaltyPhiDef}, then 
\begin{equation}\label{eqn:BoundOnPenalty}
\| v\|_{\infty} \leq \phi_M (v) \leq \frac{\log(2d)}{M} + \| v\|_{\infty}.
\end{equation}
\end{prop}


We add this penalty term to the objective function of~\eqref{prob:SDPlin} to penalize the equality constraints and then solve the problem
\begin{equation} \label{prob:LogPenalty}
\tag{SDP-LSE}
\max_X\ \ \langle C,X\rangle - \beta \phi_M(\mathcal{A}(X)-b)\ \ \ \textup{subject to}\ \ \ \begin{cases}&\textup{Tr}(X) \leq \alpha,\\
& X \succeq 0,\end{cases}
\end{equation}
where $M$ and $\beta$ are positive constants to be chosen later, and $\alpha$ is chosen such that $X\succeq 0$ and $\mathcal{A}(X)=b$ implies that $\textup{Tr}(X)\leq \alpha$. This is possible because the feasible region is assumed to be bounded.


A similar approximation and penalty technique is used by~\citet*{hazan2008sparse} to compute an approximate solution of a feasibility problem with linear inequality constraints. The objective function in this case is simply $\phi_M(\mathcal{A}(X) - b)$ and the optimal objective value is zero.

Let $(u,v) = \mathcal{B}(X) = (\langle C,X\rangle, \mathcal{A}(X))$ so that the objective function of~\eqref{prob:LogPenalty} can be expressed as $g(u,v) = u - \beta \phi_M(v-b)$. This has the same structure as~\eqref{prob:PrimalSDP}. We can now use Algorithm~\ref{Algo1} to compute an $\epsilon$-optimal solution $\widehat{X}_{\epsilon}$ of~\eqref{prob:LogPenalty}. Moreover, we will show (in Lemma~\ref{lemma:SubOptFeas}) that by choosing the parameters $M$ and $\beta$ approximately, $\widehat{X}_{\epsilon}$ is also a near-feasible, near-optimal solution of~\eqref{prob:SDPlin}. The convergence rate and the bounds on the infeasibility and objective function value at $\widehat{X}_{\epsilon}$ for~\eqref{prob:SDPlin} are given in Lemma~\ref{lemma:SubOptFeas}.

\subsubsection{Convergence of the Frank-Wolfe algorithm} \label{sec:ConvergenceFW}
We have seen that Theorem~\ref{thm:FWapproxConv} states the convergence result for Algorithm~\ref{algo:FWapprox} when applied to~\eqref{prob:BoundedMapped}. The algorithm converges to an $\epsilon$-optimal solution after $\mathcal{O}\left( \frac{C_g^u}{\epsilon}\right)$ iterations. Moreover, this convergence result also holds for Algorithm~\ref{Algo1} when applied to~\eqref{prob:LogPenalty}. We now determine an upper bound on the curvature constant $C_g$ for~\eqref{prob:LogPenalty}.

\begin{lemma}\label{lemma:boundOnCf}
An upper bound on the curvature constant $C_g$ of the concave, smooth function $g(\mathcal{B}(X)) = \langle C, X\rangle - \beta \phi_M (\mathcal{A}(X)-b)$ over the compact, convex set $\mathcal{S} = \{X\succeq 0: \textup{Tr}(X) \leq \alpha\}$ is
\begin{equation*}
C_g \leq \beta \omega M \alpha^2,
\end{equation*}
where $\omega = \max_i \lambda_{\textup{max}}(A_i)$.
\end{lemma}

\begin{proof}
Let $f(X) = -\phi_M(\mathcal{A}(X)-b)$ and let $\mathcal{S}_f = \{X\succeq 0: \textup{Tr}(X) \leq 1\}$. An upper bound on the curvature constant $C_f$ of $f$ over $\mathcal{S}_f$, given by~\citet*{hazan2008sparse} is
\begin{equation}
C_f \leq \lambda_{\textup{max}}(-\nabla^2f)\textup{diam}(\mathcal{S}_f)^2 \leq \omega M.
\end{equation}

Note that $\lambda_{\textup{max}}(-\nabla^2g) = \beta \lambda_{\textup{max}}(-\nabla^2f)$ and $\textup{diam}(\mathcal{S})^2 =  \alpha^2\textup{diam}(\mathcal{S}_f)^2$. Hence, an upper bound on the curvature constant of $g$ over $\mathcal{S}$ is \[ C_g \leq \lambda_{\textup{max}}(-\nabla^2g)\textup{diam}(\mathcal{S})^2 = \beta \alpha^2 \lambda_{\textup{max}}(-\nabla^2f)\textup{diam}(\mathcal{S}_f)^2 \leq \beta \omega M\alpha^2. \qedhere \]
\end{proof}


\subsubsection{{Optimality and feasibility results for~\eqref{prob:SDPlin}}}
Given $\widehat{X}_{\epsilon}$, an $\epsilon$-optimal solution to~\eqref{prob:LogPenalty}, we analyze its suboptimality and infeasibility
with respect to~\eqref{prob:SDPlin}. The dual of~\eqref{prob:SDPlin} is
\begin{equation}\label{prob:SDPLinDual}
\tag{DSDP}
\min_y\ \ b^Ty\ \ \ \textup{subject to}\ \ \ \mathcal{A}^*(y) - C \succeq 0.
\end{equation}

We assume that~\eqref{prob:SDPlin} is feasible and a constraint qualification holds ensuring that strong duality is satisfied, and the primal and dual problems have finite optimal values. Let ($X^{\star}_{SDP},y^{\star}_{SDP}$) be a primal-dual optimal pair. The optimality and infeasibility bounds that we obtain depend on the properties of the optimal dual solution, $y^{\star}_{SDP}$.

\begin{lemma}\label{lemma:SubOptFeas}
Let $(X^{\star}_{SDP}, y^{\star}_{SDP})$ be a primal-dual optimal solution of~\eqref{prob:SDPlin}\textendash\eqref{prob:SDPLinDual} and let $\widehat{X}_{\epsilon}$ be an $\epsilon$-optimal solution of~\eqref{prob:LogPenalty}. If $\beta > \| y^{\star}_{SDP}\|_1$, then
\begin{equation}\label{eqn:OptBound}
\langle C, X^{\star}_{SDP}\rangle - \epsilon \leq \langle C, \widehat{X}_{\epsilon}\rangle \leq \langle C, X^{\star}_{SDP}\rangle + \| y^{\star}_{SDP}\|_1 \frac{\beta\frac{\log(2d)}{M} + \epsilon}{\beta - \| y^{\star}_{SDP}\|_1},
\end{equation}
and
\begin{equation}\label{eqn:FeasBound}
\| \mathcal{A}(\widehat{X}_{\epsilon}) -b \|_{\infty} \leq \frac{\beta\frac{\log(2d)}{M} + \epsilon}{\beta - \| y^{\star}_{SDP}\|_1}.
\end{equation}
\end{lemma}

\begin{proof}
The bounds~\eqref{eqn:OptBound} and~\eqref{eqn:FeasBound} are derived in Appendix~\ref{appendix:OptFeasLin}.
\end{proof}

\begin{remark}
Lemma~\ref{lemma:SubOptFeas} shows that $\widehat{X}_{\epsilon}$ is a near-feasible point to~\eqref{prob:SDPlin} with bounded objective function value such that the infeasibility and optimality bounds depend on the dual solution, $\| y^{\star}_{SDP}\|_1$, and the parameters $\beta$ and $M$. If the parameter values $\beta$ and $M$ are specifically chosen to be $\beta = 2\| y^{\star}_{SDP}\|_1$ and $M = 2\frac{\log(2d)}{\epsilon}$, then $\widehat{X}_{\epsilon}$ satisfies $\| \mathcal{A}(\widehat{X}_{\epsilon}) -b \|_{\infty} \leq 2\epsilon$ and the objective function value is upper bounded by $\langle C, X^{\star}_{SDP}\rangle + 2\epsilon\| y^{\star}_{SDP}\|_1$.
\end{remark}

\begin{remark}
In some cases, it is difficult to produce truly feasible points from the generated near-feasible points. While in other cases, such as~\eqref{prob:relaxedMC}, it is fairly straightforward to make small modifications to a generated near-feasible point to produce a feasible point with a similar objective value. Furthermore, we will see that for~\eqref{prob:relaxedMC}, it is possible to reduce $\| y^{\star}_{SDP}\|_1$ to $\langle C, X^{\star}_{SDP}\rangle$, eliminating the unknown optimal dual variable from our bounds.
\end{remark}

\begin{comment}
\paragraph{Computing a feasible solution of~\eqref{prob:SDPlin}.}
The feasible region of~\eqref{prob:SDPlin} is the intersection of the set $\{X:\mathcal{A}(X) = b\}$ with PSD cone.
We can find a feasible point for~\eqref{prob:SDPlin}, by computing the projection of a point that is near-feasible, namely $\widehat{X}_{\epsilon}$, onto this set. To do so, it is possible to generate a sequence of points, starting from $\widehat{X}_{\epsilon}$, by alternatively projecting on the set $\{X:\mathcal{A}(X) = b\}$ and the PSD cone. If there exists a strictly feasible (Slater) point for the set $\mathcal{S} = \{X:\mathcal{A}(X) = b, X\succeq 0\}$, i.e., there exists some $X_S\in \mathcal{S}$ such that $X \succ 0$, then this sequence converges finitely to a point in $\mathcal{S}$~\cite{bauschke1996projection}. However, this method is computationally challenging since involves solving a linear system of equations and computing eigendecomposition. In general, it is not possible to compute the projection on the PSD cone in low memory.
\end{comment}


\begin{comment}
\paragraph{Feasible solution of~\eqref{prob:SDPlin}.}
The solution $\widehat{X}_{\epsilon}$ is near-feasible to~\eqref{prob:SDPlin}. In general, it might not be possible to compute a projection on the feasible set and determine the change in objective function value easily. However, for a specific case, when the cost matrix $C$ is positive semidefinite and the constraint matrices $A_i = v_iv_i^T$ have rank-1, we can find a feasible solution simply by scaling the solution. Let the constraints be given as $v_i^TX v_i \leq 1\ \ \forall i\in \{1,\dotsc,d\}$. If $\| \mathcal{A}(\widehat{X}_{\epsilon}) -b \|_{\infty} = \epsilon$, we have $v_i^T\widehat{X}_{\epsilon} v_i \leq 1+ \epsilon\ \ \forall i$. So, by scaling each entry of $\widehat{X}_{\epsilon}$ by $\frac{1}{1+\epsilon}$, we find a feasible solution of the problem. Moreover, the objective function value at this feasible point is also scaled by $\frac{1}{1+\epsilon}$ (since $C$ is positive semidefinite) and the samples scaled by $\frac{1}{\sqrt{1+\epsilon}}$ correspond to the feasible solution.
\end{comment}

\begin{comment}
For certain problems, a feasible solution is required. However, simple rounding of the covariance matrices of the Gaussian samples could potentially violate positive semidefinite constraint. We can find the solution by projecting it onto the feasible region. This is explored further for \textsc{MaxCut} problem in the next section.
\end{comment}

\section{Approximation Algorithm for \textsc{MaxCut}}\label{Sec:MaxCut}
We now apply the general framework from Section~\ref{Sec:Methodology} to give an implementation of the Goemans-Williamson approximation algorithm for \textsc{MaxCut} that uses only $\mathcal{O}(n)$ working memory. Recall from Section~\ref{Sec:MaxCutIntro} that the standard SDP relaxation of \textsc{MaxCut},~\eqref{prob:relaxedMC} is a special case of~\eqref{prob:SDPlin} with a symmetric diagonally dominant cost matrix $C$ and the constraint set $\{X\in \mathbb{S}^n: X\succeq 0,\textup{diag}(X) = \mathbbm{1}\}$. A key challenge in applying the results of Lemma~\ref{lemma:SubOptFeas} in this setting is that the suboptimality and infeasibility bounds obtained depend on the optimal solution of the dual SDP~\eqref{prob:SDPLinDual}. With no \emph{a priori} knowledge of the dual solution, selecting the value of parameter $\beta$ is difficult and the additive error term in~\eqref{eqn:OptBound} could be quite high. Furthermore, setting the value of $\beta$ to an arbitrarily large value increases the curvature constant.

For~\eqref{prob:relaxedMC}, we show how to apply Lemma~\ref{lemma:SubOptFeas} to obtain relative error bounds on the suboptimality and infeasibility of the output of Algorithm~\ref{Algo1} without \textit{a priori} knowledge of the dual optimal solution. This allows us to appropriately choose the parameters $\beta$ and $M$ in the penalized formulation~\eqref{prob:MCpenalized}. Moreover, it is possible to generate a feasible solution to \textsc{MaxCut} that nearly achieves the approximation guarantee of Goemans-Williamson method, by applying Algorithm~\ref{algo:GenSample} (see Section~\ref{sec:MCRounding}) to the Gaussian samples generated from the PSD matrix produced by Algorithm~\ref{Algo1}.

Let the constraint `$\textup{diag}(X) = \mathbbm{1}$' in~\eqref{prob:relaxedMC} be penalized with $\phi_M(\cdot)$~\eqref{eqn:penaltyPhiDef} and consider the modified problem
\begin{equation}\label{prob:MCpenalized}
\tag{MaxCut-LSE}
\max_X\ \ \langle C, X\rangle - \beta \phi_M (\textup{diag}(X) -\mathbbm{1})\ \ \ \textup{subject to}\ \ \ \begin{cases} &\textup{Tr}(X) \leq n\\
&X\succeq 0.
\end{cases}
\end{equation}

We will see that choosing $M = 4\frac{\log(2d)}{\epsilon}$ and $\beta = 4\textup{Tr}(C)$ in~\eqref{prob:MCpenalized} allows us to derive a relative error bound on the objective function value of~\eqref{prob:relaxedMC}.  The main result of this section is given as follows:

\begin{theorem}\label{thm:MaxCutResult}
Let~\eqref{prob:MCpenalized} be solved to $\epsilon\textup{Tr}(C)$-optimality using Algorithm~\ref{Algo1} with $\epsilon \in \left(0,\frac{1}{3}\right)$, $\eta =1$, $p = \frac{\epsilon}{T(n,\epsilon)}$, where $T(n,\epsilon) = 64\frac{\log(2n)n^2}{\epsilon^2}$, followed by the rounding scheme of Algorithm~\ref{algo:GenSample}. For a diagonally dominant matrix $C$, this procedure generates a binary vector $w$ that satisfies
\begin{equation}\label{eqn:MCoptimality}
\alpha_{GW}(1-3\epsilon) \textup{opt} \leq \mathbb{E}[w^TCw] \leq \textup{opt},
\end{equation}
where $\textup{opt}$ is the maximum of $w^TCw$ over the set $w \in \{-1,1\}^n$. Algorithm~\ref{Algo1} terminates after at most $T(n,\epsilon)$ iterations, where each iteration performs at most $280 \frac{n}{\epsilon}\log\left(\frac{2n}{\epsilon}\right)$ matrix-vector multiplications. The working memory required is at most $3n$ numbers.
\end{theorem}

\subsection{Relative Error Bounds on Suboptimality and Infeasibility for~\eqref{prob:relaxedMC}}
In this subsection, we derive relative error bounds on the suboptimality and infeasibility for~\eqref{prob:relaxedMC} at $\widehat{X}_{\epsilon}$, the output of Algorithm~\ref{Algo1}, when solving~\eqref{prob:MCpenalized}. This is an application of Lemma~\ref{lemma:SubOptFeas} with the key difference that any dependence on the dual optimal solution has been eliminated.

The dual of SDP relaxation of \textsc{MaxCut} is,
\begin{equation}\label{prob:MCSDPDual}
\min_y\ \ \sum_{i=1}^n y_i\ \ \ \textup{subject to}\ \ \ \textup{diag}^*(y) - C \succeq 0.
\end{equation}

Let $y^\star_{SDP}$ be an optimal solution of Problem~\eqref{prob:MCSDPDual}.
We first show that $\| y^{\star}_{SDP} \|_1$ is upper bounded by $2\textup{Tr}(C)$.

\begin{lemma}\label{lemma:MaxCutSDPbounds}
Let $(X^\star_{SDP}, y^{\star}_{SDP})$ be a primal-dual optimal pair for~\eqref{prob:relaxedMC} and its dual~\eqref{prob:MCSDPDual}. If $C$ is diagonally dominant, then
\begin{equation}
\textup{Tr}(C) \leq \langle C, X^{\star}_{SDP}\rangle = \| y^{\star}_{SDP}\|_1 \leq 2 \textup{Tr}(C).
\end{equation}
\end{lemma}

\begin{proof}
For a symmetric diagonally dominant cost matrix $C$ with nonnegative entries on the diagonal, it follows from the Gershgorin cirle theorem, that $C$ must be a PSD matrix. Since $C$ is PSD, $\textup{diag}^*(y) - C \succeq 0$ implies $y \geq 0$. Thus, the objective function of Problem~\eqref{prob:MCSDPDual} can be written as $\| y\|_1$. Moreover, the SDP relaxation of \textsc{MaxCut} satisfies Slater's condition, so
\begin{equation*}
\langle C, X^{\star}_{SDP}\rangle = \| y^{\star}_{SDP}\|_1,
\end{equation*}
for a primal-dual optimal pair $(X^{\star}_{SDP},y^{\star}_{SDP})$. To see that $\textup{Tr}(C) = \langle C, I \rangle \leq \langle C, X^{\star}_{SDP}\rangle$, we simply note that identity matrix $I$ is feasible for~\eqref{prob:relaxedMC}.

To prove $\langle C, X^{\star}_{SDP}\rangle \leq 2\textup{Tr}(C)$, we use the fact that $C$ is a diagonally dominant matrix with nonnegative entries on the diagonal. As such the matrix $2\textup{diag}^*(\textup{diag}(C)) - C$ is also symmetric diagonally dominant and has nonnegative diagonal entries. It follows that
\begin{equation}
2\textup{diag}^*(\textup{diag}(C)) - C \succeq 0.
\end{equation}
Then,
\begin{equation}
\langle C, X^{\star}_{SDP}\rangle \leq  \langle 2\textup{diag}^*(\textup{diag}(C)), X^{\star}_{SDP}\rangle = 2\textup{Tr}(C),
\end{equation}
where we used the fact that $X^{\star}_{SDP} \succeq 0$ and $\textup{diag}(X^{\star}_{SDP}) = \mathbbm{1}$.
\end{proof}

\begin{remark}
If $C \succeq 0$, but not diagonally dominant, we have $\| y^{\star}_{SDP}\|_1 \leq n \lambda_{\textup{max}}(C) \leq n\textup{Tr}(C)$, since $y = \mathbbm{1}\cdot\lambda_{\textup{max}}(C)$ is feasible for the dual problem~\eqref{prob:MCSDPDual}. So, $\textup{Tr}(C) \leq \langle C, X^{\star}_{SDP}\rangle = \| y^{\star}_{SDP}\|_1 \leq n\textup{Tr}(C)$ in that case.
\end{remark}

\subsubsection{Optimality and feasibility bounds for~\eqref{prob:relaxedMC}} \label{sec:ConvergenceFW2}
By finding a near-optimal point for~\eqref{prob:MCpenalized}, the penalized relaxation of~\eqref{prob:relaxedMC}, we can obtain a near-feasible solution to~\eqref{prob:relaxedMC} that has relative error~$\epsilon$. Note that a relative error bound is exactly what we need to obtain a multiplicative approximation guarantee for \textsc{MaxCut}.

\begin{lemma}\label{lemma:SubOptFeasMaxCut}
Let $X^{\star}_{SDP}$ be an optimal solution of ~\eqref{prob:relaxedMC} and $\widehat{X}_{\epsilon}$ be an $\epsilon\textup{Tr}(C)$-optimal solution to~\eqref{prob:MCpenalized} with $M = 4 \frac{\log(2d)}{\epsilon}$ and $\beta = 4\textup{Tr}(C)$. Then
\begin{equation}\label{eqn:MCbounds}
\langle C, X^{\star}_{SDP}\rangle(1 - \epsilon) \leq \langle C, \widehat{X}_{\epsilon}\rangle \leq \langle C, X^{\star}_{SDP}\rangle (1 + \epsilon),
\end{equation}
and
\begin{equation}\label{eqn:MCinfeasbounds}
\| \mathcal{A}(\widehat{X}_{\epsilon}) -b \|_{\infty} \leq \epsilon.
\end{equation}
\end{lemma}

\begin{proof}
This result is an application of Lemma~\ref{lemma:SubOptFeas} for specific parameter values. Substituting the values of $M$ and $\beta$, and 
using the inequality $\|y^\star_{SDP}\|_1 \leq 2\textup{Tr}(C)$ from Lemma~\ref{lemma:MaxCutSDPbounds}, we see that 
\[ \| \mathcal{A}(\widehat{X}_{\epsilon}) - b\|_{\infty} \leq  \frac{\beta\frac{\log(2d)}{M} + \epsilon \textup{Tr}(C)}{\beta - \|y_{SDP}^\star\|_1} \leq \frac{2\epsilon \textup{Tr}(C)}{4\textup{Tr}(C) - 2\textup{Tr}(C)} = \epsilon.\]
Furthermore, combining the result $\langle C,X^\star_{SDP}\rangle = \|y^\star_{SDP}\|_1$ (from Lemma~\ref{lemma:MaxCutSDPbounds}) with Lemma~\ref{lemma:SubOptFeas} gives the upper bound on $\langle C, \widehat{X}_{\epsilon}\rangle$. Finally, using the fact that $\langle C, X^\star_{SDP}\rangle \geq \textup{Tr}(C)$ and substituting in~\eqref{eqn:OptBound}, gives
\begin{equation*}
\langle C, \widehat{X}_{\epsilon}\rangle \geq \langle C, X^\star_{SDP}\rangle - \epsilon \textup{Tr}(C) \geq \langle C, X^\star_{SDP}\rangle (1-\epsilon).\qedhere
\end{equation*}
\end{proof}

\begin{comment}
The proof is similar to the proof of Lemma~\ref{lemma:SubOptFeas}.

We get $\| \mathcal{A}(\widehat{X}_{\epsilon}) -b \|_{\infty} \leq \epsilon$ by substituting the values of $\beta$, $M$ and changing $\epsilon$ to $\epsilon \textup{Tr}(C)$ in~\eqref{eqn:BoundonInfeas}. Similarly, replacing $\epsilon$ by $\epsilon\textup{Tr}(C)$ in lower bound of~\eqref{eqn:OptBound} and using lower bound on $\langle C, X^{\star}_{SDP}\rangle$ from Lemma~\ref{lemma:MaxCutSDPbounds}, we the lower bound on $\langle C, \widehat{X}_{\epsilon}\rangle$.

We replace $\| y^{\star}_{SDP} \|_1$ in proof of Lemma~\ref{lemma:SubOptFeas} by its upper bound from Lemma~\ref{lemma:MaxCutSDPbounds} and the desired upper bound on $\langle C, \widehat{X}_{\epsilon}\rangle$.
\end{comment}


\begin{remark}
If $C \succeq 0$ (not necessarily diagonally dominant), then the bounds in Lemma~\ref{lemma:SubOptFeasMaxCut} hold for $M = (n+2)\frac{\log(2d)}{\epsilon}$ and $\beta = (n+2)\textup{Tr}(C)$.
\end{remark}


\subsection{Generating a Feasible Solution to~\textsc{MaxCut}}\label{sec:MCRounding}
We now show how to adapt the rounding procedure of Goemans-Williamson to our setting (Algorithm~\ref{algo:GenSample}). The reason we need to modify the Goemans-Williamson scheme is because the zero-mean Gaussian random vector returned by Algorithm~\ref{Algo1} has covariance that is not feasible for~\eqref{prob:relaxedMC}.

\begin{algorithm}[ht] 
\SetAlgoLined
\SetKwInOut{Input}{Input}\SetKwInOut{Output}{Output}
\DontPrintSemicolon

\Input{A sample $\widehat{z}_{\epsilon} \sim \mathcal{N}(0,\widehat{X}_{\epsilon})$ and $\textup{diag}(\widehat{X}_{\epsilon})$}
\Output{A feasible solution to \textsc{MaxCut} ${w} = \textup{sign}(\overline{w})$}
\BlankLine

\SetKwFunction{GS}{GenerateSample}
\SetKwProg{Fn}{Function}{:}{\KwRet{$\textup{sign}(\overline{w})$}}

 \Fn{\GS}{
 Generate $\zeta \sim \mathcal{N}\left(0, I - \textup{diag}^*\left(\frac{\textup{diag}(\widehat{X}_{\epsilon})}{\max (\textup{diag}(\widehat{X}_{\epsilon}))}\right)\right)$\;
 Set $\overline{w} = \frac{\widehat{z}_{\epsilon}}{\sqrt{\max (\textup{diag}(\widehat{X}_{\epsilon}))}} + \zeta$
 }
\caption{\label{algo:GenSample}Generate a binary vector from a Gaussian vector}
\end{algorithm}

Algorithm~\ref{algo:GenSample} can be used to generate a feasible solution of \textsc{MaxCut} from any PSD (covariance) matrix $X$. The first step of the algorithm generates $n$ independent zero-mean random variables with covariance defined by the diagonal entries of $I - \textup{diag}^*\left(\frac{\textup{diag}(\widehat{X}_{\epsilon})}{\max (\textup{diag}(\widehat{X}_{\epsilon}))}\right)$. The random vector $\overline{w}$ in step 3 is a sum of two independent zero-mean Gaussian random vectors. The covariance of this random vector $\overline{w}$ can be stated as
\begin{equation}\label{eqn:GenerateFeasMatrix}
\overline{X} = \frac{\widehat{X}_{\epsilon}}{\max (\textup{diag}(\widehat{X}_{\epsilon}))} + \left(I - \textup{diag}^*\left(\frac{\textup{diag}(\widehat{X}_{\epsilon})}{\max (\textup{diag}(\widehat{X}_{\epsilon}))}\right)\right)
\end{equation}
so that $\overline{w} \sim \mathcal{N}(0,\overline{X})$. The matrix $\overline{X}$ is a sum of two PSD matrices and so is PSD. Moreover, $\textup{diag}(\overline{X}) = \mathbbm{1}$, so $\overline{X}$ is feasible for~\eqref{prob:relaxedMC}. We can then apply the standard analysis of the Goemans-Williamson rounding scheme to $\overline{X}$.

\paragraph{Goemans-Williamson rounding.} For a PSD matrix $C$ and a Gaussian random vector $\overline{w} \sim \mathcal{N}(0,\overline{X})$, such that $\textup{diag}(\overline{X}) = \mathbbm{1}$,~\citet*{nesterov1998semidefinite} derived a $\frac{2}{\pi}$-approximation bound,
\begin{equation}\label{eqn:GWBound}
\mathbb{E}_G[w^TCw] \geq \alpha \langle C, \overline{X}\rangle,
\end{equation}
where $w = \textup{sign}(\overline{w})$, $\alpha = \frac{2}{\pi}$ and $\mathbb{E}_G[\cdot]$ denotes the expectation over Gaussian random vectors. Moreover, if $C$ is diagonally dominant,~\citet*{GW} provide a tighter bound with $\alpha = \alpha_{GW} \approx 0.878$.

When the input of Algorithm~\ref{algo:GenSample} is an approximate solution of~\eqref{prob:relaxedMC}, we analyze the expected objective value of $\textup{sign}(\overline{w})$.

\begin{lemma}\label{lemma:MaxCutapproxSol}
Let $\epsilon \in \left(0,\frac{1}{2}\right)$, $C$ be a diagonally dominant matrix, and let $\widehat{X}_{\epsilon} \succeq 0$ satisfy the bounds given in Lemma~\ref{lemma:SubOptFeasMaxCut}. If a binary vector $w = \textup{sign}(\overline{w})$ is generated by Algorithm~\ref{algo:GenSample} with input $\widehat{z}_{\epsilon} \sim \mathcal{N}(0, \widehat{X}_{\epsilon})$, then the expected value of $w^TCw$ satisfies
\begin{equation}
\alpha_{GW} (1-2\epsilon)\textup{opt} \leq \mathbb{E}_G[w^TCw] \leq \textup{opt} \leq \langle C, X^\star_{SDP}\rangle,
\end{equation}
where $\textup{opt}$ is the optimal value of $w^TCw$ over the set $w \in \{\pm 1 \}^n$.
\end{lemma}

\begin{comment}
\begin{proof}
The objective function value of SDP relaxation of \textsc{MaxCut} at $\overline{X}$ is
\begin{align}
\langle C, \overline{X}\rangle &= \left\langle C, \frac{\widehat{X}_{\epsilon}}{\max (\textup{diag}(\widehat{X}_{\epsilon}))} + \left(I - \textup{diag}^*\left(\frac{\textup{diag}(\widehat{X}_{\epsilon})}{\max (\textup{diag}(\widehat{X}_{\epsilon}))}\right)\right)\right\rangle\\
&\geq \frac{\langle C, \widehat{X}_{\epsilon}\rangle}{\max (\textup{diag}(\widehat{X}_{\epsilon}))}\label{eqn:MCApprox1}.
\end{align}

In~\eqref{eqn:MCApprox1}, we use the fact that both $C$ and $I -  \textup{diag}^*\left(\frac{\textup{diag}(\widehat{X}_{\epsilon})}{\max (\textup{diag}(\widehat{X}_{\epsilon}))}\right)$ are PSD and their inner product is greater than 0.  

Define an event $\mathcal{E} = \{\max (\textup{diag}(\widehat{X}_{\epsilon})) \leq 1+\epsilon\}$. Note that this event occurs with probability $1-Tp$. So,
\begin{align}
\mathbb{E}\left[\left\langle C,\overline{X}\right\rangle\right] &\geq \mathbb{E}\left[\frac{\langle C, \widehat{X}_{\epsilon}\rangle}{\max (\textup{diag}(\widehat{X}_{\epsilon}))}\right]\\
&= \mathbb{E}\left[\frac{\langle C, \widehat{X}_{\epsilon}\rangle}{\max (\textup{diag}(\widehat{X}_{\epsilon}))}| \mathcal{E}\right]\textup{Pr}(\mathcal{E}) + \mathbb{E}\left[\frac{\langle C, \widehat{X}_{\epsilon}\rangle}{\max (\textup{diag}(\widehat{X}_{\epsilon}))}| \mathcal{E}^C\right]\textup{Pr}(\mathcal{E}^C)\\
&\geq \mathbb{E}\left[\frac{\langle C, \widehat{X}_{\epsilon}\rangle}{\max (\textup{diag}(\widehat{X}_{\epsilon}))}| \mathcal{E}\right]\textup{Pr}(\mathcal{E})\label{eqn:expectedMC1}\\
&\geq \frac{(1-Tp)(1-\epsilon)}{1+\epsilon}\mathbb{E}\left[\langle C, {X}^{\star}_{\textup{SDP}}\rangle\right]\label{eqn:expectedMC2}.
\end{align}

Equation~\eqref{eqn:expectedMC1} holds since both the quantities in the second term are nonnegative and~\eqref{eqn:expectedMC2} follows from Lemma~\ref{lemma:SubOptFeasMaxCut}

Let $\epsilon = Tp$. Since $\frac{1}{1+\epsilon}\geq 1- \epsilon$,
\begin{align}
\mathbb{E}\left[\left\langle C,\overline{X}\right\rangle\right] &\geq (1-\epsilon)^3\mathbb{E}\left[\langle C, {X}^{\star}_{\textup{SDP}}\rangle\right]\\
&\geq (1-3\epsilon)\langle C, {X}^{\star}_{\textup{SDP}}\rangle\label{eqn:expectedMC3}.
\end{align}

Substituting the result of~\eqref{eqn:expectedMC3} in~\eqref{eqn:GWBound} gives the desired result.
The upper bound $w^{\star T}Cw^{\star}$ on $ \mathbb{E}[w^TCw]$ is satisfied since $w$ is feasible to \textsc{MaxCut} and the relaxed bound $\langle C, X^\star_{SDP}\rangle$ follows since~\eqref{prob:relaxedMC} is a relaxation of \textsc{MaxCut}.
\end{proof}
\end{comment}

\begin{proof}
The objective function value of~\eqref{prob:relaxedMC} at $\overline{X}$ is
\begin{align}
\langle C, \overline{X}\rangle &= \left\langle C, \frac{\widehat{X}_{\epsilon}}{\max (\textup{diag}(\widehat{X}_{\epsilon}))} + \left(I - \textup{diag}^*\left(\frac{\textup{diag}(\widehat{X}_{\epsilon})}{\max (\textup{diag}(\widehat{X}_{\epsilon}))}\right)\right)\right\rangle\\
&\geq \frac{\langle C, \widehat{X}_{\epsilon}\rangle}{\max (\textup{diag}(\widehat{X}_{\epsilon}))}\label{eqn:MCApprox1}\\
&\geq \frac{1-\epsilon}{1+\epsilon}\langle C, X^{\star}_{SDP}\rangle \label{eqn:MCApprox2}\\
&\geq (1-2\epsilon)\langle C, X^{\star}_{SDP}\rangle, \label{eqn:MCApprox3}
\end{align}

where~\eqref{eqn:MCApprox1} follows from the fact that both $C$ and $I -  \textup{diag}^*\left(\frac{\textup{diag}(\widehat{X}_{\epsilon})}{\max (\textup{diag}(\widehat{X}_{\epsilon}))}\right)$ are PSD and their inner product is greater than 0,~\eqref{eqn:MCApprox2} follows from Lemma~\ref{lemma:SubOptFeasMaxCut}, and~\eqref{eqn:MCApprox3} uses the fact that $\frac{1}{1+\epsilon} \geq 1-\epsilon$ and $(1-\epsilon)^2 \geq 1-2\epsilon$. Substituting~\eqref{eqn:MCApprox3} in~\eqref{eqn:GWBound} gives the desired result.
\end{proof}

Note that the result in Lemma~\ref{lemma:MaxCutapproxSol} holds irrespective of the algorithm used to compute $\widehat{X}_{\epsilon}$. We are now in a position to prove Theorem~\ref{thm:MaxCutResult}.

\begin{proof}[Proof of Theorem~\ref{thm:MaxCutResult}]
Since we use Algorithm~\ref{Algo1} to solve~\eqref{prob:MCpenalized} with $p = \frac{\epsilon}{T(n,\epsilon)}$, the bounds in Lemma~\ref{lemma:SubOptFeasMaxCut} are satisfied with probability at least $1-\epsilon$. Thus, the bound $\langle C, \overline{X}\rangle \geq (1-2\epsilon)\langle C, X^{\star}_{SDP}\rangle$ also holds with probability at least $1-\epsilon$. Moreover, for any $\overline{X} \succeq 0$, $\langle C, \overline{X}\rangle \geq 0$, and thus, a lower bound on the expected value of $\langle C, \overline{X}\rangle$ over the random initialization of the power method, i.e., $\mathbb{E}_P[\langle C, \overline{X}\rangle]$, at each iteration is,
\begin{align*}
\mathbb{E}_P[\langle C, \overline{X}\rangle] &\geq (1-2\epsilon)\langle C, X^{\star}_{SDP}\rangle
(1-\epsilon)\\
&\geq (1-3\epsilon) \langle C, X^{\star}_{SDP}\rangle.
\end{align*}

The lower bound in~\eqref{eqn:MCoptimality} then follows from~\eqref{eqn:GWBound} because
\begin{align*}
\mathbb{E}[w^TCw] &= \mathbb{E}_P[\mathbb{E}_G[w^TCw]]\\
&\geq \alpha_{GW}\mathbb{E}_P[\langle C, \overline{X}\rangle]\\
&\geq \alpha_{GW}(1-3\epsilon)\langle C, X^{\star}_{SDP}\rangle\\
&\geq \alpha_{GW}(1-3\epsilon)\textup{opt}.
\end{align*}

\paragraph{Bound on $T$, number of iterations of Algorithm~\ref{Algo1}.}
An upper bound on the curvature constant of~\eqref{prob:MCpenalized} is $C_g^u = 16\frac{\textup{Tr}(C)\log(2n)n^2}{\epsilon}$ since $\omega = 1$ and $\alpha = n$. Algorithm~\ref{Algo1} converges to an $\epsilon \textup{Tr}(C)$-optimal solution after at most $T =\frac{2C_g^u(1+\eta)}{\epsilon \textup{Tr}(C)} - 2 \leq 64 \frac{\log(2n)n^2}{\epsilon^2}$ iterations with probability at least $1-Tp$.

\paragraph{Bound on number of iterations of power method at each $t$.}
From Lemma~\ref{lemma:PowerMethod}, the number of matrix-vector multiplications performed at iteration $t$ of Algorithm~\ref{Algo1} is at most $\frac{\lambda\alpha}{\delta}\log\left( \frac{n}{p^2}\right)$ with $\delta = \frac{1}{2}\eta \gamma_t C_g^u \approx \frac{1}{8}\epsilon\textup{Tr}(C)$ and
\begin{align*}
\lambda &= \max_i |\lambda_i (\nabla g(v_t))|\\
&= \max_i |\lambda_i (C - \beta D)|,
\end{align*}
where $D$ is a diagonal matrix with $d_{ii} \in [-1,1]$.
Thus, $\lambda \leq \lambda_{\textup{max}}(C) + \beta \leq 5\textup{Tr}(C)$.

Substituting the value of $p$, and bounds on $\lambda$ and $\delta$ in Lemma~\ref{lemma:PowerMethod}, the number of iterations performed by the power method and thus, the number of matrix-vector computations at each $t$ is bounded by $280 \frac{n}{\epsilon}\log\left(\frac{2n}{\epsilon}\right)$. Furthermore, at iteration $t$ of Algorithm~\ref{Algo1}, we keep track of the sample $z$ and $\mathcal{A}(X)$ which requires storage of $2n$ numbers. Furthermore, the working memory of the power method is $n$ numbers. This leads to a total working memory of at most $3n$ numbers.
\end{proof}

\begin{remark}
If $C\succeq 0$, but not diagonally dominant, then the suboptimality bound~\eqref{eqn:MCoptimality} holds for $M = (n+2)\frac{\log(2d)}{\epsilon}$ and $\beta = (n+2)\textup{Tr}(C)$.
However, due to the dependence of the values of parameters $\beta$ and $M$ on $n$ in this case, it takes $\mathcal{O}\left(\frac{\log(2n)n^4}{\epsilon^2}\right)$ iterations within Algorithm~\ref{Algo1} to achieve the stated guarantee.
\end{remark}
{"config": "arxiv", "file": "2006.10945/Methodology.tex"}