\shrink
\section{The Dual View: Fast Optimality Certification}
\label{sec:dualCertification}
\shrink

Despite scaling linearly in $N$, the sparse relaxation~\eqref{eq:sparseMoment} is still too large to be solved efficiently using current interior point methods (IPM)~\cite{Nesterov18book-convexOptimization} when $N > 20$. 
On the other hand, fast heuristics such as graduated non-convexity~\cite{Yang20ral-GNC} can compute globally optimal solutions to the POP~\eqref{eq:pop} with high probability of success. In this section, we show that, by taking the dual perspective of sums-of-squares (SOS) relaxations, we can develop efficient \emph{certifiers} to verify the optimality of a candidate solution $(\hatxextend,\hatf)$ for large $N$ (\eg~$N=100$), for which the SDP relaxation~\eqref{eq:sparseMoment} is not even implementable. 

\subsection{Sums-of-Squares Relaxation}
A candidate solution $(\hatxextend,\hatf)$ is globally optimal for the POP~\eqref{eq:pop} if and only if $f(\xextend) - \hatf \geq 0, \forall \xextend \in \calxextend$. 

However, testing \emph{nonnegativity} of a polynomial on a constraint set is NP-hard~\cite{Blekherman12Book-sdpandConvexAlgebraicGeometry}, so instead we test if the polynomial is SOS on the constraint set and provide a sufficient condition for global optimality.
\begin{theorem}[Sufficient Condition for Global Optimality] \label{thm:sufficientOptimalityCondition}
Given any candidate solution $(\hatxextend,\hatf)$ to the POP~\eqref{eq:pop}, if the following optimization is \emph{feasible} (\ie has at least one solution):
\bea
\text{find} & \vlambda_j^x \in \Real{\dimbasis{\dimxextend}{2}}, \vlambda_j^\theta \in \Real{\dimbasis{n}{2}}, \MS_0 \in \psd^{\dimbasis{}{\rbasisset}}, \MS_k \in \psd^{\dimbasis{\dimxextend}{1}}   \label{eq:sosfeasiblity}\\
 \subject &  \hspace{-3mm} \displaystyle f(\xextend)\! - \! \hatf \! - \!\!\!\! \sum_{h_j \in \vh^x} \!\! h_j \! \left(\monoleq{\xextend}{2}\tran\vlambda_j^x \right) \! - \!\!\!\! \sum_{h_j \in \vh^{\theta} } \!\! h_j \! \left(\monoleq{\vxx}{2}\tran \vlambda_j^\theta \right) \! = \! \monoleq{\xextend}{\rbasisset}\tran \MS_0 \monoleq{\xextend}{\rbasisset} \! + \!\! \sum_{k=1}^{\nrIneq} g_k \! \left(\monoleq{\xextend}{1}\tran \MS_k \monoleq{\xextend}{1} \right) \! ,\edit{ \forall \xextend}, \label{eq:matchingcoefficients}
\eea
then $\hatf$ (resp. $\hatxextend$) is the global minimum (resp. global minimizer) of the POP~\eqref{eq:pop}. Moreover, problem~\eqref{eq:sosfeasiblity} can be written compactly as a feasibility SDP:
\bea
\text{find} \quad \dualVar, \quad   
\subject \quad  \dualVar \in \calK \cap \calA, \label{eq:sdpfeasibility}
\eea 
where $\dualVar = [(\vlambda_1^x)\tran,\dots,(\vlambda^x_{|\vh^x|})\tran,(\vlambda_1^\theta)\tran,\dots,(\vlambda^\theta_{|\vh^\theta|})\tran,\svec{\MS_1}\tran,\dots,\svec{\MS_{\nrIneq}}\tran,\svec{\MS_0}\tran]\tran$ concatenates all variables in~\eqref{eq:sosfeasiblity}, $\calK$ defines a convex cone, and $\calA \doteq \left\{ \dualVar: \MA\dualVar = \vb \right\}$ defines an affine subspace, \edit{where $\vb$ is a vector and $\MA$ is a matrix} satisfying the \emph{partial orthogonality} property~\cite{Zheng18TAC-partialOrthogonality,Bertsimas13OMS-acceleratedSOSRelaxation}.
\end{theorem}


In the \supp, we provide a proof of Theorem~\ref{thm:sufficientOptimalityCondition}. 

Intuitively, 
 if problem~\eqref{eq:sosfeasiblity} is feasible, then for any $\xextend \in \calxextend$, the left-hand side of~\eqref{eq:matchingcoefficients} reduces to $f(\xextend) - \hatf$ (due to $h_j=0$) and the right-hand side of~\eqref{eq:matchingcoefficients} is nonnegative (due to $g_k \geq 0,\MS_0,\MS_k\succeq 0$), producing a certificate that \edit{$f(\xextend) \geq \hatf$}.
 
 The SOS relaxation~\eqref{eq:sosfeasiblity} also uses basis reduction and it is the dual of the sparse moment relaxation~\eqref{eq:sparseMoment}~\cite{lasserre10book-momentsOpt} with the constraint that $\hatf$ is the global optimum. In SDP~\eqref{eq:sdpfeasibility}, the convex cone $\calK$ corresponds to the PSD constraints in~\eqref{eq:sosfeasiblity} and the affine subspace $\calA$ corresponds to matching coefficients in the equality constraint~\eqref{eq:matchingcoefficients}.
 The partial orthogonality of $\MA$ is a property for SDPs resulting from SOS relaxations and allows efficient \emph{projection} onto the affine subspace $\calA$~\cite{Zheng18TAC-partialOrthogonality,Bertsimas13OMS-acceleratedSOSRelaxation}.


\subsection{Douglas-Rachford Splitting}
\label{sec:drs}

In this section, we propose a first-order method based on \emph{Douglas-Rachford Splitting} (\DRS)~\cite{Combettes11book-proximalSplitting,Jegelka13NIPS-DRSreflection} to solve~\eqref{eq:sdpfeasibility} at scale. 
\DRS iteratively solves~\eqref{eq:sdpfeasibility} by starting at an arbitrary initial point $\dualVar_0$, 
and performing the following  three-step updates (at each iteration $\tau \geq 0$):
\bea
(i)\  \dualVar_\tau^{\calK} = \proj_{\calK} \left(\dualVar_{\tau} \right) ,\quad  
(ii)\ \dualVar^{\calA}_{\tau} = \proj_{\calA} \left(2 \dualVar_\tau^{\calK} - \dualVar_{\tau} \right), \quad 
(iii)\ \dualVar_{\tau+1} = \dualVar_{\tau} + \gamma_{\tau} \left( \dualVar^{\calA}_{\tau} - \dualVar_\tau^{\calK} \right), \label{eq:DRSIterates}
\eea 
where $\proj_{\calK}$ (resp. $\proj_{\calA}$) denotes the orthogonal projection onto $\calK$ (resp. $\calA$) 
and $\gamma_{\tau}$ is a parameter of the algorithm.
\optional{
Towards this goal, we first rewrite problem~\eqref{eq:sdpfeasibility} equivalently as $\min_{\dualVar} \indicator_{\calA}(\dualVar) + \indicator_{\calK} (\dualVar)$, where $\indicator_{\calA}(\dualVar)$ is the indicator function for convex set $\calA$ (\ie~$\indicator_{\calA}(\dualVar) = 0$ if $\dualVar \in \calA$ and $\indicator_{\calA}(\dualVar) = \infty$ if $\dualVar \not\in \calA$). Then we apply \DRS and obtain the following result.
}{}
The rationale behind the use of \DRS to solve the feasibility SDP~\eqref{eq:sdpfeasibility} is that, although finding $\dualVar \in \calK \cap \calA$ is expensive (requires solving a large-scale SDP), finding $\dualVar \in \calK$ and $\dualVar \in \calA$ separately (\ie~projecting onto $\calK$ and $\calA$ separately) is computationally inexpensive~\cite{Henrion12handbook-conicProjection,Bauschke96SIAM-projectionFeasibility,Higham88LA-nearestSPD}.
The following 
result shows how to certify optimality using the \DRS iterations~\eqref{eq:DRSIterates}. 
\begin{theorem}[\DRS for Optimality Certification] \label{thm:DRSOptimalityCertification} 
Consider the \DRS iterations~\eqref{eq:DRSIterates}. Then the following properties hold true: 


(i) If the SDP~\eqref{eq:sdpfeasibility} is feasible, then the sequence $\{ \dualVar_\tau \}_{\tau \geq 0}$ in~\eqref{eq:DRSIterates} converges to a solution of~\eqref{eq:sdpfeasibility} when $0 < \gamma_\tau < 2$;

(ii) Let $\subopt = (\hatf - f^\star)/\hatf$ be the relative suboptimality between $\hatf$ and the global minimum $f^\star$ of the POP~\eqref{eq:pop}, then each \DRS iteration~\eqref{eq:DRSIterates} gives a valid suboptimality upper bound $\suboptbound_\tau$, \ie  $\subopt \leq \suboptbound_{\tau}$, and $\suboptbound_{\tau}$ can be efficiently computed from $\dualVar_\tau^\calA$.


\end{theorem}


A complete proof of Theorem~\ref{thm:DRSOptimalityCertification} 

is given in the \supp.

The intuition behind Theorem~\ref{thm:DRSOptimalityCertification}(i) is that, 
by using the two projections alternatively (thus, the name ``\emph{splitting}''), the \DRS iterates~\eqref{eq:DRSIterates} \edit{converge} to a solution in $\calK \cap \calA$ if the intersection is nonempty. Moreover, even if the intersection is empty (\eg~when $\hatf$ is not the global minimum), Theorem~\ref{thm:DRSOptimalityCertification}(ii) states that each \DRS iteration is still able to \emph{assess} the suboptimality of $\hatf$, 
which enables the \emph{dual certifiers} to detect wrong candidate solutions (\cf~Section~\ref{sec:experiments:certification}). 

\DRS converges faster than the vanilla \emph{alternating projections to convex sets} used in~\cite{Yang20arxiv-teaser} (\cf~\cite{Bauschke04JAT-averagedAlternatingReflections}). 
Moreover, we further boost convergence speed by initializing \DRS with an initial point $\vd_0$ computed by solving an
inexpensive SOS program with \emph{chordal sparsity}~\cite{Waki06jopt-SOSSparsity,lasserre10book-momentsOpt} (see the \supp for 
implementation details).


\optional{
In theory, one can start \DRS~\eqref{eq:DRSIterates} at any initial condition $\dualVar_0$. However, to speed up \DRS, we compute the initial point $\vd_0$ by solving a cheap SOS program with \emph{chordal sparsity}~\cite{Waki06jopt-SOSSparsity,lasserre10book-momentsOpt}.

\subsection{Chordal Sparse Initialization}
In theory, one can start \DRS~\eqref{eq:DRSIterates} at any initial condition $\dualVar_0$. However, to speed up \DRS, we compute the initial point $\vd_0$ by solving a cheap SOS program with \emph{chordal sparsity}~\cite{Waki06jopt-SOSSparsity,lasserre10book-momentsOpt}.
\begin{proposition}[Chordal Sparse Initialization] \label{prop:chordalSparseInitialization}
Define $\monoleq{\xextend}{\rbasisset_i} \doteq [1,\vxx\tran,\theta_i, \theta_i \vxx\tran]\tran \in \Real{2n+2}$, $\monoleq{\xextend}{1i} \doteq [1,\vxx\tran,\theta_i]\tran \in \Real{n+2}$, as the sparse monomial bases only in $\vxx$ and $\theta_i,i=1,\dots,N$, and define $\sosindeg{\xextend}{2\rbasisset_i},\sosindeg{\xextend}{2i}$ as the set of SOS polynomials parametrized by the sparse bases $\monoleq{\xextend}{\rbasisset_i},\monoleq{\xextend}{1i}$,~\ie~$q \in \sosindeg{\xextend}{2\rbasisset_i}$ if and only if $q= \monoleq{\xextend}{\rbasisset_i}\tran \MQ \monoleq{\xextend}{\rbasisset_i}$ for some $\MQ \in \psd^{2n+2}$ and $q \in \sosindeg{\xextend}{2i}$ if and only if $q = \monoleq{\xextend}{1i}\tran \MQ \monoleq{\xextend}{1i}$ for some $\MQ \in \psd^{n+2}$. Then $\dualVar_0$ is the solution of the SOS program:
\bea
\max & \zeta \label{eq:soschordalsparse}\\
\subject &  f(\xextend) - \zeta - \sum_{h_j \in \vh^{\vxx}} \lambda_j^{\vxx} h_j  - \sum_{h_j \in \vh^{\vtheta} }  \lambda_j^{\vtheta} h_j = s_0 + \sum_{k=1}^l s_k g_k, \\
& \lambda_j^\vxx, \lambda_j^\vtheta \text{ satisfy eq.~}\eqref{eq:lambdaSparsity}, \\
& s_0 = \sum_{i=1}^N s_{0i}, s_{0i} \in \sosindeg{\xextend}{2\rbasisset_i}, s_k = \sum_{i=1}^N s_{ki}, s_{ki} \in \sosindeg{\xextend}{2i},k=1\dots,l. \label{eq:soschordalsparseCon}
\eea
\end{proposition}
The SOS program~\eqref{eq:soschordalsparse} is the same as the SOS program~\eqref{eq:sosfeasiblity}, except that we have replaced constraint~\eqref{eq:sosSparsity} with constraint~\eqref{eq:soschordalsparseCon}, and we maximize $\zeta$ instead of requiring $\zeta = \hatf$. The advantage of SOS program~\eqref{eq:soschordalsparse} is that it breaks the large PSD constraint of size $\dimbasis{}{\rbasisset}$ (\ie~$s_0 \in \sosindeg{\xextend}{2\rbasisset}$) into $N$ small PSD constraints of fixed size $2n+2$ (\ie~$s_{0i} \in \sosindeg{\xextend}{2\rbasisset_i}$), and can scale to much larger $N$. However, the trade-off is program~\eqref{eq:soschordalsparse} is more restrictive and in general its optimum $\zeta^\star$ cannot certify the global optimality of $\hatf$ (\ie~$\zeta^\star < p^\star_{\rbasisset} \leq f^\star \leq \hatf$). Therefore, we compute $\dualVar_0$ by solving this cheap SOS program~\eqref{eq:soschordalsparse} using IPM-based SDP solvers and refine $\dualVar_0$ by running \DRS for the more powerful SOS program~\eqref{eq:sosfeasiblity}.
}{}
{"config": "arxiv", "file": "2006.06769/dual-certification.tex"}