Revisiting Spectral Bundle Methods: Primal-dual (Sub)linear Convergence Rates
The spectral bundle method proposed by Helmberg and Rendl is well established for solving large-scale semidefinite programs (SDP) thanks to its low per iteration computational complexity and strong practical performance. In this paper, we revisit this classic method showing it achieves sublinear convergence rates in terms of both primal and dual SDPs under merely strong duality, complementing previous guarantees on primal-dual convergence. Moreover, we show the method speeds up to linear convergence if (1) structurally, the SDP admits strict complementarity, and (2) algorithmically, the bundle method captures the rank of the optimal solutions. Such complementary and low rank structure is prevalent in many modern and classical applications. The linear convergent result is established via an eigenvalue approximation lemma which might be of independent interest. Numerically, we confirm our theoretical findings that the spectral bundle method, for modern and classical applications, speeds up under these conditions. Finally, we show that the spectral bundle method combined with a recent matrix sketching technique is able to solve an SDP with billions of decision variables in a matter of minutes.
\section{Analysis}\label{sec: analysis}
In this section, we present and derive our sublinear convergence guarantees for both Block-Spec and HR-Spec as well as our improved local linear convergence for Block-Spec under extra condition of strict complementarity (defined momentarily).
Key structural lemmas for our proofs are discussed in Section \ref{sec: analtical conditon}. 
Next, in Section \ref{sec: Relating primal and dual convergence}, we describe the relationship between primal and dual convergence for bundle methods. From this, we are able to focus our analysis on dual convergence only and yet conclude convergence for both. 


Let us first define the standard conditions of strong duality and strict complementarity. 

\paragraph{Strong duality} Strong duality in this paper means that the solution set $\xsolset$ and $\ysolset$ are nonempty and 
compact, and 
there is a solution pair $(\xsol,\ysol)\in \xsolset\times \ysolset$ such that the duality 
gaps is zero:
\[
\pval := \inprod{-C}{\xsol} = \inprod{-b}{\ysol}=:\dval.
\]

Note that we require $\xsolset$ and $\ysolset$ to be nonempty and compact instead of just  $\pval =\dval$. Such 
condition is ensured under Slater's conditions for both \eqref{p} and \eqref{d}. Next, we spell out strict complementarity. 

\paragraph{Strict Complementarity} \cite[Definition 4]{alizadeh1997complementarity} A pair of primal dual solution $(\xsol,\ysol)$ with dual slack matrix $\zsol = C-\Ajmap(\ysol)$ is \emph{strict complementary} if 
\[
\rank(\xsol)+\rank(\zsol)=\dm. 
\]
If \eqref{p} and \eqref{d} admit such pair, we say the pair \eqref{p} and \eqref{d} (or simply \eqref{p}) satisfies 
strict complementarity. Such condition is satisfied for a generic SDP \cite{alizadeh1997complementarity} and also many structural SDPs \cite{ding2020regularity}.


From these two conditions, we prove the following sublinear guarantees in Section \ref{sec: SublinearRates}\footnote{The proof of these rates rely on existing $\bigO \left(\frac{\log \frac{1}{\epsilon}}{\epsilon^{3}}\right)$ rates for the bundle method. Improvements in the underlying convergence theory for classic bundle methods would immediately give improvements in the convergence rates presented here.}. Denote 
$D_{\xsolset} :\,=\sup_{\xsol\in\xsolset} \nucnorm{\xsol}$, the largest nuclear norm of the primal solution set.
\begin{theorem}\label{thm: sublinearates}
	Suppose strong duality holds. Given any $\beta\in(0,1)$, $\bar{r}\geq 1$, $\rho_t = \rho>0$,
	$\alpha \geq 2D_{\xsolset}$, and 
	$z_0=y_0\in \RR^{\ncons}$, 
	there is a constant $\kappa$, such that for any $\epsilon\in (0,\frac{1}{2})$, both  Block-Spec and HR-Spec
	produce a solution pair $X_t$ and $y_t$ with
	\[
	t \leq \kappa \frac{ \log (\frac{1}{\epsilon})}{\epsilon^3}
	\]
	that satisfies $F(y_t)-F(\ysol)\leq \epsilon$, and 
	\begin{align*}
	\text{approximate primal feasibility: }& \quad\|b - \Amap X_{t}\|^2 \leq \epsilon, \quad  X_{t}\succeq 0,\\
	\text{approximate dual feasibility: }& \quad \lambda_{\min}(C - \Ajmap y_{t})\geq -\epsilon,\\
	\text{approximate primal-dual optimality: }&\quad  |\langle b, y_t \rangle - \langle C, X_t\rangle| \leq \sqrt{\epsilon}.
	\end{align*}
	Moreover, if strict complementarity holds, then 
	the iteration number $t$ is bounded by 
	\[
	t\leq \kappa \frac{\log(\frac{1}{\epsilon})}{\epsilon}.
	\]
\end{theorem}

We improve on these results by showing Block-Spec locally converges linearly whenever both assumptions of strong duality and strict complementarity hold. 
Additionally, we require $\bar{r}\geq r_d$, where $r_d =\max_{y\in \ysolset}\dim(\nullspace(Z(\ysol)))$, 
the largest dimension of the null space of dual slack matrices. 
The linear convergence is local because it  
happens after $y_t$ is close to $\ysolset$.  Upon first read, 
one can treat the dual solution as unique to ease the understanding. 
In this setting, $r_d = \dim(\nullspace(Z(\ysol))$.

We are now ready to state the our improved convergence theorem, which is proven in Section \ref{sec: Proof of Local linear convergence of Block bundle method}.


\begin{theorem}\label{thm: linear convergence of Block SBM under the extra condition strict complementarity}
	Suppose strong duality and strict complementarity holds. 
	Then under proper selection of $\rho$ and $T_0>0$, for any $\beta\in [0,\frac{1}{2}]$ and $\bar{r}\geq r_d$, after at most
	$T_0$ many steps, Block-Spec will subsequently only take descent steps and converge linearly to an optimal solution.
	Consequently, there is a constant $\varkappa$, such that for any $\epsilon\in(0,\frac{1}{2})$, Block-Spec will produce a solution pair $X_t$ and $y_t$ with
	\[
	t \leq T_0 + \varkappa \log (\frac{1}{\epsilon})
	\]
	that satisfies $F(y_t)-F(\ysol)\leq \epsilon$, and 
	\begin{align*}
	\text{approximate primal feasibility: }& \quad\|b - \Amap X_{t}\|^2 \leq \epsilon, \quad  X_{t}\succeq 0,\\
	\text{approximate dual feasibility: }& \quad \lambda_{\min}(C - \Ajmap y_{t})\geq -\epsilon,\\
	\text{approximate primal-dual optimality: }&\quad  |\langle b, y_t \rangle - \langle C, X_t\rangle| \leq \sqrt{\epsilon}.
	\end{align*}
\end{theorem}

The constants $\kappa$ and 
$\varkappa$ in Theorem
\ref{thm: sublinearates} and 
\ref{thm: linear convergence of Block SBM under the extra condition strict complementarity} 
can depend on \eqref{p}, the choice of $\alpha,\beta,$ and $\rho$, and the initialization $X_0,z_0,y_0$, but 
does not depend on $\epsilon$. The choice of $T_0$ and $\rho$ can be found in Section \ref{sec: Proof of Local linear convergence of Block bundle method}.
\subsection{Structural Lemmas}\label{sec: analtical conditon}

We introduce a growth lemma that is central to our proof of faster convergence rate. We call $\gamma_1$ the 
\emph{quadratic growth parameter} when $\zeta_1=2$ in the following lemma. 
\begin{lemma}[Quadratic Growth]\cite[Section 4]{sturm2000error} \label{lem: qg}
Suppose strong duality holds for \eqref{p} and \eqref{d}, then  there are some $\rho_1,\rho_2$, such that for any fixed $\epsilon>0$, there are some $\gamma_1,\gamma_2>0$ such that for all $y$ with $F(y)\leq \inprod{-b}{\ysol}+ \epsilon$, and all $X\succeq 0$ with 
$|\inprod{C}{X}- \inprod{C}{\xsol}|\leq \epsilon$ and $\twonorm{\Amap{X}-b}\leq \epsilon$:
	\[
\dist^{\zeta_1}(y,\ysolset) \leq  \gamma _1(F(y)-F(\ysol)),\quad \dist^{\zeta_2}(X,\xsolset)\leq \gamma_2\left(\abs{\inprod{C}{X}-\inprod{C}{\xsol}}+\twonorm{\Amap{X}-b}\right).
\]
If in addition, strict complementary holds for some primal dual solution pair $(\xsol,\ysol)\in 
\xsolset\times \ysolset$, then $\zeta_1=\zeta_2=2$ regardless of $\epsilon$. 
\end{lemma}
\begin{proof}
	The result in \cite[Section 4]{sturm2000error} requires the set  $S_1=\{y\mid F(y)\leq \inprod{-b}{\ysol}+ \epsilon\}$, and the set  $S_2=\{X\mid X\succeq 0,
	|\inprod{C}{X}- \inprod{C}{\xsol}|\leq \epsilon,\,\text{and}\, \twonorm{\Amap{X}-b}\leq \epsilon\}$ being compact. 
	Using \cite[Theorem 7.21]{ruszczynski2006nonlinear}, the optimization problem 


	$\min_{X\succeq 0} g(X):= \inprod{C}{X} +\gamma \twonorm{\Amap X-b}$


	has the same solution set as 
	the primal SDP \eqref{p} for some large $\gamma >0$ . Thus the compactness of the set $S_1$ and $S_2$ is ensured by Lemma \ref{lem: sublevelset} below, and the proof is completed.
\end{proof}

\begin{lemma}[Compact sublevel set]\label{lem: sublevelset}
	If a convex lower semicontinuous function $f(x): \RR^\dm \rightarrow \RR\cup \{\infty\}$ has a compact nonempty solution set, then all of its sublevel set is compact. 
\end{lemma}
\begin{proof}
	Suppose for some $L\in \RR$, then the closed sublevel set $S_L = \{x\mid f(x)\leq L\}$ is unbounded. Then there is a unit direction vector $\gamma \in \RR$ such that for all $x \in S_L,\alpha \geq 0$, $x+\alpha \gamma \in S_L$.  This in particular violates the fact the solution set is bounded and the proof is completed. 
\end{proof}


\subsection{Relating primal and dual convergence} \label{sec: Relating primal and dual convergence}
In this section, we present the convergence measures of the iterates 
$X_t$ in \eqref{eq: X_tdefinition} and dual iterates $y_t$ in terms of the convergence of $F(y_t)$. Specifically, 
we consider primal feasibility $X\succeq 0$ and $\Amap X=b$, dual feasibility 
$Z(y)\succeq 0$, and primal-dual optimality $|\inprod{C}{X}-\inprod{b}{y}|$.  
Such reduction allows us to focus on the dual penalization objective $F(y_t)$ in the next two sections. 

To facilitate the presentation, we denote 
$
D_{\xsolset} :\,=\sup_{\xsol\in\xsolset} \nucnorm{\xsol},
$
and 
$
D_y \geq \sup_{F(y)\leq F(y_0)} \norm{y}.
$
Also we introduce the shorthand $\bar{F}_t:\, = \bar{F}_{(V_t,\bar{X}_t)}$. 
Let us first deals with primal feasibility. 
\begin{lemma}[Primal Feasibility]\label{lem:primal-feas}
	At every descent step $t$, we have approximate primal feasibility
	$$ X_{t}\succeq 0, \quad\text{and}\quad \|b - \Amap X_{t}\|^2\leq \frac{2\rho_t}{\beta}(F(y_{t}) - F(\ysol )). $$
\end{lemma}
\begin{proof}
	The $X_{t+1} =\eta^\star_{t}\bar{X}_t+V_t S^\star _{t}V_t^\top$ in \eqref{eq: X_tdefinition}.
	Since $\eta\geq 0$ and $S^\star_t\succeq 0$ by construction in \eqref{eq: subproblemAMSpectralbundleMethod}, $X_t$ is 
	positive semidefinite. 
   
    From the optimality of \eqref{eq: subproblemAMSpectralbundleMethod}, and the descent step assumption, we reach that 
    \begin{equation}\label{eq: linearFeasibilityEq1}
    -b + \Amap(X_t)  = \rho_t (y_{t}-y_{t+1}).
    \end{equation}
    Hence $\twonorm{-b + \Amap(X_t)}^2  = \rho_t^2\twonorm{ (y_{t}-y_{t+1})}^2$. The distance traveled during any descent step can be bounded by the objective value gap as
	$$ \frac{\rho_t}{2}\|y_{t+1} - y_t\|^2 \leq F(y_{t})-\bar F_{t}(y_{t+1}) \leq \frac{F(y_{t})-F(y_{t+1})}{\beta} \leq \frac{F(y_{t})-F(\ysol)}{\beta}$$
	where the first inequality uses the fact that $z_{t+1}$ minimizes $\widetilde F^{t}(\cdot)+\frac{\rho_t}{2}\|\cdot-y_t\|$ and the second inequality uses the definition of a descent step.

	Combining this with our feasibility bound completes our proof.
\end{proof}

Next, we consider dual feasibility. 
\begin{lemma}[Dual Feasibility]\label{lem:dual-feas}
	At every descent step $t$, with the choice $\alpha \geq 2 D_{\xsolset}$, we have approximate dual feasibility
	$$ \lambda_{\min}(C - \Ajmap y_{t+1})\geq \frac{-(F(y_{t}) - F(\ysol))}{D_{\xsolset}}.$$
\end{lemma}
\begin{proof}
	Standard strong duality $(\inprod{C}{\xsol}=\inprod{b}{\ysol}\implies \inprod{\xsol}{Z(\ysol)}=0)$ 
    show for any $\xsol\in\xsolset$,
	\begin{align*}
	\langle b, y_{t+1}-\ysol \rangle & = \langle \Amap\xsol, y_{t+1}-\ysol \rangle\\
	& \leq \langle \xsol, \Ajmap(y_{t+1}-\ysol)\rangle = \langle \xsol, Z(\ysol)-Z(y_{t+1})\rangle\\
	& \leq -\nucnorm{\xsol}\min\{\lambda_{\min}(C - \Ajmap y_{t+1}),0\}.
	\end{align*}
	Recalling our assumption that $\alpha \geq 2D \geq 2\nucnorm{\xsol}$ and $\ysol$ indeed solves \eqref{eq: penaltySDP}, 
	yield the claimed feasibility bound
	\begin{align*}
	F(y_{t}) - F(\ysol) \geq F(y_{t+1}) - F(\ysol) & = \langle -b, y_{t+1}-\ysol \rangle -\alpha\min\{\lambda_{\min}(C - \Ajmap y_{t+1}),0\}\\
	& \geq -\|\xsol\|\min\{\lambda_{\min}(C - \Ajmap y_{t+1}),0\}. \qedhere
	\end{align*}
\end{proof}

Finally, we consider the primal-dual optimality.
\begin{lemma}[Primal-Dual Optimality]\label{lem:optimality}
	At every descent step $t$,  with the choice $\alpha \geq 2 D_{\xsolset}$, we have approximate primal-dual optimality bounded above by
	$$\langle b, y_{t+1} \rangle - \langle C, X_t\rangle \leq \frac{\alpha}{D_{\xsolset}}(F(y_{t}) - F(\ysol)) + \sqrt{\frac{2\rho_t}{\beta}(F(y_{t}) - F(\ysol))}\ D_y.$$
	and below by
	$$
	\langle b, y_{t+1} \rangle - \langle C, X_t\rangle \geq -\frac{1-\beta}{\beta}(F(y_{t}) - F(\ysol)) - \sqrt{\frac{2\rho_t}{\beta}(F(y_{t}) - F(\ysol))}\ D_y.$$
\end{lemma}
\begin{proof}
	The standard duality analysis shows the primal-dual objective gap equals
	\begin{align*}
	\langle b, y_{t+1} \rangle - \langle C, X_t\rangle &= \langle \Amap X_t, y_{t+1} \rangle - \langle C, X_t\rangle + \langle b-\Amap X_t, y_{t+1}\rangle\\
	&= \langle X_t, \Ajmap y_{t+1} - C\rangle + \langle b-\Amap X_t, y_{t+1}\rangle.
	\end{align*}
	Notice that the second term here is bounded above and below as
	$$ |\langle b-\Amap X_t, y_{t+1}\rangle| \leq \sqrt{\frac{2\rho_t}{\beta}(F(y_{t}) - F(\ysol))}\ \|y_t\|\leq \sqrt{\frac{2\rho_t}{\beta}(F(y_{t}) - F(\ysol))}\ D_y$$
	by Lemma~\ref{lem:primal-feas}. 
	Hence we only need to show that the first term also approaches zero (that is, we approach holding complementary slackness).
	
	An upper bound on this inner product follows from Lemma~\ref{lem:dual-feas} as
	$$ 
	\langle X_t, \Ajmap y_{t+1} - C \rangle \leq -\nucnorm{X_t}\lambda_{\min}(C - \Ajmap y_{t+1}) \leq \frac{\nucnorm{X_t}(F(y_{t})-F(\ysol))}{D_{\xsolset}}.
	$$
	Combining the above with $\tr(X_t) \leq \alpha$ by construction, we have  
	$$
	\langle b, y_{t+1} \rangle - \langle C, X_t\rangle \leq \frac{\alpha}{D_{\xsolset}}(F(y_{t}) - F(\ysol)) + \sqrt{\frac{2\rho_t}{\beta}(F(y_{t}) - F(\ysol))}\ D_y.
	$$

	A lower bound on this inner product follows as

	\begin{align*}
	\frac{1-\beta}{\beta}(F(y_{t}) - F(y_{t+1})) 
	&\overset{(a)}{\geq} F(y_{t+1}) - \bar F_{t}(y_{t+1})\\
	& \overset{(b)}{=} -\alpha\min\{\lambda_{\min}(C-\Ajmap y_{t+1}),0\} + \langle C, X_{t}\rangle -\langle \Amap X_{t}, y_{t+1}\rangle\\
	& \geq \langle X_{t}, C - \Ajmap y_{t+1}\rangle,
	\end{align*}
	where step $(a)$ follows from the definition of a descent step, and steo $(b)$ follows 
	from the definition of $\bar{F}_t$ and the optimality of $z_{t+1}= y_{t+1}$ and $X_t$ in 
	\eqref{eq: subproblemAMSpectralbundleMethod}. Hence
	\begin{align*}
	\langle b, y_{t+1} \rangle - \langle C, X_t\rangle &\geq -\frac{1-\beta}{\beta}(F(y_{t}) - F(\ysol)) - \sqrt{\frac{2\rho_t}{\beta}(F(y_{t}) - F(\ysol))}\ D_y.\qedhere
	\end{align*}

\end{proof}
\paragraph{Conclusion of primal dual convergence relattion} 
In conclusion, for an $y_t$ from the descent step satisfying $F(y_t)-F(\ysol)\leq \epsilon$,
assuming $\rho$, $\beta$, and $y_0$ are 
$\epsilon$ independent, then both primal feasibility (in terms of $\norm{b-\Amap X_t}$)
and primal-dual optimality (in terms of $\langle b, y_{t+1} \rangle - \langle C, X_t\rangle$) is $\bigO(\sqrt{\epsilon})$, while dual feasibility 
(in terms of $\lambda_{\min}(C - \Ajmap y_{t+1})$) is $\bigO(\epsilon)$.

\subsection{Proof of Theorem \ref{thm: sublinearates}}\label{sec: SublinearRates}
We shall utilize results in \cite{Du2017,grimmer2019general,kiwiel2000efficiency}. These results
state that a proximal bundle method admits a $y_t$ from descent step satisfying $F(y_t)-F(\ysol)\leq \epsilon$
within $\bigO(\frac{\log (\frac{1}{\epsilon})}{\epsilon^3})$ 
steps for general $F$, and  $\bigO(\frac{\log (\frac{1}{\epsilon})}{\epsilon})$
steps when $F$ has quadratic growth, a consequence of strict complementarity as shown 
in Lemma \ref{lem: qg}. Once we establish these convergence rates for Block-Spec, and 
HR-Spec, then Lemma \ref{lem:primal-feas}, \ref{lem:dual-feas}, and \ref{lem:optimality} imply the 
results for primal feasibility, dual feasibility, and primal-dual optimality. 

Recall 
in Section \ref{sec: Block spectral bundle method}, we claimed
the following inequalities \eqref{eq: subgradientAndAggregation},
\begin{equation*}
\bar{F}_{\texttt{simple}}(y)\overset{(a)}{\leq} \bar{F}_{(V_{t+1},\bar{X}_{t+1})}(y) \overset{(b)}{\leq} F(y).
\end{equation*}
Here the inequality $(b)$ is satisfied by construction. 
Once we establish the inequality $(a)$, we can immediately employ the results in \cite{Du2017,grimmer2019general}
\footnote{Both proofs actually only show convergence for a proximal bundle method with cut aggregation 
$\bar{F}_{\texttt{simple}}$ or the method with a multiple yet \emph{finite} cuts. As 
remarked by the authors of \cite[Section 3]{Du2017}, `` All the results hold true for
both versions, because the analysis of the method with multiple cuts uses the version
with cut aggregation anyway", a detailed check yields that \eqref{eq: subgradientAndAggregation} 
(which appears as the inequality (12) in \cite{Du2017}) is the 
only ingredient needed 
to prove the convergence rate while the finiteness of the cuts is not needed. A future work \cite{DG2020} of the second author further confirms this fact.}. To avoid extra complication,
we only show the inequalities for Block-Spec and defer the derivation for HR-Spec to 
Section \ref{sec: Model comparison for HR bundle method} in the appendix. 
Recall we write $\bar{F}_{t+1}:\,=\bar{F}_{(V_{t+1},\bar{X}_{t+1})}$ for short.

Let us first define $\bar{F}_{\texttt{simple}}$. Denote $g_t \in \partial F(z_{t+1})$ and $s_{t+1}= -\rho_t (z_{t+1}-y_t)$. The simplest model 
$\bar{F}_{\texttt{simple}}$ (for bundle method with cut aggregation)
in \cite[Section 2.2]{Du2017} for the $t+1$ step, which is a maximum of two lower bounds based on the aggregation $s_{t+1}$ and a subgradient $g_t$, takes form  
\begin{equation}\label{eq: simplemodel}
\bar{F}_{\texttt{simple}}(y) = \max\{F(z_{t+1})+\inprod{g_{t+1}}{y-z_{t+1}}, F_{t}(z_{t+1})+\inprod{s_{t+1}}{y-z_{t+1}}\}.
\end{equation}
By the optimality condition of \eqref{eq: subproblemAMSpectralbundleMethod}, and definition of $X_t,z_{t+1}$, we know that 
\begin{align}
-b + \Amap(X_t) & = \rho_t (y_t-z_{t+1})=s_{t+1} \label{eq: optimalityConditionEqsubproblemAMSpectralbundleMethod}\\ 
\bar{F}_{t}(z_{t+1}) & = \inprod{-b}{z_{t+1}} +\inprod{X_t}{\Amap^*(z_{t+1})-C}.\label{eq: F_{t+1}(z_{t+1})} 
\end{align}
Let the $v_+$ be the top eignevector of $\lambda_{\max}(\Amap^* z_{t+1}-C)$ if $\lambda_{\max}(\Amap^* z_{t+1}-C)>0$ 
and be zero otherwise. Using \eqref{eq: optimalityConditionEqsubproblemAMSpectralbundleMethod}, \eqref{eq: F_{t+1}(z_{t+1})},
and the definition of $v_+$, the simple model  $\bar{F}_{\texttt{simple}}$ can be rewritten as 
\begin{equation}\label{eq: simplemodelLambdaform}
\begin{aligned} 
\bar{F}_{\texttt{simple}}(y)
& = \max\{ \inprod{-b}{y}+ \inprod{\alpha v_+v_+^\top}{\Amap^*y-C},  \inprod{-b}{y}+ \inprod{X_t}{\Amap^*y-C} \}. 
\end{aligned}
\end{equation}
Recall the $t$-th spectral set $\mathcal{W}_{t+1}$ in \eqref{eq: tthspectralset} is defined as
\[
\mathcal{W}_{t+1}:=\{\eta \bar{X}_{t+1} + V_{t+1}SV_{t+1}^\top \mid \eta \geq 0,\;S\in \symMat_+^{r_{t+1}},\;\text{and}\;\eta\alpha  + \tr(S)\leq \alpha \},
\]
which is the constraint set of Problem \eqref{eq: subproblemAMSpectralbundleMethod}.  We find that 
${X}_t=\bar{X}_{t+1}$ (due to the update of Block-Spec) and the subgradient $\alpha v_+v_+^\top \in  \partial \max\{\lambda_{\max}(\Amap^* (z_{t+1})-C),0\}$ are both in $\mathcal{W}_{t+1}$:
\begin{equation}\label{eq: aggregatebarxinthenewsetAmnesiaMethod}
\bar{X}_{t+1}={X}_t\in \mathcal{W}_{t+1}\quad\text{and}\quad \alpha  v_1v_1^\top \in \mathcal{W}_{t+1}.
\end{equation}
Hence \eqref{eq: subgradientAndAggregation} is established, and our proof is complete. 

\subsection{Local linear convergence of Block bundle method}\label{sec: Proof of Local linear convergence of Block bundle method}
First, we discuss estimates of $T_0$, $\rho$ and the contraction factor, and the condition $\bar{r}\geq r_d$ for local linear convergence.

Let us introduce some notations to facilitate the discussion: the gap parameter $\delta := \inf_{\ysol \in \ysolset}\max_{r\leq r_d} \lambda_{r}(-Z(\ysol ))-\lambda_{r+1}(-Z(\ysol))$ and
the quadratic growth 
parameter $\gamma>0$ from Lemma \ref{lem: qg}.
The gap parameter $\delta$ is nonzero 
from the definition of $r_d$, the compactness of $\ysolset$, and continuity of the 
function $\max_{r\leq r_d} \lambda_{r}(Z(\cdot ))-\lambda_{r+1}(Z(\cdot))$. When the dual solution is unique, we have
$\delta =  \lambda_{r_d}(-Z(\ysol))-\lambda_{r_d+1}(-Z(\ysol))$. 

\paragraph{Estimate of $T_0$, $\rho$, and contraction factor} The number $T_0$ is chosen so that $Z(y_t)$ is $\frac{\delta}{3}$ close to the solution set 
$Z(\ysolset)=\{Z(\ysol)\mid \ysol \in \ysolset\}$. Using Theorem \ref{thm: sublinearates}, we know 
$T_0$ satisfies 
\begin{equation} \label{eq: T0estimate}
T_0\leq \kappa \frac{\log\left(\frac{\delta^2\gamma}{9\opnorm{\Amap^*}^2}\right)}{\frac{\delta^2\gamma}{9\opnorm{\Amap^*}^2}}.
\end{equation} 
In our proof in Section \ref{sec: Proof of the quadratic accurate model}, the condition on regularization parameter $\rho$ 
is that
\begin{equation}\label{eq: rhoEstimate}
\rho \geq 4\opnorm{\Amap^*}^2\max\{\frac{72\sup_{\ysol \in \ysolset}\opnorm{2Z(\ysol)}}{\delta^2},\frac{9(8\sqrt{2}+16)}{\delta}\}.
\end{equation}
The contraction quantity is the distance to the dual solution set, $\dist(y_t,\ysolset)$, and the factor is $\sqrt{\frac{\rho}{2\gamma +\rho}}$, see 
\eqref{eq: approximationModelStep4} in Section \ref{sec: Linear convergence under a quadratic accurate model}.

\paragraph{The condition $\bar{r}\geq r_d$} Modern applications \cite{candes2009exact,candes2013phaselift,recht2010guaranteed,ding2020regularity} of \eqref{p} actually shows that $\xsol$ is unique, admits rank $\rsol :\,=\rank(\xsol)\ll \dm$, and satisfies strict complementarity under certain structural probabilistic assumptions. If in addition, dual solution is unique, then we only need $\bar{r}\geq r_d=\rsol$. Hence, the subproblem \eqref{eq: subproblemSpectralbundleMethod} can be solved efficiently for these problems. Also $\bar{r}\geq \rsol$ should be required even from an eigenvalue computational perspective as the bottom $\rsol$ eigenvalues of the slack $Z(y_t)$ start to coalesce once $y_t$ is close to $\ysolset$. Moreover, as we observe in numerics in Section \ref{sec: numerics}, even if there are multiple dual solutions, $\bar{r}\geq \rsol $ is still enough to guarantee quick convergence while $\bar{r} <\rsol$ introduces slow convergence.

\paragraph{Proof outline} The proof is split into two steps. In Section \ref{sec: Linear convergence under a quadratic accurate model}, 
We first show  that linear convergence 
can be proved under the claim that 
the model $\bar{F}_{(V_t,\bar{X}_t)}(z)$ is quadratic accurate, 
see \eqref{eq: quadraticAccurateModel} for the definition. 
Next, in Section \ref{sec: Proof of the quadratic accurate model}, 
we show the claim indeed holds under strong duality, strict complementarity, and $\bar{r}\geq r_d$.

\subsubsection{Linear convergence under a quadratic accurate model}\label{sec: Linear convergence under a quadratic accurate model}
Let us first clarity the notion quadratic accurate model. Our model $\bar{F}_{t}(z):=\bar{F}_{(V_t,\bar{X}_t)}(z)$ is quadratic 
accurate if there is some $\rho>0$ (independent of $t$) such that for any $z\in \RR^{\ncons}$,
\begin{equation}\label{eq: quadraticAccurateModel}
\begin{aligned}
\bar{F}_{t}(z)\leq F(z)\leq \bar{F}_{t}(z) + \frac{\rho}{2}\twonorm{z-y_t}^2.
\end{aligned}
\end{equation}  
Under this condition the method with $\rho_t =\rho$ always take a descent step for $\beta\leq \frac{1}{2}$. We verify the model is quadratic accurate, 
i.e., \eqref{eq: quadraticAccurateModel} holds, in the next section.

Now suppose \eqref{eq: quadraticAccurateModel} is satisfied. We know the minimizer 
$z_{t}^\star$ of $\bar{F}_{t}(z) + \frac{\rho}{2}\twonorm{z-y_t}^2$ satisfies that for any $z\in \RR^\ncons$
\begin{equation}\label{eq: threePointInequality}
\begin{aligned} 
\bar{F}_t(z_t^\star) +\frac{\rho}{2}\twonorm{z_t^\star-y_t}^2 + \frac{\rho}{2}\twonorm{z^\star_t-z}^2\leq \bar{F}_t(z) + \frac{\rho}{2}\twonorm{z- y_t}^2,
\end{aligned} 
\end{equation} 
as $\bar{F}_t(z)+\frac{\rho}{2}\twonorm{z-y_t}^2$ is $\rho$-strongly convex. By setting $z=y_t$ and \eqref{eq: quadraticAccurateModel}, we find that 
\begin{equation}\label{eq: approximationModelStep1}
F(y_t)-\bar{F}_t(z_t^\star)\geq \rho \twonorm{z_t^\star-y_t}^2\geq 0,\quad \text{and}\quad \frac{\rho}{2}\twonorm{z_t^\star-y_t}^2 \leq F(y_t)-F(z_{t}^\star ).
\end{equation}
Next, using \eqref{eq: quadraticAccurateModel} in the following step $(a)$ and 
\eqref{eq: approximationModelStep1} in the following step $(b)$, we find that
\begin{equation}\label{eq: descentStepASBMLocalLinearConvergence}
\begin{aligned}
\beta \left(F(y_t)-\bar{F}_t(z_t^\star)\right)\leq  \frac{1}{2}\left(F(y_t)-\bar{F}_t(z_t^\star)\right) 
& \overset{(a)}{\leq}\frac{1}{2}\left(F(y_t)-F(z_t^\star)\right)+\frac{\rho}{4}\twonorm{y_t-z_t^\star}^2\\
&\overset{(b)}{\leq}  F(y_t)-F(z_{t}^\star).
\end{aligned} 
\end{equation}
Hence, we see the method will indeed take a descent step and $z_t^\star$ is in the sublevel set defined by $\{y\mid F(y)\leq F(y_0)\}$. 


By setting $z=\ysol$ for an arbitrary $\ysol\in \ysolset$ in \eqref{eq: threePointInequality}, and using 
\eqref{eq: quadraticAccurateModel}, we find that 
\begin{equation}\label{eq: approximationModelStep2}
\begin{aligned} 
F(z^\star_t)\leq F(\ysol) +\frac{\rho}{2}\left(\twonorm{\ysol -y_t}^2-\twonorm{z_t^\star-\ysol}^2\right). 
\end{aligned} 
\end{equation} 
Now recall the quadratic growth of $F$  (derived from Lemma \ref{lem: qg}) that there is a $\gamma>0$ such that for all $z\in \{y\mid F(y)\leq F(y_0)\}$, 
\[
F(z)-F(\ysol)\geq \gamma \disttwonorm^2 (z,\ysolset).
\]
Hence combining this with \eqref{eq: approximationModelStep2}, we find that 
\begin{equation}\label{eq: approximationModelStep3}
\begin{aligned} 
            \gamma \disttwonorm^2 (z_t^\star,\ysolset) &\leq \frac{\rho}{2}\left(\twonorm{\ysol -y_t}^2-\twonorm{z_t^\star-\ysol}^2\right) \\ 
 \implies  \left(\gamma +\frac{\rho}{2}\right)\disttwonorm^2 (z_t^\star,\ysolset)  & \leq \gamma \disttwonorm ^2(z^\star_t,\ysolset) + \frac{\rho}{2}\twonorm{z_t^\star-\ysol}^2\leq \frac{\rho}{2}\twonorm{\ysol -y_t}^2 \\ 
 \implies \disttwonorm^2 (z^\star_t, \ysolset) & \leq \frac{\rho}{2\gamma +\rho} \twonorm{\ysol-y_t}^2. 
\end{aligned} 
\end{equation}
Now take $\ysol$ to be the point in $\ysolset$ nearest to $y_t$, we find that
\begin{equation}\label{eq: approximationModelStep4}
\disttwonorm (z_t^\star,\ysolset)  \leq \sqrt{\frac{\rho}{2\gamma +\rho}} \disttwonorm (y_t,\ysolset). 
\end{equation}
This shows that the new iterate $y_{t+1}=z_{t+1}=z_t^\star$ (the method takes a descent step as we just argued in \eqref{eq: descentStepASBMLocalLinearConvergence}) is getting closer to the solution set $\ysolset$ geometrically with a factor of $\sqrt{\frac{\rho}{2\gamma +\rho}} <1$.  Next we prove that the model assumption \eqref{eq: quadraticAccurateModel} indeed holds. 
\subsubsection{Proof of the quadratic accurate model} \label{sec: Proof of the quadratic accurate model}
Withoutloss of generality, we may suppose our current $z_t$ is $y_t$. Let us first define the $r$-th spectral plus set of a matrix $X\in\symMat^{\dm}$ with $\lambda_{r}(X)-\lambda_{r+1}(X)>0$ as  $\faceplus r(X):=\{VSV^{\top}\mid\tr(S)\leq1,S\succeq0,S\in\symMat^{r}\}.$ Here $V\in \RR^{\dm \times r}$ is the matrix formed by orthonomal eigenvectors of $X$ corresponding to its 
$r$ largest eigenvalues. We now utilize the following lemma, proved in Section \ref{sec: proofOfImportantLemma} in the appendix, to establish that our model $\bar{F}_{t}$ is quadratic accurate.
\begin{lemma}\label{lem: importantLemmaQuadraticAccurateModel}
	Suppose $X\in\symMat^{\dm}$ has eigenvalues $\lambda_{r}(X)-\lambda_{r+1}(X)=\delta$ and denote the $\Lambda_{r,\dm}=\max\{|\lambda_{r+1}(X)|,|\lambda_\dm (X)| \}$. Then for any $Y\in \symMat^{\dm}$, the quantity $f_X(Y):\,=\max\{\lambda_1(Y),0\}-\max_{W\in \faceplus r (X)}\inprod{W}{Y}$ satisfies that 
	\begin{equation}
	\begin{aligned} \label{eq: quadraticAccurateModelLemma}
	0\leq f_X(Y) \leq \frac{8\fronorm{Y-X}^{2}\Lambda_{r,\dm}}{\delta^{2}}+\frac{(8\sqrt{2}+16)\fronorm{Y-X}^{2}}{\delta}.
	\end{aligned} 
	\end{equation}
\end{lemma}

\paragraph{Unique solution case} Let us first suppose the dual solution $\ysol$ is unique. In this case, if $y_t$ is close to $\ysol$ so that $\opnorm{Z(y_t)-\zsol}\leq \frac{\delta}{3}$ where $\delta$ is the $r_d$-th eigengap of $-\zsol$, $\delta = \lambda_{r_d}(-\zsol)-\lambda_{r_d+1}(-\zsol)>0$. Then from
Weyl's inequality, we know the $r_d-$th eigengap of $-Z(y_t)$, $\lambda_{r_d}(-Z(y_t))-\lambda_{r_d+1}(-Z(y_t))$, is at least $\frac{\delta}{3}$, and $\opnorm{Z(y_t)}\leq 
2\opnorm{\zsol}$. Denote $V\in \RR^{\dm \times r_d}$ to be the matrix formed 
by the eigenvectors corresponding to the $r_d$ largest eigenvalue of $-Z(\ysol)$. Hence we find that 
\begin{equation}
\begin{aligned} 
F(y)-\bar{F}_t(y)& =\alpha \max\{\lambda_{\max}(-Z(y)),0\} -  \max_{\eta \alpha +\tr(S)\leq \alpha , \eta \geq 0, S\in \symMat_+^{\bar{r}}}\inprod{\eta\bar{X}+V_tSV_t^\top}{-Z(y)} \\ 
&\overset{(a)}{\leq } \alpha \max\{\lambda_{\max}(-Z(y)),0\} -  \max_{\tr(S)\leq \alpha, S\in \symMat_+^{r_d}}\inprod{VSV^\top}{-Z(y)} \\
& = \alpha \left( \max\{\lambda_{\max}(-Z(y)),0\} - \alpha \max_{W\in \faceplus {r_d} (-Z(y_t))}\inprod{W}{-Z(y)} \right) \\
& \overset{(b)}{\leq }  \frac{72\fronorm{Z(y_t)-Z(y)}^{2}\opnorm{2\zsol}}{\delta^{2}}+\frac{9(8\sqrt{2}+16)\fronorm{Z(y_t)-Z(y)}^{2}}{\delta}\\
&\leq 2\opnorm{\Amap^*}^2\max\{\frac{72\opnorm{2\zsol}}{\delta^2},\frac{9(8\sqrt{2}+16)}{\delta}\}\twonorm{y-y_t}^2,
\end{aligned} 
\end{equation}
Here $(a)$ is because $\bar{r} \geq r_d$ by assumption and $(b)$ is due to Lemma \ref{lem: importantLemmaQuadraticAccurateModel}.
Combining the fact that $F(y)-\bar{F}_t(y)$ as $\bar{F}_t$ serves as a lower bound from the construction, we see the model $\bar{F}_t$ is indeed quadratic accurate with $\rho =  4\opnorm{\Amap^*}^2\max\{\frac{72\opnorm{2\zsol}}{\delta^2},\frac{9(8\sqrt{2}+16)}{\delta}\}$ in \eqref{eq: quadraticAccurateModel}.

\paragraph{Multiple dual solutions case} What if $\mathcal{Y}_\star$ has multiple points? Recall we define $\delta$ as
\begin{equation}
\begin{aligned} \label{eq: definitionOfdelta}
\delta = \inf_{\ysol \in \ysolset}\max_{r\leq r_d} \lambda_{r}(-Z(\ysol))-\lambda_{r+1}(-Z(\ysol)).
\end{aligned} 
\end{equation}
The gap $\delta$ is nonzero from the definition of $r_d$, the compactness of $\ysolset$, and continuity of the 
function $\max_{r\leq r_d} \lambda_{r}(\zsol(\cdot))-\lambda_{r+1}(\zsol(\cdot))$. Hence if $\dist(Z(y_t),Z(\ysolset))$ is 
less than a third of $\delta$, then there is an $r$ and $\ysol\in \ysolset$, such that $-Z(\ysol)$ is no more than 
$\frac{\delta}{3}$ away from $Z(y_t)$, and has $\lambda_{r}(-\zsol(y))-\lambda_{r+1}(-\zsol(y))\geq \frac{\delta}{3}$. 
Hence we can repeat previous argument for the case of unique dual solution and replace $r_d$ and $\opnorm{\zsol}$ by $r$ and $2\sup_{\ysol \in \ysolset}\opnorm{Z(\ysol)}$ respectively. Hence we 
see that the model $\bar{F}_t$ is indeed quadratic accurate in \eqref{eq: quadraticAccurateModel} with 
$
\rho =4\opnorm{\Amap^*}^2\max\{\frac{72\sup_{\ysol \in \ysolset}\opnorm{2Z(\ysol)}}{\delta^2},\frac{9(8\sqrt{2}+16)}{\delta}\}
$ 
as stated in \eqref{eq: rhoEstimate}.



  


 
 
  


\section{Defining Spectral Bundle Methods} \label{sec:def}
In this section, we lay out the main ideas of spectral bundle methods in Section
\ref{sec: Ingredients of Spectral bundle method}. 
In Section \ref{sec: SpectralbundleMethodExactDefinition}, we 
first present our Block spectral bundle method, which we believe is notationally simpler, and then the HR 
spectral bundle method. We compare their differences in Section \ref{sec: HR spectral bundle method},
and conclude the section by a discussion in solving a common small-scale subproblem of Block-Spec and HR-Spec 
in Section \ref{sec: ImportantsubproblemSolver}. 
\subsection{Ingredients of Spectral bundle method}\label{sec: Ingredients of Spectral bundle method}
In this section, we describe three main ingredients of spectral bundle method: spectral lower bound, regularization, 
and aggregation. 
\subsubsection{Subgradient lower bound and spectral lower bound}\label{sec: Subgradient lower bound and spectral lower bound}
Let us first recall the subgradient lower bound of bundle methods ~\cite[Section 7.4]{ruszczynski2011nonlinear}. 

\paragraph{Subgradient lower bound} For any convex function $F$, the set of lower minorants at a point $y$ defines the subdifferential 
$\partial F(y) = \{g\in\RR^\ncons \mid F(y') \geq F(y)+\langle g, y'-y\rangle  \text{ for all } y' \}$. Each vector $g$ is referred to as a subgradient.
The bundle method builds a model function $\bar{F}(y)$ of $F$ via a maximum of lower minorants of $F$ at certain points $\mathcal{Y}_I = \{y_i\}_{i=1}^{|I|}$:
\[
\bar{F}_{\mathcal{Y}_I}(y) :\,= \max_{i=1,\dots,|I|}\left( F(y_i)+\inprod{g_i}{y-y_i}\right), \quad g\in\partial F(y_i)
\]

If we now substitute the special form of $F$ in \eqref{eq: penaltySDP}, we find that $\bar{F}_{\mathcal{Y}_I}(y)$ can be rewritten as 
\begin{equation}\label{eq: barFfiniteeigenvectors}
\bar{F}_{\mathcal{Y}_I}(y) :\,= \max_{i=1,\dots,|I|}\inprod{-b}{y}+\inprod{\alpha v_iv_i^\top}{\Amap^*y-C} 
\end{equation}
where each $v_i,\;i=1,\dots,|I|$ is a top eigenvector of $\Amap^*y_i-C$ if $\Amap^*y_i-C$ has a positive eigenvalue and is $0$ otherwise.

\paragraph{Spectral lower bound} The key idea of spectral bundle method is maximizing over \emph{infinitely} many lower minorants indexed by an 
SDP representable set,
 which 
might be called spectral lower bound. To state the idea more precisely, first note the fact that for any $Z\in \symMat^\dm$:
\[
\max \{\lambda_{\max}(Z),0\} = \max_{\inprod{X}{I}\leq 1,X\succeq 0}\inprod{X}{Z}.
\]
Hence we may rewrite $F$ as 
\begin{equation}\label{eqn: dualpenaltyobjectiverewrite}
F(y) = \max_{\inprod{X}{I}\leq \alpha , X\succeq 0} \inprod{-b}{y} +\inprod{X}{\Amap^*y-C}.
\end{equation}
Of course, this form is no easier to solve than the original penalized form \eqref{eq: penaltySDP}. However, we can replace the constraint set $\{\inprod{X}{I}\leq \alpha , X\succeq 0\}$ by a smaller convex set. One choice is that we compute a matrix $V\in \RR^{n\times r}$ for some small value $r$ with orthonormal columns, i.e., $V^\top V= I\in \RR^{r\times r}$.
Then we form the model based on $V$:
 \begin{equation}\label{eqn: modelFoverV}
 \bar{F}_{V}(y) :\,= \max_{\inprod{S}{I}\leq \alpha , S\in \symMat_+^r} \inprod{-b}{y} +\inprod{VSV^\top}{\Amap^*y-C}.
 \end{equation}
This function $\bar{F}_{V}$ should serve as a better approximation of $F$ than $\bar{F}_{\mathcal{Y}(I)}$ in \eqref{eqn: dualpenaltyobjectiverewrite} 
so long as the column span of $V$ contains the span of  $v_1$ to $v_{|I|}$.

How to choose $V$? The most natural choice 
is to let $V$ be the $r$ orthonormal eigenvectors of $\Amap^*y-C$ for our current iterate $y$. This is in fact the main idea of the block spectral bundle method. 
However, the choice of $V$ for the HR spectral bundle method is more complicated. We detail the procedures in choosing $V$ in Section
\ref{sec: Block spectral bundle method} for Block-Spec
and 
Section \ref{sec: HR spectral bundle method} for HR-Spec.
 \subsubsection{Aggregation and Regularization} 
Having explained the key idea of spectral bundle method, we now explain the two ingredients of proximal bundle method \cite{Du2017} 
that is also effective for spectral bundle method: 
regularization and aggregation. 

\paragraph{Regularization} Directly minimize the model $\bar{F}_V$ in \eqref{eqn: modelFoverV} in each iteration would correspond to a cutting plane method, which are known to be unstable. Hence the approach used by bundle methods is to add a regularization term $\frac{\rho}{2}\norm{y-\hat{y}^2}$ to \eqref{eqn: modelFoverV}, where $\hat{y}$ is a reference point and $\rho>0$ is a regularization parameter:
\[
\tilde{F}_{V,\hat{y},\rho}(y) := \bar{F}_V(y) + \frac{\rho}{2}\norm{y-\hat{y}}^2.
\] 
Here, the regularization $\rho$ can be considered as a step size parameter. The larger the $\rho$, the smaller the step size. 

\paragraph{Aggregation} We now describe the final ingredient: aggregation. While running a bundle method, rather than storing all the past subgradients $\alpha v_iv_i^T$ for the model $\bar F$, we can collect a weighted average of these top eigenvectors from past iterates as a matrix $\bar{X}$ such that 
\[
\bar{X}\in \symMat_+^\dm, \quad \text{and}\quad \nucnorm{\bar{X}}=\tr(X)\leq \alpha.
\]
We then build model using this matrix $\bar{X}$ along with a matrix $V\in \RR^{\dm \times r}$ with orthonormal columns: 
\begin{equation}\label{eq: modelaggreagte}
\bar{F}_{(V,\bar{X})}(y) :\,= \max_{\eta \alpha+\tr(S)\leq \alpha , \eta \geq 0, S\in \symMat_+^r} \inprod{-b}{y} +\inprod{\eta\bar{X}+VSV^\top}{\Amap^*y-C}.
\end{equation}
The regularized version with a reference point is then 
\begin{equation}\label{eq: modelregularizedfull} 
\tilde{F}_{(V,\bar{X},\hat{y},\rho)}(y) := \bar{F}_{(V,\bar{X})}(y) + \frac{\rho}{2}\norm{y-\hat{y}}^2.
\end{equation} 
We specify the exact choice of $\bar{X}$ required by each spectral bundle method in the next section.

\subsection{The Spectral bundle Methods}\label{sec: SpectralbundleMethodExactDefinition}
Having seen the three ingredients, we define the spectral bundle methods in this section. 
We shall first display the Block-Spec due to its simpler form and then 
present the HR-Spec. 
 
\subsubsection{Block spectral bundle method}\label{sec: Block spectral bundle method} 
Block-Spec starts with two initial iterates $z_0=y_0\in \RR^\ncons$. Here $z_0$ is the always exploring iterate 
and $y_0$ is the reference iterate which does not move unless the model $\bar{F}_{V,\bar{X}}$ in \eqref{eq: modelaggreagte} is accurate enough (see step 2 below). 
We also intialize a primal iterate $\bar{X}_0=0$, an integer $\bar{r}>0$, and a matrix $V_0 \in \real^{\dm \times r}$ which is the collection of 
top $\bar{r}$ eigenvectors of $\Amap^*z_{t+1}-C$. Finally, we pick a fraction $\beta \in (0,1)$ which represents whether the 
quality of current model $\bar{F}_{V,\bar{X}}$ is satisfactory or not. Block-Spec iterates the following: 
\begin{enumerate}
	\item \textbf{Solving the regularized model:} Pick $\rho_t>0$ and solve $\min_{z} \;\tilde{F}_{(V_t,\bar{X}_t,y_t,\rho_t)} (z)$. Concretely, the minimization 
	problem is 
	\begin{equation}\label{eq: subproblemAMSpectralbundleMethod}
	\min_{z} \quad \max_{\eta \alpha +\tr(S)\leq \alpha , \eta \geq 0,S\in \symMat_+^{\bar{r}}} \inprod{-b}{z} +\inprod{\eta\bar{X}_t+V_t S V_t^\top}{\Amap^*z-C}+\frac{\rho_t}{2}\norm{z-y_t}^2. 
	\end{equation}
	Obtain $z^\star_{t}, \eta^\star_{t}$, and $S^\star_{t}$ solving the minimax problem above. Define  $z_{t+1}= z^\star_t$ and 
	\begin{equation} \label{eq: X_tdefinition}
	X_t= \eta^\star_{t}\bar{X}_t+V_t S^\star _{t}V_t^\top
	\end{equation}
	\item \textbf{Testing the model based on objective value decrease:} If $F(z_{t+1}) \leq F(y_t) - \beta\left(F(y_t) - \bar{F}_{(V_t,\bar{X}_t)}(z_{t+1})\right)$,
	then we set $y_{t+1} = z_{t+1}$. In this case, we call this step a \emph{descent step}. Otherwise, we set $y_{t+1} = y_{t}$, and 
	this step is called a \emph{null step}.
	\item \textbf{Aggregation:} We update the aggregate matrix $\bar{X}_t$ by 
	\begin{equation}\label{eq: aggregationUpdateOurMethod}
	\bar{X}_{t+1} = X_t = \eta^\star_t\bar{X}_t+V_t S^\star_t V_t^\top.
	\end{equation}
	\item \textbf{Update the bundle:} We update $V_t$ by computing the top $\bar{r}$ orthonormal eigenvectors $v_1,\dots,v_{\bar{r}}$ of $\Amap^*z_{t+1}-C$, and set 
	\begin{equation}
	V_{t+1}=[v_1,\dots ,v_{\bar{r}}].
	\end{equation}
\end{enumerate}
Notice that we did not specify how to find $\eta_t^\star,S_t^\star$ and $z^\star_t$, and the choice of $\rho_t$. We describe a way of solving the minimax subproblem 
\eqref{eq: subproblemAMSpectralbundleMethod} in Section \ref{sec: ImportantsubproblemSolver}. 
From Theorem \ref{thm: sublinearates} and \ref{thm: linear convergence of Block SBM under the extra condition strict complementarity} 
in Section \ref{sec: analysis}, we are guided to set $\rho$ to be constant.

To avoid extra computation in the evaluation of $\bar{F}_{(V_t,\bar{X}_t)}$ in step 2, 
we note that 
$
\bar{F}_{(V_t,\bar{X}_t)}(z_{t+1}) = \tilde{F}_{(V_t,\bar{X}_t,y_t,\rho_t)}(z_{t+1})-\frac{\rho_t}{2}\twonorm{z_{t+1}-y_t}^2.
$ 

Finally, we remark an important property of Block-Spec: 
the new model $\bar{F}_{(V_{t+1},\bar{X}_{t+1})}$ lower bounds $F(y)$ by construction, 
and is always better than the simple model $\bar{F}_{\texttt{simple}}$ , defined in \eqref{eq: simplemodel}, of the bundle method with cut aggregation \cite[Section 2.2]{Du2017}: 
for all $y\in \RR^{\ncons}$
\begin{equation}\label{eq: subgradientAndAggregation}
\bar{F}_{\texttt{simple}}(y)\leq \bar{F}_{(V_{t+1},\bar{X}_{t+1})}(y)\leq F(y).
\end{equation}
We prove this inequality in Section \ref{sec: SublinearRates} which is the key insight for 
establishing convergence of sublinear rates.

\subsubsection{HR spectral bundle method}\label{sec: HR spectral bundle method} The HR bundle method 
proceeds the same as the Block bundle method in the first two steps but differs at the third and fourth 
step due to extra complexity of aggregation and orthogonalization. 
\begin{enumerate}
	\item [3'.] \textbf{Aggregation:} We do an eigenvalue decomposition of $S^\star_t= Q_1\Lambda_1 Q_1^\top + Q_2 \Lambda_2 Q_2^\top$, where $\Lambda_1$ consists of the largest 
	$\bar{r}-1$ eigenvalues, and $\Lambda_2$ consists of the rest eigenvalues. We now update $\bar{X}_t$ by 
	\begin{equation}\label{eqn: updatebarX_t+1FC}
	\bar{X}_{t+1}= \eta^\star_t \bar{X}_t +V_tQ_2\Lambda_2Q_2^\top V_t^\top. 
	\end{equation}
	The rational behind this update is that we keep the important information in the matrix $V_t$ which can improve 
	the model accuracy of $\bar{F}$ in the next round.
	\item [4'.] \textbf{Update the bundle:} Next we update $V_t$. We first compute a top eigenvector $v_{t+1}$ of $\Amap^*z_{t+1}-C$. Finally, we set $V_{t+1}$ consists of an orthonormal basis of the 
	column space of $[V_tQ_1;v_{t+1}]$, e.g., using QR factorization:
	\begin{equation}\label{eqn: updateVFC}
	QR = [V_tQ_1;v_{t+1}],\quad V_{t+1}=Q.
	\end{equation}
\end{enumerate}

In step 3', we could choose $r$ eigenvectors with $0\leq r\leq  \bar{r}-1$. Such choice reduce the variable size $S$ in \eqref{eq: subproblemAMSpectralbundleMethod} to $\mathbb{S}^{r+1}$. 
The inequality \eqref{eq: subgradientAndAggregation} continues to hold for HR bundle method as well even under this choice of $r$. 
We defer the detailed derivation to Section \ref{sec: Model comparison for HR bundle method} in the appendix.

\subsubsection{A discussion on the differences of FC and Blockbundle methods}\label{sec: A discussion on the differences of FC and Blockbundle methods}
Here we discuss the difference in updating the matrix $V$ between the two spectral bundle methods. Denote the current iterate as $T>1$.
\begin{enumerate}
	\item HR-Spec computes a matrix $V$ using past top eigenvectors $\{v_t\}_{t=1}^{T}$ of iterates $\{\Amap^*y_t-C\}_{t=1}^{T}$. The orthoganalization step ensures that the trace
	 of $VSV^\top$ is exactly the same as $S$. 
	\item  Block-Spec instead uses the current iterate $z_T$ and compute $V$ as a block of $\bar{r}$ top eigenvectors  of 
	$C-\Amap^*z_T$ without referencing any past eigenvectors. Hence the name block.
\end{enumerate} 
Despite the simple form of Block-spec, the block procedure is advantageous when evaluating top eigenvectors becoming harder when $y_T$ is close to the solution set $\mathcal{Y}_\star$ as the eigenvalues of $Z(y_T)$ begin coalescing. In such setting, it is more computational efficient to compute a block of eigenvectors instead of just the top one. In fact, a block Lanczos method is actually employed in the original paper \cite{helmberg2000spectral} in computing eigenvectors, even though only one is needed. 

\paragraph{A hybrid approach} We note it is possible to combine these two ideas when forming $V$. For example, one could update the aggregation $\bar{X}$ according to \eqref{eqn: updatebarX_t+1FC}. Then 
introduce the $\bar{r}$ top eigenvectors $v_1,\dots,v_{\bar{r}}$ of $\Amap^*z-C$ and update $V$ using the following computation:
\begin{equation}\label{eq: combinedupdateMethod}
QR = [VQ_1; v_{1},\dots,v_{\bar{r}}],\quad V=Q.
\end{equation}
Here $QR$ stands for the $Q,R$ factors of QR factorization. Of course, in this case, we should be able to store $2\bar{r}$ many $\dm$-dimensional vectors. The properties described in 
\eqref{eq: subgradientAndAggregation} continue to hold. Since the model $F_{V,\bar{X}}$ is better than the model in HR and Block bundle methods, the theoretical 
results displayed in Section \ref{sec: analysis} continues to hold for this hybrid approach.


\subsubsection{Solving the minimax subproblem \eqref{eq: subproblemAMSpectralbundleMethod}}\label{sec: ImportantsubproblemSolver}
The last step in fully specifying the spectral bundle method is to discuss how to solve the minimax subproblem  \eqref{eq: subproblemAMSpectralbundleMethod}. Define the $t$-th spectral set $\mathcal{W}_t$ as 
\begin{equation}\label{eq: tthspectralset}
\mathcal{W}_t = \{\eta \bar{X}_{t} + V_{t}SV_{t}^\top \mid \eta \geq 0,\;S\in \symMat_+^{\bar{r}},\;\text{and}\;\eta \alpha + \tr(S)\leq \alpha \}.
\end{equation}
Hence we may rewrite the minimax problem  \eqref{eq: subproblemAMSpectralbundleMethod} as 
\begin{equation}\label{eq: subproblemSpectralbundleMethod}
\min_{z} \max_{X\in \mathcal{W}_t} \inprod{-b}{z} +\inprod{X}{\Amap^*z-C}+\frac{\rho_t}{2}\norm{z-y_t}^2. 
\end{equation}
Using the strong duality of the above problem, we can interchange the min and max, and achieve 
\begin{equation}\label{eq: subproblemSpectralbundleMethodsolveStep1}
\begin{aligned} 
 &\min_{z} \max_{X\in \mathcal{W}_t} \inprod{-b}{z} +\inprod{X}{\Amap^*z-C}+\frac{\rho_t}{2}\norm{z-y_t}^2\\
= &\max_{X\in \mathcal{W}_t} \min_{z}  \inprod{-b}{z} +\inprod{X}{\Amap^*z-C}+\frac{\rho_t}{2}\norm{z-y_t}^2.
\end{aligned} 
\end{equation}
By completing the square, we find that the inner minimization is achieved only when $z =y_t + \frac{1}{\rho_t}\left(b-\Amap X\right),$ and 
the minimax problem reduces to 
\begin{equation}\label{eq: subproblemSpectralbundleMethodsolveStep2}
\begin{aligned} 
 &\max_{X\in \mathcal{W}_t} \min_{z}  \inprod{-b}{z} +\inprod{X}{\Amap^*z-C}+\frac{\rho_t}{2}\norm{z-y_t}^2\\ 
=&\max_{X\in \mathcal{W}_t} \inprod{-b}{y_t} +\inprod{X}{\Amap^*y_t-C} -\frac{1}{2\rho_t}\twonorm{b-\Amap X}^2. \\
= &-\min_{X\in \mathcal{W}_t} \inprod{b}{y_t} +\inprod{X}{C-\Amap^*y_t} +\frac{1}{2\rho_t}\twonorm{b-\Amap X}^2. \\ 
\end{aligned} 
\end{equation}
Interestingly, the last minimization problem in \eqref{eq: subproblemSpectralbundleMethodsolveStep2} 
is the augmented Lagrangian problem of Problem \eqref{p} with 
the decision variable $X$ restricted to $\mathcal{W}_{t+1}$ instead of $\symMat_+^\dm$.
So in essence, spectral bundle methods are secretly solving an augmented Lagrangian during its iterations. 

To solve the augmented Lagrangian problem in \eqref{eq: subproblemSpectralbundleMethodsolveStep2}, we may further 
rewrite it as 
\begin{equation}\label{eq: subproblemSpectralbundleMethodsolveStep3}
\begin{aligned} 
\min_{(\eta,S)\in \mathcal{S}_t}  
f_t(\eta,S),
\end{aligned} 
\end{equation}
where $f_t(\eta,S)$ and $\mathcal{S}_t$ are defined as 
\begin{equation}\label{eq: subproblemSpectralbundleMethodsolveStep4}
\begin{aligned} 
f_t(\eta,S) &:\;=\inprod{b}{y_t} +\inprod{\eta \bar{X}_t+V_tSV_t^\top}{C-\Amap^*y_t} +\frac{1}{2\rho_t}\twonorm{b-\Amap
	\left( \eta \bar{X}_t+V_tSV_t^\top \right)}^2, \\ 
\mathcal{S}_t & := \{{S\succeq 0,\;\eta \geq 0,\;\tr(S)+\alpha\eta \leq \alpha}\}.
\end{aligned} 
\end{equation}

\paragraph{Solving \eqref{eq: subproblemSpectralbundleMethodsolveStep3}} Since we can compute the gradient of $f_t$ easily and the projection to the constraint set $\mathcal{S}_t$ (after 
a proper scaling)
can be done with time complexity $\bigO(\bar{r}^3)$ (see Section \ref{sec: projecttoS_t} in 
Appendix for detail).  Hence, we may use the 
accelerated projected gradient method ~\cite{nesterov2013introductory}. Alternatively, we may use the interior 
point method described in \cite[Section 6]{helmberg2000spectral}, which has at least $\bigO({\bar{r}}^6)$ time complexity 
due to inverting an ${\bar{r}}^2 \times {\bar{r}}^2$ matrix. Moreover, in the special case $\bar{r}=1$, we have only two variables 
$(S,\eta)\in \RR^{2}$ with constraints $S\geq 0,\eta\geq 0$ and $S+\eta \leq \alpha$, and a quadratic objective $f_t$. Explicit formulas for the optimal $S_t^\star$ and $\eta_t^\star$ 
can be derived easily to avoid numerical optimization in this situation.

\paragraph{Storage concerns} We note that just for the purpose of computing $S^\star_t$ and $\eta^\star_t$, one needs not to store $\bar{X}_t$ but only need to store $\Amap(\bar{X}_t)$ and 
$\inprod{C}{\bar{X}_t}$ as we may write $f_t$ as 
\[
f_t(\eta,S) :\;=\inprod{b}{y_t} +\eta \inprod{ \bar{X}_t}{C} -\inprod{\Amap(V_tSV_t^\top)}{y_t} +\frac{1}{2\rho_t}\twonorm{b-\eta(\Amap
	\bar{X}_t)- \Amap(V_tSV_t^\top)}^2.
\]
The updates of $\Amap(\bar{X}_t)$ and 
$\inprod{C}{\bar{X}_t}$ are also easy given the low rank updates of $\bar{X}_t$ in \eqref{eqn: updatebarX_t+1FC} and \eqref{eq: aggregationUpdateOurMethod}. Keeping only 
$\Amap(\bar{X}_t)$ and 
$\inprod{C}{\bar{X}_t}$ is advantageous when $\Amap$ and $\inprod{C}{\cdot}$ can be quickly applied to low rank matrices. Moreover,
one can recover the matrix $\bar{X}_t$ without the need of storing $\bar{X}_t$ for both Block-Spec and HR-Spec 
using the matrix sketching idea in \cite{tropp2017practical}. 
We explain the sketching procedures more carefully in Section \ref{sec: matrixSketching} in the appendix. 
\paragraph{Summary of solving \eqref{eq: subproblemAMSpectralbundleMethod}} We summarize how to solve the minimax problem and obtain $z_t^\star$, $\eta_t^\star$, and $S^\star_t$ in the following two steps: 
\begin{enumerate}
	\item Solve \eqref{eq: subproblemSpectralbundleMethodsolveStep3} for general $\bar{r}>1$ via the accelerated gradient method or the interior point method, and obtain its optimal $S_t^\star$ and $\eta_t^\star$. The problem can be solved quickly for small $\bar{r}$ and admits explicit solution formulas for $\bar{r}=1$. 
	\item Compute $z_t^\star = y_t +\frac{1}{\rho_t}\left(b-\Amap(\eta_t^\star \bar{X}_t +V_tS_t^\star V_t^\top)\right)$.
\end{enumerate}






\section{Discussion}\label{sec: discussion}
In this paper, we give sublinear convergence rates of the classical spectral bundle method proposed by Helmberg and Rendl ~\cite{helmberg2000spectral}. We also 
develop Block-spec which not only enjoys the same rate, but also speed up to linear convergence with proper parameter choice and 
structural assumptions. We verify our theoretical results via numerical experiments. 

We here display a few potential future and (hopefully) interesting directions:
\begin{itemize}
	\item \textbf{Handling other constraints:} The problem format \eqref{p} only has equality constraints. Incorporating inequality 
	constraints and certain norm constraints such as $\norm{\Amap X-b}\leq \varepsilon$ for some $\varepsilon>0$ might be beneficial 
	for other applications of SDP such as 
	stochastic block model with more than 2 blocks \cite{amini2018semidefinite} and noisy matrix completion \cite{candes2010matrix}. 
	It seems 
	straightforward to extend the current framework to these new settings by introducing additional dual variables or analyzing 
	new dual objective.
	\item \textbf{Converging to the relative interior of the dual solution set $\ysolset$:} In Theorem \ref{thm: linear convergence of Block SBM under the extra condition strict complementarity},
	the rank estimate $\bar{r}$ needs to satisfy $\bar{r}\geq r_d$ instead of $\bar{r}\geq \rank(\xsol)$ assuming uniqueness of the primal solution. Though
	the quantity $r_d$ can be indeed larger than $\rank(\xsol)$ as shown in \cite[Theorem 5.1]{ding2020regularity}, $\bar{r}\geq \rank(\xsol)$ already ensures quick 
	convergence in our numerics. By speculating the proof of 
	Theorem \ref{thm: linear convergence of Block SBM under the extra condition strict complementarity}, 
	the linear convergence can be proved assuming Block-Spec converges to a dual solution that is in the relative interior of $\ysolset$. And 
	this is indeed what we observed after examining the dual slack matrices (not shown here). Of course, this 
	cannot be guaranteed by current algorithm design. Hence we ask the question that whether 
	 it is possible to design an algorithm that has low per iteration complexity and always converges to the relative interior of the optimal solution set.
	 \item \textbf{Adaptive choice of $\rho$ and $\bar{r}$:} Currently, the choice of $\rho$ and $\bar{r}$ are set to be a constant. Is it possible to have an algorithmic 
	 adaptive way in setting $\rho$ and $\bar{r}$? Does such choice hurt or enhance primal convergence?
\end{itemize}
\section{Introduction}
We consider solving semidefinite problems of the following form:
\begin{equation}\label{p} \tag{P}
\begin{aligned}
& \underset{X\in\symMat^{\dm}\subset \RR^{\dm \times \dm}}{\text{maximize}}
& & \langle -C, X \rangle \\
& \text{subject to}
& & \Amap X = b \\
&&& X \succeq 0,
\end{aligned}
\end{equation}
whose the decision variable $X$ is required to be symmetric and  positive definite: $X\succeq 0$. The  problem data comprises a symmetric cost matrix $C\in \symMat^{\dm}\subset \RR^{\dm \times \dm}$, a linear map $\Amap: \symMat^{\dm} \rightarrow \RR^{\ncons}$, and a right hand side vector $b \in \RR^{\dm}$. 
The task of solving \eqref{p} is often approached by considering its dual problem:
\begin{equation}\label{d} \tag{D}
\begin{aligned} 
& \underset{y\in\RR^{\ncons }}{\text{minimize}}
& & \langle -b, y \rangle \\
& \text{subject to}
& & \Amap^*y \preceq C.
\end{aligned}
\end{equation}

Semidefinite programming occurs at the heart of many important large-scale problems (for example, 
matrix completion \cite{candes2009exact}, max-cut \cite{goemans1995improved}, community detection \cite{bandeira2018random}, and phase 
retrieval \cite{candes2013phaselift}). A huge body of works has devoted to solve the SDP problem \eqref{p} \cite{todd2001semidefinite,nesterov1989self,nesterov1994interior,alizadeh1995interior,burer2003nonlinear,glowinski1975approximation,helmberg2000spectral,boyd2011distributed,friedlander2016low,friedlander2016low,yurtsever2019conditional,renegar2014efficient,ding2019optimal}. We refer the reader 
to \cite{monteiro2003first}, \cite[Section 2]{ding2019optimal}, and \cite[Section 3 and 4]{majumdar2019survey} for surveys of methods for solving the SDP problem. 

Among the methods, the spectral bundle method proposed by Helmberg and Rendl~\cite{helmberg2000spectral} is of particular interest of this paper due to its merits of low per iteration complexity and faster convergence in practice. These two properties make it ideal for solving large-scale problems where a high iteration cost may make computing even a single iteration prohibitively slow.
Interestingly, instead of dealing with \eqref{p} and \eqref{d} directly, Helmberg and Rendl's method considers an equivalent penalization yet unconstrained form of \eqref{d}:
for all sufficiently large $\alpha$, e.g., larger than the trace of any optimal solution of \eqref{p} \cite[Lemma 6.1]{ding2019optimal} \footnote{The lemma in \cite{ding2019optimal} actually 
	requires the primal solution to be unique. However, a closer look at the proof reveals that such 
	condition can be replaced by $\alpha >\sup_{\xsol \in \xsolset}\tr(\xsol)$.}, \eqref{d} is equivalent to (in the sense of same optimal value and solution set)\footnote{Actually, the HR 
	spectral bundle method requires the trace of every feasible $X$ for \eqref{p} to be the same and 
	deals with the eigenvalue  instead of the maximum of the eigenvalue and zero in \eqref{eq: penaltySDP}.
	However, the method can be translated word to word to the general situation where different feasible $X$ admits different trace.}
\begin{equation}\tag{\texttt{pen}-D}
\begin{aligned}\label{eq: penaltySDP}
& \underset{y\in\RR^{\ncons}}{\text{minimize}}
& & F(y) :\,= \langle -b, y \rangle +\alpha \max\{\lambda_{\max}(\Ajmap y-C), 0\}.
\end{aligned}
\end{equation}

In Section~\ref{sec: HR spectral bundle method}, we formally define the Helmberg and Rendl's Spectral Bundle Method, which we denote by HR-Spec. The main idea behind HR-Spec is to approximate the nonnegative eigenvalue function, $\alpha \max\{\lambda_{\max}(\Ajmap y-C), 0\}$, by a maximum of lower bound indexed by a small SDP representable set. Then using this model of the objective iterative improvement can be achieved by running the proximal bundle method~\cite{Du2017}.
Importantly, using a small SDP representable set makes the subproblems involved with the approximation of $\alpha \max\{\lambda_{\max}(\Ajmap y-C), 0\}$ and the bundle method easier to solve.

We defer the exact details constructing this model to Section~\ref{sec: Subgradient lower bound and spectral lower bound} and of the algorithm to Section~\ref{sec: HR spectral bundle method}. The per iteration complexity of such methods can be much lower than that of many ADMM type methods \cite{boyd2011distributed}, which requires projection to the SDP cone, $\symMat_+^\dm$, costing $\bigO(\dm^3)$ in general. 
Indeed, in Section~\ref{sec: ImportantsubproblemSolver}, we discuss the computation advantages of this approach and show that with proper choice of the parameters, only one eigenvector computation is required per iteration. 


HR-Spec has received considerable attention since it was first proposed (over 500 citations) and has developed various variants later on \cite{helmberg2002spectral,apkarian2008trust,helmberg2014spectral}. However, despite the success of this method, its convergence theory has left open two important issues. 
First, there is no convergence rate guarantees of HR-Spec. The main guarantee, Lemma 5 in \cite{helmberg2000spectral}, only ensures the dual iterate of \eqref{p} converges in terms of its objective value.
Second, there is no convergence guarantee on the primal side. In \cite{helmberg2000spectral}, only a paragraph in \cite[Section 8]{helmberg2000spectral} discusses how one may reconstruct the primal iterates without any guarantees and numerical verification on the quality of such primal reconstruction.  



\paragraph{Our contribution} We close these open areas in the theory for the spectral bundle method and propose a new variant of the method which enjoys greatly improved convergence rates: 
\begin{itemize}
	\item \textbf{Convergence rates for HR spectral bundle method}: In Theorem 
	\ref{thm: sublinearates}, we show that the HR-Spec admits 
	convergence rate $\bigO(\frac{\log(1/\epsilon)}{\epsilon^3})$ in terms of dual objective, and $\bigO(\frac{\log(1/\epsilon)}{\epsilon^6})$ in terms of primal convergence merely under strong duality. 
	With the extra condition strict 
	complementarity described in Section \ref{sec: analtical conditon}, dual objective and primal convergences 
	speed up to $\bigO(\frac{\log(1/\epsilon)}{\epsilon})$ and $\bigO(\frac{\log(1/\epsilon)}{\epsilon^2})$ respectively. 

	\item \textbf{Development of Block spectral bundle method (Block-Spec)}: We develop a new variant of HR-Spec, called 
	Block spectral bundle method. It not only enjoys 
	the same per iteration complexity and the theoretical guarantees of HR-Spec, but also enjoys \emph{linear convergence}
	whenever strict 
	complementarity holds, see Theorem \ref{thm: linear convergence of Block SBM under the extra condition strict complementarity}.
\end{itemize} 
We finally remark that both HR-Spec and Block-Spec are compatible with the matrix sketching idea in \cite{tropp2017practical} and hence 
can avoid storing a matrix $X$ with $\dm^2$ entries during its iteration and hence achieve the so-called storage optimality discussed in \cite[Section 1.2]{ding2019optimal}. An earlier attempt 
by the authors \cite{ding2019bundle} has explored this idea. Since storage is not a central topic of the paper, we defer the sketching procedure to Section \ref{sec: matrixSketching} in the appendix. 


\paragraph{Paper organization} The rest of the paper is organized as follows. In Section \ref{sec:def}, we describe 
the main ingredients of spectral bundle methods (spectral lower bounds, regularization, and aggregation) and formally define the Block-Spec and HR-Spec. In Section \ref{sec: analysis}, we present the convergence rates for both HR-Spec and Block-Spec. In Section \ref{sec: numerics}, we demonstrate that both methods are effective in solving modern problems admitting low rank solutions and valid our theory, finding when the conditions of our linear convergence theory are met, Block-Spec substantially speeds up.

\paragraph{Notation} We denote the sets of optimal solutions for \eqref{p} and \eqref{d} by $\xsolset$ and $\ysolset$ respectively. We equip 
$\symMat^{\dm}$ and $\real^{\ncons}$ with the trace inner product and dot product respectively, and denoted both as $\inprod{\cdot}{\cdot}$. The induced norms are both
denoted as $\norm{\cdot}$. For a symmetric matrix $A\in \symMat^\dm$, we write its eigenvalues as 
$\lambda_{\max}(A)=\lambda_{1}(A)\geq \dots \geq \lambda_{\dm}(A)$. The matrix 
operator norm (maximum of singular value), 
Frobenius norm, and nuclear norm (sum of singular values) are denoted as $\fronorm{\cdot}$, $\opnorm{\cdot}$, and $\nucnorm{\cdot}$ 
respectively. The dual slack matrix $Z(y)$ for each $y\in \RR^{\ncons}$ is defined as $Z(y):\, = C-\Amap^*y$. 
The operator norm of $\Amap^*$ is defined as $\opnorm{\Amap^*}=\max_{y\in \RR^{\ncons},\twonorm{y}\leq 1}\fronorm{\Amap^*y}$. For a closed set $\mathcal{X}\subset \RR^{\ncons}$ and a point $z\in \RR^{\ncons}$, we define the distance of $z$ to  $\mathcal{X}\subset \RR^{\ncons}$ as $\dist(z,\mathcal{X})=\inf_{x\in \mathcal{X}}\norm{x-z}$.


\section{Numerics}\label{sec: numerics}



In this section, we give preliminary numerics demonstrating that (i) both methods are able to 
solve modern problems  of interesting size, 
matrix completion and max-cut, and (ii) \emph{more importantly}, the rates given 
in Section \ref{sec: analysis} are not artifact of the proof but are what
observed in the experiments.

\paragraph{Experiments setup} The two SDP for matrix completion and max-cut are displayed in Table \ref{tb: SpecBundleMethodperformance}. 
For max-cut, $L$ is the Laplacian of the graph G1 in \cite{Gset} with $800$ vertices. 
For matrix completion,
$\Omega$ is the set of indices of the observed entries of the underlying rank $3$ matrix
$\trux \in \RR^{400\times 400}$. Here $\trux = WW^\top$ where $W\in \RR^{400\times 3}$ with each 
entry following Rademacher distribution. Each entry of $\trux$ is observed with 
probability $0.1$. Both problem has decision variable size $800\times 800$, and 
satisfy strong duality, strict complementarity, and has a unique primal solution 
$\xsol$ as verified in \cite{ding2020regularity}. 
We set $\alpha =2\dm$, $\rho=0.5,\beta = 0.25$ for Max-Cut and $\alpha = 4 \nucnorm{\trux},\rho =5,\beta= 0.25$ for matrix completion. 
Both methods start the initial $X_0,y_0,z_0$ at zero. Subproblem \eqref{eq: subproblemAMSpectralbundleMethod} is solved via Mosek \cite{mosek2010mosek}.
We run 100 iterations for matrix completion, and 200 iterations for max-cut.
The optimal value $\pval$ and primal solution $\xsol$ for max-cut is obtained through Mosek \cite{mosek2010mosek}. 
We set 
$\pval = 2\nucnorm{\trux}$ and $\xsol = \begin{bmatrix}
\trux & \trux \\
\trux & \trux 
\end{bmatrix}$ for matrix completion. Such choice of $\xsol$ indeed solves matrix completion SDP 
with high probability \cite{candes2009exact}. 
Let the rank of the optimal solution be $\rsol = \rank(\xsol)$, which is $13$ for max-cut and 
$3$ for matrix completion. We set the rank estimate $\bar{r}= \rsol-1$, $\rsol$, and $\rsol+1$.

\paragraph{Experiment Results} The experiment results are shown in Table \ref{tb: SpecBundleMethodperformance}, the accuracy of the last iterates, and Figure 
\ref{fig: penalizedObjectiveValue}, the evolution of dual objective value $F$ in \eqref{eq: penaltySDP}. 
The dual optimality (dual opt.), primal optimality (primal opt.), and primal feasibility (primal feas.) 
are defined as $\frac{F(y)-\dval}{\dval}$, $\abs{\frac{\inprod{C}{X}-\pval}{\pval}}$, and $\frac{\norm{\Amap X-b}}{\norm{b}}$ respectively. 

\paragraph{Solving problem to moderate accuracy} As shown in Table \ref{tb: SpecBundleMethodperformance},  whenever $\bar{r}\geq \rsol$,
HR-Spec is able to solve both problems in terms of dual optimality 
 to moderate low accuracy ($10^{-3}\sim 10^{-4}$), while Block-Spec is able 
 to solve both problems in terms of dual optimality to moderate high accuracy ($10^{-7}\sim 10^{-8}$). We note that 
 Block-spec actually achieves $10^{-4}\sim 10^{-5}$ accuracy within 50 iterations for both problems 
 as shown in Figure \ref{fig: penalizedObjectiveValue}. We include additional 
 50 iterations of matrix completion and 150 iterations of max-cut to see how HR-spec
 performs. Since most of the time, it only achieves $10^{-2}$ accuracy after 50 iterations. 
 In terms of primal iterate $X_t$, primal feasibility is usually worse than dual optimality by one or two order of magnitude, while
 primal optimality is usually on the same order. The fact primal feasibility is a bit worse might 
 be expected from our lemmas in Section \ref{sec: Relating primal and dual convergence}, as we indicate that bound 
 for primal feasibility can be worse than dual optimality. 
 
 \paragraph{Sublinear and linear rates} As shown in Figure \ref{fig:figure Matrix Completion}, HR-Spec may exhibit sublinear 
 convergence rate for the case $\bar{r}=2,3$. In Figure \ref{fig:figure Max-cut}, the method also converges 
 slowly indicating a potential sublinear rates. Block-spec does converge very fast in both case whenever $\bar{r}\geq \rsol$ 
 as expected from our Theorem \ref{thm: linear convergence of Block SBM under the extra condition strict complementarity}. It flats 
 out after reaching $10^{-8}$ accuracy, which we suspect is due to the inaccuracy in the eigenvalue solver or the choice of $\rho$. Block-spec 
 also has slower convergence whenever $\bar{r}<\rsol$ indicating the condition $\bar{r}\geq \rsol$ is crucial for quick convergence.
 
\begin{table}
\begin{tabular}{clllll}
	\hline 
   Problem & $\bar{r}$ & Method & Dual Opt. & Primal Opt. & Primal Feas.\\ 
	\hline 
\textbf{Matrix Completion}                                 	 & \multirow{2}{*}{2} & HR & $0.6177$   & 0.2621 &  0.3707 \\ 
\multirow{5}{*}{$
	\displaystyle
	\begin{aligned}
	&{\text{max}}
	& & -\inprod{I}{W_1}-\inprod{I}{W_2} \\
	& \text{s.t.}
	& &  X_{ij} = \trux_{ij},\, (i,j)\in \Omega \\
	&  & &
	\begin{bmatrix}
	W_1 & X \\ X^\top & W_2
	\end{bmatrix}\succeq 0.
	\end{aligned}
	$}                                                        &   &  B & 0.07236 & 0.01155 &  $0.01308$ \\ 
                                                                    & \multirow{2}{*}{3} & HR & 0.1898 & $1.2146 \times 10^{-11}$ & 0.1448  \\ 
                                                                	&  & B & $2.788 \times  10^{-7}$& $8.441\times 10^{-7}$ & $	1.098\times 10^{-5}$ \\ 
                                                                    & \multirow{2}{*}{4}  & HR &  $7.0455\times 10^{-7}$ & $4.782\times 10^{-11}$  & $1.547\times 10^{-4}$  \\ 
	                                                                &  &  B  & $4.967\times 10^{-8}$ & $2.823\times 10^{-8}$ &  $9.154\times 10^{-6}$\\ 
	\hline \hline 

\textbf{Max-cut}                                 	& \multirow{2}{*}{12} &  HR & $8.740 \times 10^{-4}$ & $6.092\times 10^{-4}$ & 0.04898  \\ 
	\multirow{5}{*}{$
\displaystyle
\begin{aligned}
&\mbox{max} & & \inprod{L}{X}\\
&\mbox{s.t.} & & \diag(X) =  \ones \\
& & & X\succeq 0
\end{aligned}
$
}
&  &  B &  $6.670\times 10^{-6}$ &  $3.550\times 10^{-6}$& $6.092 \times 10^{-4}$ \\ 
& \multirow{2}{*}{13} &  HR & $5.357\times 10^{-4}$ & $2.185\times 10^{-4}$ & 0.03988 \\ 
&  & B  & $6.234\times 10^{-10}$ & $6.988\times 10^{-9}$& $9.772\times 10^{-6}$  \\ 
& \multirow{2}{*}{14} & HR & $1.278\times 10^{-4}$ & $1.574\times 10^{-4}$ & 0.01574\\ 
&  & B  &  $1.207\times 10^{-8}$ & $2.281\times 10^{-7}$  & $2.733\times 10^{-5}$\\ 
	\hline 
\end{tabular} 
\caption{This table shows the final accuracy of $y_t$ and $X_t$ for block spectral bundle method (B) and 
	HR spectral bundle method (HR) after $100$ iteration for Matrix completion, and $200$ iterations for 
Max-cut problem with varying rank estimate $\bar{r}$. The dual optimality (dual opt.), primal optimality (primal opt.), and primal feasibility (primal feas.) 
are defined as 
 $\frac{F(y)-\dval}{\dval}$, $\abs{\frac{\inprod{C}{X}-\pval}{\pval}}$, and $\frac{\norm{\Amap X-b}}{\norm{b}}$ respectively. The set 
$\Omega$ is the set of indices for observed entries of the ground truth $\trux$ for matrix completion. The matrix $L$ is the 
Laplacian matrix of the graph.}\label{tb: SpecBundleMethodperformance}
\end{table} 

\begin{figure}[H]
	\begin{subfigure}[(a)]{.5\textwidth}
		\centering
		\includegraphics[width=0.8\linewidth, height= 0.25\textheight]{MatrixCompletionPlot.pdf}
		\caption{Matrix Completion}
		\label{fig:figure Matrix Completion}
	\end{subfigure}%
	\hspace{5pt}
	\begin{subfigure}[(b)]{.5\textwidth}
		\centering
		\includegraphics[width=0.8\linewidth, height= 0.25\textheight]{Max-cutPlot.pdf}
		\caption{Max-cut}
		\label{fig:figure Max-cut}
	\end{subfigure}
	\caption{The evolution of the relative penalized dual objective value $\abs{\frac{F(y_t)-F(\ysol)}{{F(\ysol)}}}$. Here B is Block-Spec shown 
		as the solid line, and HR is HR-spec, shown as the dotted line.}
	\label{fig: penalizedObjectiveValue}
\end{figure}
\section{Lemmas, proofs, and procedures for Section \ref{sec:def}}\label{sec: detailsForSectionref{sec:def}}
\section{Inequalities \eqref{eq: subgradientAndAggregation} of HR bundle method}\label{sec: Model comparison for HR bundle method}
We aim to show that the property \eqref{eq: subgradientAndAggregation} continues to hold for HR bundle method. 
This is verified if we establish the following membership relation:
\begin{equation*} 
X_t\in \mathcal{W}_{t+1}\quad\text{and}\quad \alpha  v_1v_1^\top \in \mathcal{W}_{t+1}.
\end{equation*}
We consider the general situation that we may choose $r$ many eigenvectors of $S_t^\star$ instead 
of $\bar{r}-1$ as stated in step 3' in Section \ref{sec: Block spectral bundle method}. 

The membership of $X_t$ in $\mathcal{W}_{t+1}$ is by construction. To see this,
recall the number of column of $V_tQ_1$ is $r$ as we pick $r$ eigenvectors of 
$S_t^\star$. Denote the left $\dm\times r$ matrix of $R$ by $R_{:,1:r}$. 
We have $V_tQ_1  = Q R_{:,1:r}$.
Then we can set $\eta = 1$, and the left top $r\times r$ submatrix of $S$ as $S_{1:r,1:r}=R_{:,1:r}\Lambda_1R_{:,1:r}^\top$, and the rest of $S$ to be $0$. Then we have 
$V_{t+1}SV_{t+1}^\top = V_tQ_1\Lambda_1 Q_1^\top V_t^\top$. This choice of $\eta$ and $S$ 
makes $X_t = \eta \bar{X}_{t+1} + V_{t+1}SV_{t+1}^\top $ due to the updating scheme \eqref{eqn: updatebarX_t+1FC} of $\bar{X}_{t+1}$ 
and definition of $X_t$ in \eqref{eq: X_tdefinition}. This choice of $S$ is feasible because $S\succeq 0$ as $\Lambda_1\succeq 0$, and 
\begin{equation}
\begin{aligned}  \label{eq: X*intheNewSettrace}
\eta \tr (\bar{X}_{t+1})+\tr(S) & \overset{(a)}{=} \tr(\eta^\star _t\bar{X}_t)+ \tr(V_tQ_2\Lambda_2Q_2^\top V_t^\top) + \tr(V_tQ_1\Lambda_1 Q_1^\top V_t^\top)\\
&\leq \eta_t^\star\alpha +\tr(V_tS^\star_tV_t^\top) \overset{(b)}{\leq} \alpha.
\end{aligned} 
\end{equation} 
Here in step $(a)$ we use the definiton of $\bar{X}_{t+1}$ and $\tr(S)= \tr(V_{t+1}SV_{t+1})$ 
because $V_{t+1}$ has orthonormal columns. Step $(b)$ is due to $V_t$ has orthonormal columns 
and the iterate $\eta^\star _t$ and $S^\star _t$ satisfies the constraint $\eta^\star \alpha  +\tr (S^\star _t)\leq \alpha$ by construction.

\section{Projecting to a scaled $\mathcal{S}_t$}\label{sec: projecttoS_t}
Recall we are trying to solve 
\[
\min_{(\eta,S)\in\mathcal{S}_t} f_t(\eta,S)
\]
\begin{equation} 
\begin{aligned} 
	f_t(\eta,S) &:\;=\inprod{b}{y_t} +\inprod{\eta \bar{X}_t+V_tSV_t^\top}{C-\Amap^*y_t} +\frac{1}{2\rho_t}\twonorm{b-\Amap
		\left( \eta \bar{X}_t+V_tSV_t^\top \right)}^2. \\ 
	\mathcal{S}_t & := \{{S\succeq 0,\;\eta \geq 0,\;\tr(S)+\alpha \eta \leq \alpha}\}
\end{aligned} 
\end{equation}
After a proper scaling of $S$, we may consider
the constraint set as 
\[ 
\tilde{\mathcal{S}}= \{S\in \symMat_+^k,\, \eta\geq 0,\, \tr(S)+\eta \leq 1\},
\]
and a new objective $\tilde{f}_t(\eta,S) = f(\eta,\alpha S)$.

Let us explain how to project to the set $\tilde{\mathcal{S}}$ for arbitrary 
$(\eta_0,S_0)\in \RR\times \symMat^{\bar{r}}$ and get its projection $(\eta^\star,S^\star)$. The procedures are as follows:
\begin{enumerate}
	\item Compute the eigenvalue decomposition of $S_0 = V\Lambda_0 V^\top$,
	where $\Lambda_0\in \symMat^{\bar{r}}$ is a diagonal matrix with diagonal $\vec{\lambda}_0=( \lambda_1,\dots,\lambda_{\bar{r}})$.
	\item Compute $(\eta^\star, \vec{\lambda}^\star) = \arg\min_{\eta+\sum_{i=1}^{\bar{r}}\lambda_i \leq 1,\;\eta \geq 0,\;\lambda_i\geq 0}\twonorm{(\eta_0,\vec{\lambda}_0)-(\eta,\vec{\lambda})}$.
	\item Form $S^\star = V\diag(\vec{\lambda^\star})V^\top$.
	Here $\diag(\lambda)$ forms a diagonal matrix with the vector $\lambda$ on the diagonal.

\end{enumerate}
The main computational cost is the eigenvalue decomposition which requires $\bigO(\bar{r}^3)$ time. The second 
step requires projection to the convex hull of probability simplex and the origin, which can be done 
in $\bigO(k\log k)$ time \cite{wang2013projection}. The correctness of the procedure
can be verified as in \cite[Lemma 3.1]{allen2017linear}
and \cite[Lemma 6]{garber2019convergence}.

\section{Matrix Sketching} \label{sec: matrixSketching}
When $m$ is on the order $\dm$, it is beneficial to use the matrix sketching 
idea developed in  \cite{tropp2017practical} to 
achieve storage reduction. As noted in \ref{sec: ImportantsubproblemSolver},
if we store $\Amap (X_t)= d_t$ and $c_t = \inprod{C}{X_t}$ at each iteration, then we have no problem in 
doing the small-scale SDP \eqref{eq: subproblemSpectralbundleMethodsolveStep3}. If
$\Amap$ and inner product with $C$ can be applied to low rank matrices efficiently, then 
updating $z_t$ and $c_t$ is not hard  
due to linearity of our updating scheme $\bar{X}_{t+1} = X_t = \eta^\star_t\bar{X}_t+V_t S^\star_t V_t^\top$ in 
\eqref{eq: aggregationUpdateOurMethod} for Block-spec, and  
the updating scheme of $\bar{X}_t$ in \eqref{eqn: updatebarX_t+1FC} for HR-spec.
 
Let us first explain how to  omit storing the iterate $\bar{X}_{t+1}=X_t = \eta^\star_t\bar{X}_t+V_t S^\star_t V_t^\top$
for Block-spec. 
We first draw two matrices with independent standard normal  entries
\begin{equation}
\begin{aligned}
\Psi \in \real^ {n \times k} \quad \text{with} \quad k=2r+1; \\
\Phi \in \real ^{l\times \dm} \quad \text{with} \quad l=4r+3.
\end{aligned}\nonumber
\end{equation} 
Here $r$ is chosen by the user. It either represents the estimate of the true rank of the primal solution or the user's computational budget in dealing with larges matrices.

We use $Y^C_t$ and $Y^R_t$ to capture the column space  and the row space  of $X_t$:
\begin{equation}\label{eqn: onetimeSketch}
Y^C_t =\bar{ X}_t\Psi \in \real^{\dm \times k},\qquad  Y^R_t =\Phi \bar{X}_t \in \real^{l\times n} .
\end{equation}
Hence we initially have $Y^C_0=0$ and $Y^R_0=0$.
Notice that we should not store matrix $X_t$ in this case. However, using the updating scheme
$X_{t+1} = \eta^\star_t\bar{X}_t+V_t S^\star_t V_t^\top$. 
where $V_t\in \RR^{\dm \times \bar{r}}$ and $S_t^\star \in \symMat^{\bar{r}}$. 

In this setting, $Y^C_{t+1}$ and $Y^R_{t+1}$ can be directly computed as
\begin{align}
Y^C_{t+1} = V_tS^\star _t(V_t^\top\Psi)+ \eta_t^\star  Y^C_t \in \real^{\dm \times k }, \label{eq:sketch-primal-update1}\\
Y^R_{t+1} = (\Psi V_t)S_t^\star V_t^\top+ \eta_t^\star  Y^R_t\in \real^{l\times n} . \label{eq:sketch-primal-update2}
\end{align}
This observation allows us to form the sketch $Y^C_t$ and $Y^R_t$ from the stream of updates. 

We then reconstruct $X_t$ and get the reconstructed matrix $\hat{X}_t$ by
\begin{align}\label{eqn: reconstructionsketch}
Y^C_t = Q_tR_t, \quad  B_t = (\Phi Q_t)^{\dagger} Y^R_t, \quad \hat{X}_t= Q_t[B_t]_r,
\end{align}
where $Q_tR_t$ is the $QR$ factorization of $Y^C_t$ and $[\cdot]_r$ returns the best rank $r$ approximation in Frobenius norm.
Specifically, the best rank $r$ approximation of a matrix $Z$ is $U\Sigma V^\top $,
where $U$ and $V$ are right and left singular vectors corresponding to the $r$
largest singular values of $Z$ and $\Sigma$ is a diagonal matrix with $r$ largest singular values of $Z$. In actual implementation, we 
may only produce the factors $(QU, \Sigma, V)$ defining $\hat{X}_T$ in the end instead of reconstructing $\hat{X}_t$ in every iteration. We refer the reader to \cite[Theorem 5.1]{tropp2017practical} for the theoretical guarantees on the reconstruction matrix $\hat{X}_t$. 

Hence we can avoid the \emph{forming a new iteratre} procedure in Block-spec. We remark that 
the reconstructed matrix $\hat{X}_t$ is not necessarily positive semidefinite. However, this suffices for the purpose of finding a matrices close to $X_t$. More sophisticated procedure is available for producing a positive semidefinite approximation of $X_t$ \cite[Section 7.3]{tropp2017practical}.  

For HR-spec, using the definition of $\bar{X}_t$ in \eqref{eqn: updatebarX_t+1FC}, we can sketch 
it via the same procedure as $X_t$ for Block-spec as above using two matrices $Y_{t}^C$ and 
$Y_t^R$. To construct the converging primal iterate $X_t$, we note from the 
definition of $X_t$ in  \eqref{eq: X_tdefinition}: $X_{t} = \eta^\star_{t+1}\bar{X}_{t}+V_{t+1} S^\star_{t+1} V_{t+1}^\top$,
we can sketch $X_t$ by constructing $\tilde{Y}_{t+1}^C$ and $\tilde{Y}_{t+1}^R$
 as the  $Y_{t+1}^C$ and $Y_{t+1}^R$ in \eqref{eq:sketch-primal-update1}
and \eqref{eq:sketch-primal-update2}, and reconstruct it using \eqref{eqn: reconstructionsketch} 
from $\tilde{Y}_{t+1}^C$ and $\tilde{Y}_{t+1}^R$.

\section{Proof of Lemma \ref{lem: importantLemmaQuadraticAccurateModel}}\label{sec: proofOfImportantLemma}
In this section, we provide a proof of Lemma \ref{lem: importantLemmaQuadraticAccurateModel}. In addition
to the Frobenius norm bound, we also provide an operator two norm bound \eqref{eq:qudraticaccuracyopnorm}. 

Recall the assumption that $\lambda_{r}(X)-\lambda_{r+1}(X)=\delta>0$ for
some $\delta>0$. Let $V\in\real^{n\times r}$ be an orthonormal matrix
formed by the $r$ eigenvectors corresponding to the top $r$ eigenvalues.
Hence if $V$ is an representation, then any $VO$ with $O\in\real^{r\times r}$
orthonormal is also a representation. Recall $r$-th spectral plus set $\faceplus r(X):=\{VSV^{\top}\mid\tr(S)\leq1,S\succeq0,S\in\symMat^{r}\}.$
Note that any $VO$ produce the same $\faceplus r(X)$. 

Now for any $Y\in\symMat^{\dm},$ since $\max\{\lambda_{1}(Y),0\}=\max_{W\succeq0,\tr(W)\leq1,}\inprod WY,$
we see the following always holds (as $\{W|W\succeq0,\tr(W)\leq1\}\supset\faceplus r(X)$):
\begin{align*}
\max\{\lambda_{1}(Y),0\} & \geq\max_{W\in\faceplus r(X)}\inprod WY.
\end{align*}
Define the error $f_X(Y)$ as 
\begin{align*}
f_{X}(Y) & =\lambda_{1}(X)-\max_{W\in\faceplus r(X)}\inprod WY.
\end{align*}
We always have $f_{X}(Y)\geq0$ as previously argued. If $\lambda_{1}(Y)<0$,
then $\max\{\lambda_{1}(Y),0\}=0$ and hence $Y\preccurlyeq0$. Thus
the approximation $\max_{W\in\faceplus r(X)}\inprod WY=0$ as well.
Hence we may only consider the case $\lambda_{1}(Y)>0$ in the following.
Let $v$ be the eigenvector with two norm $\twonorm v=1$ corresponding
to the largest eigenvalue $\lambda_{1}(Y)$, then 
\begin{align*}
f_{X}(Y) & =\lambda_{1}(Y)-\max_{W\in\faceplus r(X)}\inprod WY  =\min_{W\in\faceplus r(X)}\inprod{vv^{\top}-W}Y\\
         & =\min_{W\in\faceplus r(X)}\underbrace{\inprod{vv^{\top}-W}{Y-X}}_{T_{1}}+\underbrace{\inprod{vv^{\top}-W}X}_{T_{2}}.
\end{align*}
To analyze $T_{1}$ and $T_{2}$, we define some notaion first. Denote
$V'\in\mathbb{\mathbb{R}}^{\dm\times r}$ to be the orthonormal matrix
formed by the eigenvectors corresponding to the top $r$ eigenvalues
of $Y$. We assume the first column of $V'$ is $v.$ Also denote
$F\in\real^{n\times(n-r)}$ to be an orthonormal matrix fomed by the
rest eigenvectors of $X.$ So the eigenvalue decomposition of $X$
is $X=V\Lambda_{1}V^{\top}+F\Lambda_{2}F^{\top},$ for some diagonal
$\Lambda_{1}\in\symMat^{r}$ and $\Lambda_{2}\in\symMat^{(\dm-r)\times(n-r)}$.

Let us now analyze the $T_{2}$ term first. We may choose $W_{1}=VV^{\top}vv^{\top}VV^{\top}$
here. With such choice, $T_{2}$ equals the following:
\begin{align*}
T_{2} & =\inprod{vv^{\top}-W}X  =\inprod{vv^{\top}}X-\inprod WX\\
      & \overset{(a)}{=}\inprod{vv^{\top}}{V\Lambda_{1}V^{\top}+F\Lambda_{2}F^{\top}}-\inprod{VV^{\top}vv^{\top}VV^{\top}}{V\Lambda_{1}V^{\top}+F\Lambda_{2}F^{\top}}\\
      & \overset{(b)}{=}\inprod{vv^{\top}}{V\Lambda_{1}V^{\top}}+\inprod{vv^{\top}}{F\Lambda_{2}F^{\top}}-\inprod{vv^{\top}}{V\Lambda_{1}V^{\top}}\\
      & \overset{(c)}{=}\inprod{V'e_{1}(V'e_{1})^{\top}}{F\Lambda_{2}F^{\top}}.
\end{align*}
Here we use the eigenvalue decomposition of $X$ in step $(a)$. The
step $(b)$ uses the fact that $V$ has orthonormal columns and that
$V^{\top}F=0$ as the columns are orthonormal. Step $(c)$ uses the
fact that $v$ is the first column of $V'$.

Let the 
matrix $O^{\star}\in\real^{r\times r}:O^{\star}\in\arg\min_{OO^\top =I}\fronorm{VO-V'}$.
Below, we denote $VO^{\star}$ to be $V$ to save some notation. Let
the error of $V$ and $V'$ be $E=V'-V.$ Also denote $e_{1}\in\real^{r}$
to be the vector with first entry $1$ and all other entries $0$. 
Let us now bound $T_{2}=\inprod{V'e_{1}(e_{1}V')^{\top}}{F\Lambda_{2}F^{\top}}$:
\begin{equation} \label{eq:T2maxlambda}
\begin{aligned}
T_{2} 
& =\inprod{(V+E)e_{1}e_{1}^{\top}(V+E)^{\top}}{F\Lambda_{2}F^{\top}
& \overset{(a)}{=}&\inprod{Ee_{1}e_{1}^{\top}E^{\top}}{F\Lambda_{2}F^{\top}}
\\
&\overset{(b)}{\leq}\nucnorm{Ee_{1}e_{1}^{\top}E^{\top}}\opnorm{F\Lambda_{2}F^{\top}}
& \overset{(c)}{=}&\opnorm{Ee_{1}e_{1}^{\top}E^{\top}}\opnorm{F\Lambda_{2}F^{\top}
\\
& \overset{(d)}{\leq}\opnorm{Ee_{1}}^{2}\opnorm{\Lambda_{2}
& \overset{(e)}{\leq}&\opnorm E^{2}\opnorm{\Lambda}
\end{aligned}
\end{equation} 
 Here we use the fact $V^{\top}F=0$ in step $(a)$. Step $(b)$ is
due to the Holder's inequality. Step $(c)$ uses the fact that for
rank $1$ matrix, the Frobenius norm is the same as its operator norm.
Step $(d)$ uses the submultiplicity of operator two norm. The last
step $(e)$ uses the fact operator norm of $e_{1}$ is $1$. 

Let us now turn back to analyze $T_{1}$. With the choice of $W=VV^{\top}vv^{\top}VV^{\top}$,
the difference $vv^{\top}-W$ is 
\begin{align*}
vv^{\top}-W & =V'e_{1}(V'e_{1})^{\top}-VV^{\top}vv^{\top}VV^{\top}\\
 & =V'e_{1}(V'e_{1})^{\top}-VV^{\top}V'e_{1}(V'e_{1})^{\top}VV^{\top}\\
 & =(V+E)e_{1}e_{1}^{\top}(V+E)^{\top}-VV^{\top}(V+E)e_{1}e_{1}^{\top}(V+E)^{\top}VV^{\top}\\
 & =Ee_{1}e_{1}^{\top}V^{\top}+Ve_{1}e_{1}^{\top}E^{\top}+Ee_{1}e_{1}E^{\top}\\
 & -VV^{\top}Ee_{1}e_{1}^{\top}V^{\top}-Ve_{1}e_{1}^{\top}E^{\top}VV^{\top}-VV^{\top}Ee_{1}e_{1}^{\top}EVV^{\top}.
\end{align*}
Hence, using the fact the nuclear norm of rank one matrix is the same
as operator norm, the nuclear norm of $vv^{\top}-W$ is bounded by
\begin{align*}
\nucnorm{vv^{\top}-W} & \leq\opnorm{Ee_{1}e_{1}^{\top}V^{\top}}+\opnorm{Ve_{1}e_{1}^{\top}E^{\top}}+\opnorm{Ee_{1}e_{1}E^{\top}}\\
 & +\opnorm{VV^{\top}Ee_{1}e_{1}^{\top}V^{\top}}+\opnorm{Ve_{1}e_{1}^{\top}E^{\top}VV^{\top}}+\opnorm{VV^{\top}Ee_{1}e_{1}^{\top}EVV^{\top}}\\
 & \overset{(a)}{\leq}4\opnorm E+2\opnorm E^{2}.
\end{align*}
Here in step $(a)$, we use the fact that $\opnorm{e_{1}e_{1}^{\top}}\leq1$
and $\opnorm V\leq1.$ Thus the first term $T_{1}$ is bounded by
\begin{equation}\label{eq: T1maxlambda}
T_{1} 
=\inprod{vv^{\top}-W}{Y-X}
\overset{(a)}{\leq}\nucnorm{vv^{\top}-W}\opnorm{Y-X}
\leq\left(4\opnorm E+2\opnorm E^{2}\right)\opnorm{Y-X}
\end{equation}

Now combining (\ref{eq: T1maxlambda}) and (\ref{eq:T2maxlambda}),
we found that $f_{X}(Y)$ is upper bounded by 
\begin{align*}
f_{X}(Y) & \leq\opnorm{\Lambda_{2}}\opnorm E^{2}+\left(4\opnorm E+2\opnorm E^{2}\right)\opnorm{Y-X}.
\end{align*}

Let us consider two cases:
\begin{enumerate}
\item First is the case of Frobenius norm. The Frebenius bound \cite[Theorem 2]{yu2015useful} of $E$ says
that 
\begin{align*}
\fronorm E & \leq\frac{2\sqrt{2}\fronorm{Y-X}}{\delta}.
\end{align*}
Hence in this case, we have for all $Y\in\symMat^{\dm}$
\begin{align*}
f_{X}(Y) & \leq\frac{8\fronorm{Y-X}^{2}\opnorm{\Lambda_{2}}}{\delta^{2}}+\frac{8\sqrt{2}\fronorm{Y-X}^{2}}{\delta}+\frac{16\fronorm{Y-X}^{3}}{\delta^{2}}.
\end{align*}
\item Second is the operator norm. We have the operator norm of $E$ bounded
by 
\begin{align*}
\opnorm E & \leq\frac{2\sqrt{2}\sqrt{r}\opnorm{Y-X}}{\delta}.
\end{align*}
In this case, the function $f_{X}(Y)$ is upper bouned by 
\begin{align}
f_{X}(Y) & \leq\frac{8r\opnorm{Y-X}^{2}\opnorm{\Lambda_{2}}}{\delta^{2}}+\frac{8\sqrt{2}\sqrt{r}\opnorm{Y-X}^{2}}{\delta}+\frac{16r\opnorm{Y-X}^{3}}{\delta^{2}}.\label{eq:fxYopeartornormboundglobalwithr}
\end{align}
If $\opnorm{Y-X}\leq\delta$, then using it for the term $\frac{16r\opnorm{Y-X}^{3}}{\delta^{2}}$,
we have 
\begin{align}
f_{X}(Y) & \leq\frac{8r\opnorm{Y-X}^{2}\opnorm{\Lambda_{2}}}{\delta^{2}}+\frac{(8\sqrt{r}+16r)\opnorm{Y-X}^{2}}{\delta}.\label{eq:fxYoperatornormboundwithr}
\end{align}
\end{enumerate}
Still, we have not reach a global quadratic accurate model. Let us show that the function $f_{X}(Y)$ is bounded by a linear
difference always. We decompose $f_{X}(Y)$ into two terms: 
\begin{align*}
f_{X}(Y) & =\max\{\lambda_{1}(Y),0\}-\max\{\lambda_{1}(V^{\top}YV),0\}\\
 & =\underbrace{\max\{\lambda_{1}(Y),0\}-\max\{\lambda_{1}(X),0\}}_{T_{1}}+\underbrace{\max\{\lambda_{1}(X),0\}-\max\{\lambda_{1}(V^{\top}YV),0\}}_{T_{2}}.
\end{align*}
 For the term $T_{1}$, we note the function $\max\{x,0\}$ for any
$x\in\real$ is $1$-Lipschitz with respect to the norm
$|x|.$ Thus the term $T_{1}$ is bounded by 
\begin{equation*}
|T_{1}|
\leq|\lambda_{1}(Y)-\lambda_{1}(X)|

 \leq\opnorm{Y-X}.
\end{equation*}

For the second term $T_{2}$, we note that $\lambda_{1}(X)=\lambda_{1}(V^{\top}XV$)
because of the definition of $V$. Hence, use the same reasoning,
we have 
\begin{equation*}
|T_{2}| 
\leq\opnorm{V^{\top}XV-V^{\top}YV}

 \leq\opnorm{X-Y}.
\end{equation*}
 where the last line is due to submultiplicity of operator two norm.
Hence, we see the error function $f_{X}(Y)$ is always bounded by
\begin{align*}
|f_{X}(Y)| & \leq2\opnorm{X-Y}.
\end{align*}

The inequality (\ref{eq:fxYoperatornormboundwithr}) tells us that
when $\opnorm{X-Y}\leq\delta$, we have $f_{X}(Y)\leq\frac{8r\opnorm{Y-X}^{2}\opnorm{\Lambda_{2}}}{\delta^{2}}+\frac{(8\sqrt{2r}+16r)\opnorm{Y-X}^{2}}{\delta}.$
Now if $\opnorm{X-Y}\geq\delta$, then definitely $\frac{2\opnorm{X-Y}^{2}}{\delta}\geq2\opnorm{X-Y}$.
Hence, the model $\max_{W\in\faceplus r(X)}\inprod WY$ is always
qudratic accurate: for all $Y\in\symMat^{\dm},$
\begin{align}
\text{0\ensuremath{\leq f_{X}\left(Y\right)\leq}} & \min\left\{ \frac{8r\opnorm{Y-X}^{2}\opnorm{\Lambda_{2}}}{\delta^{2}}+\frac{(8\sqrt{2r}+16r)\opnorm{Y-X}^{2}}{\delta},\frac{2\opnorm{X-Y}^{2}}{\delta}\right\} .\label{eq:qudraticaccuracyopnorm}
\end{align}
The same argument applies to the Frobeius norm case, and we reach 
\begin{align}\label{eq: qudraticaccuracyFrobeniusNorm}
f_{X}(Y) & \leq\frac{8\fronorm{Y-X}^{2}\opnorm{\Lambda_{2}}}{\delta^{2}}+\frac{(8\sqrt{2}+16)\fronorm{Y-X}^{2}}{\delta}.
\end{align}
2020-08-18T02:22:55.000
2008
2008.07067
en
https://arxiv.org/abs/2008.07067
The spectral bundle method proposed by Helmberg and Rendl is well established for solving large-scale semidefinite programs (SDP) thanks to its low per iteration computational complexity and strong practical performance. In this paper, we revisit this classic method showing it achieves sublinear convergence rates in terms of both primal and dual SDPs under merely strong duality, complementing previous guarantees on primal-dual convergence. Moreover, we show the method speeds up to linear convergence if (1) structurally, the SDP admits strict complementarity, and (2) algorithmically, the bundle method captures the rank of the optimal solutions. Such complementary and low rank structure is prevalent in many modern and classical applications. The linear convergent result is established via an eigenvalue approximation lemma which might be of independent interest. Numerically, we confirm our theoretical findings that the spectral bundle method, for modern and classical applications, speeds up under these conditions. Finally, we show that the spectral bundle method combined with a recent matrix sketching technique is able to solve an SDP with billions of decision variables in a matter of minutes.
Optimization and Control (math.OC)
Revisiting Spectral Bundle Methods: Primal-dual (Sub)linear Convergence Rates
Qwen/Qwen-72B
1. YES
2. YES