\section{Introduction}
\label{sec:intro}

\index{\textbf{SCP}, side-chain positioning}
\index{\textbf{IQP}, integer quadratic problem}
\index{\textbf{DNN}, doubly nonnegative relaxation}
The \emph{protein} \textdef{side-chain positioning (\textbf{SCP})} problem is one of the most important subproblems of the protein structure prediction problem.
We use a doubly nonnegative, \textbf{DNN}, relaxation and apply
a variation of the \textdef{Peaceman-Rachford splitting method (\textbf{PRSM})}
to \emph{exactly} solve an integer formulation of this NP-hard problem.

The applications of \textbf{SCP}\, extend to ligand binding \cite{LoogerLorenL2003Cdor, LAUDET200237} and protein-protein docking with backbone flexibility \cite{WangChu2007PDwB, MarzeNicholasA2018Efbp}.
A protein is a macromolecule consisting of a long main chain backbone that provides a set of anchors for a sequence of amino acid side-chains. The backbone is comprised of a repeating triplet of atoms (nitrogen, carbon, carbon) with the central carbon atom being designated as the alpha carbon.  An amino acid side-chain is a smaller (1 to 18 atoms) side branch that is anchored to an alpha carbon.  The positions of the atoms in a side-chain can be established by knowing the 3D position of its alpha carbon and the dihedral angles defined by atoms in the side-chain.  The number of dihedral angles varies from 1 to 4 depending on the length of the side-chain. This is true for 18 of the 20 amino acids with glycine and alanine being exceptions because their low atom counts preclude dihedral angles.

It has been observed that the values of dihedral angles are not uniformly distributed. They tend to form clusters with cluster centers that are equally separated (+60, 180, -60). Consequently, if the dihedral angles are unknown we at least have a reasonable estimate of their values by appealing to these discretized values.  With this strategy being applied, a side-chain with one dihedral angle would have three possible sets of positions for its atoms. We refer to each set of atomic positions as a rotamer.  A side-chain with two dihedral angles will have 3 times 3 or 9 different arrangements of the atoms (i.e. 9 rotamers).  Three dihedral angles will result in 27 rotamers and four dihedral angles will give 81 rotamers.

In the \textbf{SCP}\, problem we are given a fixed backbone and a designation of the amino acid type for each alpha carbon.  To solve the problem it is required that each amino acid is assigned a particular rotameric setting with the objective of avoiding any collisions with neighbouring amino acids that are given their rotameric settings. Avoiding collisions will lower the overall energy of the protein and, in fact, even with all possible collisions circumvented we want to have an energy evaluation that is minimal.

The \textbf{SCP}\, problem has been proven to be NP-hard \cite{Akutsu1997180}.
The nature of the \textbf{SCP}\, problem has motivated the development of many heuristic based algorithms \cite{CanutescuAdrianA2003Agaf, DesmetJohan1992Tdet, Bahadur2004ProteinSP, SamudralaR1998Dosc, BowerMichaelJ1997Pops, XuJinbo2006Faaa} and many of these approaches rely on the graph structure of the problem.
Other approaches for solving \textbf{SCP}\, problems have been proposed.
These range from probabilistic approaches \cite{LeeChristopher1994PPME, HolmLiisa1991Dafg, ShenkinPeterS1996Paeo},
integer programming \cite{MR1942452, KingsfordCarletonL2005Saas,
AlthausErnst2002Acat}, to
semidefinite programming \cite{MR2098320,ForbesVrisWolk:11}.
Our approach is based on a semidefinite programming relaxation.
Given a rotamer library, the \textbf{SCP}\, problem can be formulated as an
\textdef{integer quadratic problem (\textbf{IQP})}.
We then obtain a \textdef{semidefinite programming (\textbf{SDP})} relaxation
to the \textbf{IQP } via a lifting of variables and \textdef{facial reduction (\textbf{FR})}.
We finally obtain a \textdef{doubly nonnegative (\textbf{DNN})} relaxation
by adding nonnegativity constraints and  some additional constraints to the \textbf{SDP}\, that help strengthen our relaxation.

The facial reduction originating from the \textbf{SDP}\, relaxation delivers a natural splitting of variables.
This elegant splitting of variables fits into the framework of the splitting methods.
The framework gives an efficient procedure of engaging constraints that
are hard to process simultaneously, see
e.g.,~\cite{GHILW:20,LiPongWolk:19,OliveiraWolkXu:15}.
We solve the \textbf{DNN}\, relaxation using a variation of the so-called \textdef{Peaceman-Rachford splitting method (\textbf{PRSM})}.
Using the \textbf{PRSM}, we examine the strength of our approach in the numerical experiments.
\index{\textbf{DNN}, doubly nonnegative relaxation}
\index{\textbf{FR}, facial reduction}
\index{\textbf{SDP}, semidefinite programming}

\subsection{Notation}

We let $\mathbb{R}^n, \mathbb{R}^{m \times n}$ denote the standard real Euclidean spaces;
$\mathbb{S}^n$ denotes the Euclidean space of $n$-by-$n$ real symmetric matrices;
$\mathbb{S}_+^n$ ($\mathbb{S}_{++}^n$, resp) denotes the cone of $n$-by-$n$ positive semidefinite (definite, resp) matrices.
We write $X\succeq 0$ if $X\in \mathbb{S}_+^n$ and $X\succ 0$ if $X\in \mathbb{S}_{++}^n$.
Given $X \in \Rnn$, we use $\trace (X)$ to denote the trace of $X$, and
$\diag(X)$ to denote the vector formed from the diagonal entries of $X$.
Then $\Diag(v) = \diag^*(v)$ is the adjoint linear transformation that
forms the diagonal matrix from the vector $v$.
We use $\range(X)$ and $\Null(X)$ to denote the range of $X$ and the null space of $X$, respectively.
Given two matrices $X,Y \in \mathbb{R}^{m \times n}$, $X \circ Y$ denotes the
element-wise, or Hadamard, product of $X$ and $Y$.
We use $\<X,Y\> = \trace (XY^T)$ to denote the usual trace inner product.
We let $\bar{e}_{n}$ denote the $n$-dimensional vector with each entry set to $1$ and let $\bar{E}_{n}$ denote the $n$-by-$n$ matrix of all ones. We omit the subscript when the meaning is clear.

\subsection{Contributions and Outline}

We begin by presenting the model formulation in~\Cref{sec:Model_Derivation}.
We formulate the \textbf{SCP}\, problem as an integer quadratic program, \textbf{IQP},
and obtain the \textbf{SDP}\, and \textbf{DNN}\, relaxations.
The derivation for the \textbf{SDP}\, relaxation is first presented in~\cite{ForbesVrisWolk:11} via Lagrangian relaxation.
Here, we present a much simpler formulation for the \textbf{SDP}\, relaxation via direct lifting of the variables.
In \Cref{sec:Algorithm} we present a variation of the \textbf{PRSM } for solving the \textbf{DNN}\, relaxation as well as our strategies to obtain upper and lower bounds to the \textbf{SCP}\, problem.
In~\Cref{sec:Numerical_Experiments} we examine the strength of our approach by presenting the performance using the real world data from the Protein Data Bank\footnote{\url{https://www.rcsb.org/}}.
The numerical experiments demonstrate that our approach \emph{provably}
solves all our instances to the global optimum of the hard \textbf{SCP}.
Moreover, we significantly reduce the computation times
in~\cite{ForbesVrisWolk:11}.

\section{Model Derivation}
\label{sec:Model_Derivation}

The goal of this section is to obtain the \textbf{DNN}\, relaxation.
We start by presenting a mathematical formulation of the \textbf{SCP}\,
as an \textbf{IQP } in~\Cref{sec:Problem_Formulation}.  We then
derive its \textbf{SDP}\, relaxation in~\Cref{sec:SDP_Relaxation}.
We continue with identifying redundant constraints in \textbf{IQP } and in the
\textbf{SDP}\, relaxation in order to obtain a complete (stable) model.

\subsection{Problem Formulation as \textbf{IQP}}
\label{sec:Problem_Formulation}

\index{${\mathcal V}  _i$}
\index{side-chain positioning, \textbf{SCP}}
\index{\textbf{SCP}\,, side-chain positioning}
\index{$m_i$}
\index{$n_0$, total number of rotamers}

We are given a collection of disjoint sets ${\mathcal V}  _i, i = 1,\ldots,p$.
Each set ${\mathcal V}  _i$ has $m_i$ members, $|{\mathcal V}  _i|=  m_i$, with total
$n_0 = \sum_{i=1}^p m_i$ and \textdef{${\mathcal V}  = \cup_{i=1}^p {\mathcal V} _i$}.
We call each set ${\mathcal V} _i$ a rotamer set and its members are rotamers.
The \textdef{protein side-chain positioning problem} seeks to select \emph{exactly one} rotamer $v_{i}$ from each set ${\mathcal V} _{i}$, in order to minimize the
sum of the weights (energy) on the edges between chosen rotamers, and
the energy between each chosen rotamer and the
backbone,~see~\Cref{fig:problem_illustration}\footnote{
${\mathcal V} _i$ indicates the $i$-th rotamer set and $v_i^j$ indicates the $j$-th candidate in the $i$-th rotamer set ${\mathcal V} _i$.}.
\begin{figure}[h!]
\centering
\hspace{0.8cm}
\includegraphics[trim= 0 110 0 100,clip,height=5cm]{problem_picture_ver2.eps}
\vspace{-0.3cm}
\caption{A Diagram of the Protein Side-Chain Positioning Problem}
\label{fig:problem_illustration}
\end{figure}
We denote the edge weights between two distinct rotamers (nodes) $u\neq v$ by
the matrix entries $E_{uv}$; while the diagonal entries
$E_{uu}$ denote the weight between the rotamer $u$ and the
backbone. This yields a symmetric matrix $E$, where $E_{uv}=\infty$ if both
rotamers $u,v$ are in the same set. We note that the multiplication
 $0\infty = 0$ when adding up the weights (energies).
Alternatively, we can set these weights to $0$ and add a constraint to
choose exactly one rotamer from each set, which is what we do.
Thus each diagonal block of $E$, of size $m_{i}$, can be assumed to be
a diagonal matrix. We can make this simplification without loss of generality since we are looking to only choose one rotamer per set ${\mathcal V} _{i}$.
\index{$E$, edge weights}
\index{edge weights, $E$}
\index{${\mathcal V} $}

We are looking to solve the following binary quadratic program over
the indicator vector $x$:
\begin{equation}
\label{prob:original form}
\begin{array}{rl}
\min \  & \sum\limits_{u,v} E_{uv} x_{u} x_{v} \\
\text{s.t.} \ & \sum\limits_{u \in {\mathcal V}  _{k}} x_{u} = 1, \ \ k= 1, \dots, p \\
& x=(x_u)\in \{0,1\}^{n_0}.
\end{array}
\end{equation}
We use the block diagonal MATLAB function
\index{$A$}
\index{$n_0$, total number of rotamers}
\index{total number of rotamers, $n_0$}
\begin{equation}
\label{eq:def_A}
A =
\textdef{$\blkdiag$}(  \bar{e}_{m_{1}}^{T} , \bar{e}_{m_{2}}^{T} ,\cdots
,\bar{e}_{m_{p}}^{T} )
\in \mathbb{R}^{p\times n_0},
\end{equation}
to model the \textdef{row sum constraints} that yield exactly one
rotamer is chosen for each set.
We can rewrite the program \eqref{prob:original form} as follows:
\begin{equation}
\label{problem:IQP}
(\textdef{\textbf{IQP}}) \qquad
\begin{array}{ccl}
p^*_{\textbf{IQP}} := &  \min \ & x^{T}Ex\\
&\text{ s.t. }  & Ax = \bar{e}_{p}  \\
&& x = \begin{bmatrix} v_{1}^{T}  & v_{2}^{T} & \dots & v_{p}^{T} \end{bmatrix}^T  \in \{ 0,1\}^{n_{0}} \\
&& v_{i} \in \{ 0,1 \}^{m_{i}}, \ i = 1,\dots, p.
\end{array}
\end{equation}
\index{$p^*_{\textbf{IQP}}$, optimal value of \textbf{IQP }}
\index{optimal value of \textbf{IQP }, $p^*_{\textbf{IQP}}$}

\subsection{\textbf{SDP}\, Relaxation}
\label{sec:SDP_Relaxation}

Given $x\in \mathbb{R}^{n_0}$, we lift to symmetric matrix space
using the rank-one \emph{lifted matrix}
\index{$Y_x$}
\[
Y_x
:= \begin{bmatrix} 1\\x  \end{bmatrix}
\begin{bmatrix} 1\\x  \end{bmatrix}^T
= \begin{bmatrix} 1 & x^T \\x & xx^T  \end{bmatrix} \in \mathbb{S}^{n_0+1}.
\]
In this section we obtain the following \textbf{SDP}\, relaxation to the discrete
optimization problem in~\eqref{problem:IQP}:
\begin{equation}
(\textdef{\textbf{SDP}}) \qquad
\label{eq:SDP_FR_orig}
\begin{array}{rl}
\min\limits_{R,Y} & \trace(\hat{E}Y) \\
& G_{\hat {\mathcal J} }(Y) = E_{00} \\
& Y = VRV^T \\
& R \in \mathbb{S}^{n_0+1-p}_{+} ,  \\
\end{array}
\end{equation}
where
\index{$\hat{E}$}
 \index{$E_{00}$}
\[
\hat{E} = \blkdiag(0,E)  \in \mathbb{S}^{(n_0+1)}, \ \
E_{00} = e_0 e_0^T \in \mathbb{S}^{n_0+1},
\]
\index{$e_0$, the first unit vector}
\index{the first unit vector, $e_0$}
$e_0$ is the first unit vector, and the (gangster) linear operator
$G_{\hat {\mathcal J} }$ and matrix $V$ are defined below.
A variant of the relaxation \eqref{eq:SDP_FR_orig} is proposed
by~\cite{ForbesVrisWolk:11} via
Lagrangian relaxation to \eqref{problem:IQP}.
Here we present a simpler derivation of the model \eqref{eq:SDP_FR_orig} via
the standard direct lifting.

\subsubsection{Gangster Constraint $G_{\hat {\mathcal J} }(Y) = E_{00}$}
\label{sec:GangsterConst}

Given a matrix $Y \in \mathbb{S}^{n_0}$, we define the set of indices
\index{${\mathcal J} $}
\[
\hspace{-0.3cm}
{\mathcal J}  : =
\left\{
\left( \sum_{i=1}^j{m}_{i-1}+ k  ,\ \sum_{i=1}^j{m}_{i-1} + \ell \right) :
j \in \{1,\ldots,p-1\}, \  k,\ell \in \{ 2,\ldots,m_i-1 \}, \ k \ne \ell
\right\} .
\]
Here, $m_i$ is the cardinality of rotamer set ${\mathcal V}  _i$,
and $m_0 = 0$.
In other words,
${\mathcal J} $ is the set of off-diagonal indices of the $m_i$-by-$m_i$ diagonal blocks of $Y\in \mathbb{S}^{n_0}$.
Note that these indices correspond to exactly
\[
Y_{uv}=x_ux_v=0, u\neq v, u,v\in {\mathcal V} _i,
\]
i.e.,~the constraint on two distinct rotamers in the same rotamer set.

With the above set of indices, by abuse of notation, we view and
define the mapping
\index{$G_{{\mathcal J} }$, gangster operator}
\index{gangster operator, $G_{{\mathcal J} }$}
\[
G_{{\mathcal J} } : \mathbb{S}^{n_0} \to \mathbb{R}^{|{\mathcal J} |}  \ \ \text{ by }
G_{{\mathcal J} }(Y) = \left( Y_{ij} \right)_{ij\in {\mathcal J} } .
\]
The map $G_{{\mathcal J} }$ can also be viewed as the operator on $\mathbb{S}^{n_0}$
defined by
$G_{{\mathcal J} }(Y) = (A^T A - I) \circ Y $ with $A$ defined in \eqref{eq:def_A}.
Recall $\circ$ is the Hadamard, element-wise, matrix product.
In plain words, $G_{{\mathcal J} }(Y)$ is the projection that
chooses elements of $Y$ corresponding to the index set ${\mathcal J} $. The term
\emph{gangster} refers to the fact that the constraint
$G_{{\mathcal J} }(Y)=0$ sets (shoots) these elements to be zero.

We now define the set of indices
\index{$\hat{{\mathcal J} }$}
\[
\hat{{\mathcal J} } := \{(0,0)\} \cup {\mathcal J} .
\]
We define the analogous gangster operator $G_{\hat{{\mathcal J} }} $ with $\hat{{\mathcal J} }$:
\index{$G_{\hat{{\mathcal J} }} $, gangster operator}
\index{gangster operator, $G_{\hat{{\mathcal J} }} $}
\[
G_{\hat{{\mathcal J} }} : \mathbb{S}^{n_0+1} \to \mathbb{R}^{|\hat{{\mathcal J} }|}  \ \ \text{ by } G_{\hat{{\mathcal J} }}(Y) = \left( Y_{ij} \right)_{ij\in \hat{{\mathcal J} }} .
\]
This yields the \textdef{gangster constraint} in projection and operator
equivalent forms, respectively,
\[
G_{\hat{{\mathcal J} }} (Y) = e_0, \quad G_{\hat{{\mathcal J} }} (Y) = E_{00}.
\]

\subsubsection{Facial Reduction}
\label{sec:FRconst}

We now derive the constraint $Y = VRV^T$ where $ R \in
\mathbb{S}^{n_0+1-p}_{+}$. If $x$ satisfies the linear equalities, then
\index{$K$}
\[
\begin{array}{rl}
Ax = \bar{e}_p
& \implies \begin{bmatrix} 1\\x \end{bmatrix}^T \begin{bmatrix} - \bar{e}_{p}^T\\ A^T\end{bmatrix} = 0 \\
& \implies
\begin{bmatrix} 1\\x \end{bmatrix} \begin{bmatrix} 1\\x \end{bmatrix}^T \begin{bmatrix} - \bar{e}_{p}^T\{\mathcal A}^T\end{bmatrix} \begin{bmatrix} -\bar{e}_{p}^T\{\mathcal A}^T\end{bmatrix}^T  = 0 \\
& \implies \left\<
\underbrace{
\begin{bmatrix} 1\\x \end{bmatrix} \begin{bmatrix} 1\\x \end{bmatrix}^T }_{=:Y_x} \ ,
\underbrace{
\begin{bmatrix} - \bar{e}_{p}^T\{\mathcal A}^T\end{bmatrix} \begin{bmatrix} -\bar{e}_p^T \\ A^T\end{bmatrix}^T }_{=:K}  \right\rangle =0.
\end{array}
\]
Since both arguments in the last inner product are positive semidefinite,
we have $KY_x =0$ and $K$ provides an \textdef{exposing vector} for the feasible
set, i.e.,~$Y$ feasible implies that $\range(Y) \subseteq \Null(K)$.
Therefore, we find a full column rank matrix
$V \in \mathbb{R}^{(n_0+1)\times (n_0+1-p) }$ such that
\[
\range(V) = \Null\left(
 \begin{bmatrix} -\bar{e}_p^T \\ A^T\end{bmatrix}^T \right)
\qquad (=\Null(K)=\range(Y)).
\]
Since $A$ is full row rank, we get that $\rank(K) = p$ and $Y$ feasible
implies
\[
Y \in V  \mathbb{S}^{n_0+1-p}_{+} V^T.
\]
This is the well-known facial reduction technique,
e.g.,~\cite{DrusWolk:16}.
In particular, for our purposes we choose $V$ with normalized columns.

We note that $\rank( Y_x)=1$.
Dropping the rank restriction on the variable $Y_x$, we have the
constraints $Y = VRV^T$ and $ R \succeq 0$; and this completes the
derivation of the relaxation in \eqref{eq:SDP_FR_orig}.

\subsection{\textbf{DNN}\, Relaxation}
\label{sec:DNN_Relaxation}

In \Cref{sec:SDP_Relaxation} we presented a \textbf{SDP}\, relaxation to \eqref{problem:IQP} via direct lifting.
In this section we complete our relasxation by adding additional constraints to \eqref{eq:SDP_FR_orig}.

The theorem below illustrates two additional properties of the model \eqref{eq:SDP_FR_orig}.
\begin{thm}
\label{thm:arrow_redundancy}
Suppose that
\[
Y \in \mathbb{S}^{n_0+1},\, R \in \mathbb{S}^{n_0+1-p}_{+} \text{  with  }
Y=VRV^T, \, G_{\hat{{\mathcal J} }}(Y) = E_{00}.
\]
Then the following hold.
\begin{enumerate}
\item The first column of $Y$ is equal to the diagonal of $Y$.
\item $\tr (R) = 1+p$.
\end{enumerate}
\end{thm}

\begin{proof}
We recall that $\range(V) = \Null (K) = \Null \left( \begin{bmatrix} -\bar{e}_p & A \end{bmatrix} \right)$.
Hence we have
\[
\begin{bmatrix} -\bar{e}_p & A \end{bmatrix} Y
= \begin{bmatrix} -\bar{e}_p & A \end{bmatrix} VRV^T = 0 RV^T = 0.
\]
We then exploit the structure of
$\begin{bmatrix} -\bar{e}_p & A \end{bmatrix} Y$.
We first partition $Y$ as follows:
\begin{equation}
\label{eq:block_notation}
Y =
\begin{bmatrix}
1 & Y_{10}^T & Y_{20}^T & \cdots & Y_{p0}^T  \\
Y_{10} & Y_{11}   & Y_{12}   & \cdots & Y_{1p} \\
\vdots & \vdots   & \vdots   &  \vdots & \vdots \\
Y_{p0} & Y_{p1}   & Y_{p2}   & \cdots & Y_{pp} \\
\end{bmatrix}  \in \mathbb{S}^{n_0+1} ,
\end{equation}
where $Y_{ii} \in \mathbb{S}^{m_i}$, $Y_{ij} \in \mathbb{R}^{m_i \times m_j}$, $Y_{i0} \in \mathbb{R}^{m_i}, \ \forall i,j \in \{1,\ldots,p\}$.
We use $ Y_{ij}^{\textbf{col}\, \ell}$ to denote the $\ell$-th column of the
$(i,j)$-th block of $Y$ and $Y_{i0,\ell}$ to denote the
$\ell$-th coordinate of the vector $Y_{i0}$.
\index{$ Y_{ij}^{\textbf{col}\, \ell}$}
Then expanding $\begin{bmatrix} -\bar{e}_p & A \end{bmatrix} Y$ with the block representation \eqref{eq:block_notation} gives
\[
\begin{bmatrix} -\bar{e}_p & A \end{bmatrix} Y
=
\begin{bmatrix} a_0 & A_1 & \cdots & A_p \end{bmatrix} \in \mathbb{R}^{p\times (n_0+1)},
\]
where
\[
a_0 =
\begin{bmatrix}
-1 + \bar{e}^T Y_{10} \\
-1 + \bar{e}^T Y_{20} \\
\cdots \\
-1 + \bar{e}^T Y_{p0} \\
\end{bmatrix}  \in \mathbb{R}^{p} ,
\]
and, for each $i\in \{1,\ldots, p\}$,
\[
A_i =
\begin{bmatrix}
-Y_{i0,1} + \bar{e}^T Y_{1i}^{\textbf{col}\, 1} & -Y_{i0,2} + \bar{e}^T Y_{1i}^{\textbf{col}\, 2}
& \cdots & -Y_{i0,m_i} + \bar{e}^T Y_{1i}^{\textbf{col}\, m_i}  \\
\vdots &\vdots &\ddots &\vdots  \\
-Y_{i0,1} + \bar{e}^T Y_{ii}^{\textbf{col}\, 1} & -Y_{i0,2} + \bar{e}^T Y_{ii}^{\textbf{col}\, 2}
& \cdots & -Y_{i0,m_i} + \bar{e}^T Y_{ii}^{\textbf{col}\, m_i}  \\
\vdots &\vdots &\ddots &\vdots  \\
-Y_{i0,1} + \bar{e}^T Y_{1p}^{\textbf{col}\, 1} & -Y_{i0,2} + \bar{e}^T Y_{1p}^{\textbf{col}\, 2}
& \cdots & -Y_{i0,m_i} + \bar{e}^T Y_{1p}^{\textbf{col}\, m_i}  \\
\end{bmatrix}  \in \mathbb{R}^{ p \times m_i}.
\]
Setting $A_i =0, \ \forall i\in \{1,\ldots, p\}$, we see that
\[
Y_{i0,j} = \bar{e}^T Y_{ii}^{\textbf{col}\, j}, \ \forall i \in \{1,\ldots , p \}, \ j \in \{1,\ldots , m_i\}.
\]
Since $G_{\hat{{\mathcal J} }}(Y) = E_{00}$ holds, we see that
\[
\diag(Y_{ii}) = Y_{i0}, \ \forall i \in \{1,\ldots, p\}.
\]
Therefore, we conclude that the first column and the diagonal of $Y$ are identical.

We now show that $\tr (R) = 1+p$.
Setting $a_0 = 0$ gives
\[
\bar{e}^TY_{i0} = 1, \ \forall i \in \{1,\ldots, p\}.
\]
Since $\diag(Y_{ii}) = Y_{i0}, \ \forall i \in \{1,\ldots, p\}$, we must have that $\trace(Y_{ii}) = 1$, $\forall i \in \{1,\ldots, p\}$.
Hence $Y=VRV^T$ gives
\[
1+p = \tr (Y) = \tr(VRV^T) = \tr (R),
\]
where the last equality holds since $V^TV = I$.
\end{proof}

We recall that the original model \eqref{problem:IQP} has the binary constraint on its variable $x$.
We also recall that the direct lifting yields a variable of the form
$\begin{bmatrix} 1 & x^T \\x & xx^T  \end{bmatrix} \in \mathbb{S}^{n_0+1}$.
Hence, we may strengthen our model by including the constraint $ Y_{i,j} \in [0,1], \ \forall i,j$.

We define the sets
\index{${\mathcal Y} $}
\index{${\mathcal R} $}
\[
{\mathcal Y}  := \{ Y \in \mathbb{S}^{n_0+1} :  G_{\hat{{\mathcal J} }}(Y) = E_{00}, \ 0\le Y \le 1 \}, \]
\[
{\mathcal R}  := \{ R \in \mathbb{S}^{n_0+1-p} : R\succeq 0, \ \tr(R) = p+1 \}.
\]
By including additional constraints $\tr (R) = 1+p$ and $0\le Y \le 1$ to the model \eqref{eq:SDP_FR_orig}, we complete our model, \textbf{DNN}\, \  relaxation to \eqref{problem:IQP}:
\index{optimal value of \textbf{DNN}\, relaxation, $p^*_{\textbf{DNN}}$}
\index{$p^*_{\textbf{DNN}}$, optimal value of \textbf{DNN}\, relaxation}
\begin{equation}
\label{eq:SDP_FR}
(\textdef{\textbf{DNN}}) \qquad
\begin{array}{rlc}
p^*_{\textbf{DNN}} = \min\limits_{R,Y} & \tr (\hat{E}Y) \\
& Y = VRV^T   \\
& Y\in {\mathcal Y}  \\
& R \in {\mathcal R}  .\\
\end{array}
\end{equation}
The \textbf{DNN}\, relaxation has a linear objective with an \emph{onto} linear equality
constraint, and compact, convex, feasible set constraints.
The first-order optimality conditions for \eqref{eq:SDP_FR} are
\begin{equation}
\label{eq:optcondition}
\begin{array}{lll}
	0&\in  -V^TZV+{\mathbb N}_{{\mathcal R} }(R),
            &\quad \text{(dual $R$ feasibility)} \\
	0&\in  \hat{E}+Z+{\mathcal N} _{{\mathcal Y} }(Y),
            &\quad \text{(dual $Y$ feasibility)} \\
	Y & =   \widehat{V}R\widehat{V}^T,
              \quad R \in {\mathcal R} , \,  Y \in {\mathcal Y} ,
            & \quad \text{(primal feasibility)}
\end{array}
\end{equation}
where
${\mathcal N} _{{\mathcal R} }(R), {\mathcal N} _{{\mathcal Y} }(Y)$ are the \textdef{normal cones} and $Z$ is a Lagrange multiplier associated with the constraint $Y=VRV^T$.
\Cref{thm:EQzero} below states that some elements of the optimal dual multiplier $Z^*$ are known.
\begin{theorem}
\label{thm:EQzero}
Let $(R^*,Y^*)$ be an optimal pair for~\eqref{eq:SDP_FR}, and let
\[
\textdef{${\mathcal Z} _A$} :=\left\{Z\in\mathbb{S}^{n_0+1}: Z_{i,i} = -(\hat{E})_{i,i},
\ Z_{0,i}=Z_{i,0}= -(\hat{E})_{0,i} , \ i=1,\ldots ,n_0 \right\}.
\]
Then there exists $Z^*\in{\mathcal Z} _A$ such that $(R^*,Y^*,Z^*)$
solves \cref{eq:optcondition}.
\end{theorem}
\begin{proof}
The proof uses the optimality conditions \eqref{eq:optcondition} and  \Cref{thm:arrow_redundancy}.
The proof can be found in~\cite[Theorem 2.11]{GHILW:20}.
\end{proof}

\section{The Algorithm}
\label{sec:Algorithm}

In this section we present the algorithm for solving the problem \eqref{eq:SDP_FR}.
For $\beta >0$, we define the augmented Lagrangian ${\mathcal L} _A$ of the model \eqref{eq:SDP_FR}:
\index{${\mathcal L} _A(R,Y,Z)$, augmented Lagrangian}
\index{augmented Lagrangian, ${\mathcal L} _A(R,Y,Z)$}
\begin{equation}
\label{eq:aug_Lagrangian}
{\mathcal L} _A(R,Y,Z) :=
\<\hat{E},Y\> + \<Z,Y-VRV^T\> + \frac{\beta}{2}\|Y-VRV^T\|_F^2.
\end{equation}
We define the projection operator ${\mathcal P}_{{\mathcal Z} _A}(Z)$ onto the set ${\mathcal Z} _A$, where ${\mathcal Z} _A$ is defined in \Cref{thm:EQzero}.
\index{${\mathcal P}_{\mathcal{Z}_A}(Z)$}
In other words, the projection operator ${\mathcal P}_{{\mathcal Z} _A}(Z)$ sets the  first column and the first row of $Z$ to be 0, the diagonal elements of $Z$ to be the diagonal of $\hat{E}$, except for the $(0,0)$-th entry.

\index{\textbf{PRSM}, Peaceman-Rachford splitting method}
We use \textbf{rPRSM}, a variation of the strictly contractive
Peaceman-Rachford splitting method (\textbf{PRSM}) to solve the model
\eqref{eq:SDP_FR} and is presented in \Cref{algo:PRSM_algorithm}.
\begin{algorithm}[h!]
\caption{\textbf{rPRSM } \cite{GHILW:20} for solving \eqref{eq:SDP_FR} }
\label{algo:PRSM_algorithm}
\begin{algorithmic}
\STATE \textbf{Initialize:} $Y^0 \in \mathbb{S}^{n_0+1}$, $Z^0\in \mathcal{Z}_A$, $\beta\in [1,\infty), \gamma \in (0,1)$
\WHILE {termination criteria are not met}
\STATE $R^{k+1} = \mathop{\rm argmin}\limits_{R\in {\mathcal R} } {\mathcal L} _A (R,Y^k,Z^k)$
\STATE $Z^{k+\frac{1}{2}} = Z^{k} + \gamma  \beta \cdot {\mathcal P}_{\mathcal{Z}_A}\left( Y^{k}- V R^{k+1} V^T \right)$
\STATE $Y^{k+1} = \mathop{\rm argmin}\limits_{Y\in \mathcal{Y} } {\mathcal L} _A (R^{k+1},Y,Z^{k+\frac{1}{2}}) $
\STATE $Z^{k+1} = Z^{k+\frac{1}{2}} + \gamma  \beta \cdot {\mathcal P}_{\mathcal{Z}_A}\left( Y^{k+1}- V R^{k+1} V^T \right)$
\ENDWHILE
\end{algorithmic}
\end{algorithm}
We note that the ordinary \textbf{PRSM } updates the dual multipliers without the projection operator, i.e., ${\mathcal P}_{\mathcal{Z}_0} = I$.
We leave the details of the convergence proof of \textbf{rPRSM } scheme to~\cite[Theorem 3.2]{GHILW:20}.

\begin{remark}
The projection on the dual multipliers $Z$ is motivated from an endeavour to have better dual multipliers at each iteration.
We recall that some of the elements of the optimal dual multipliers are known by \Cref{thm:EQzero}.
The algorithm fixes these known elements to be the optimal elements at every iteration.
\end{remark}

\begin{remark}
The model \eqref{eq:SDP_FR} can be solved by using a standard \textbf{SDP}\, solver.
\cite{ForbesVrisWolk:11} engages the nonnegativity of each element of $Y$ using cutting planes.
However, this approach becomes more computationally challenging as the number of cutting planes increases.
Splitting methods engage the polyhedral constraints $0 \le Y \le 1$ in an economic manner.
We incorporate the positive semidefinite constraint and the nonnegativity constraint very efficiently.
We deal with the positive semidefinite and trace
constraint in the $R$-subproblem,
and then deal with the interval and gangster constraints in the $Y$-subproblem.
\end{remark}

\subsection{Update Formulae}
\label{sec:Update_Formulae}

In this section we present the formulae for the $R$ and $Y$ updates in \Cref{algo:PRSM_algorithm}.
The update rules are discussed in~\cite{GHILW:20}. We include the formulae for the sake of completeness.

\subsubsection{$R$-Update}
\label{sec:R_update}

In this section we present the update rule for the $R$-subproblem.
The formula for the $R$-subproblem, with ${\mathcal L} _A$ defined in \eqref{eq:aug_Lagrangian}, is as follows:
\[
\begin{array}{rcl}
R^{k+1}
& = & \mathop{\rm argmin}\limits_{R\in {\mathcal R} } {\mathcal L} _A(R,Y^{k},Z^{k} ) \\
& = & \mathop{\rm argmin}\limits_{R\in {\mathcal R} } \left\| Y^{k} - VRV^T + \frac{1}{\beta}Z^{k} \right\|_F^2 \\
& = & \mathop{\rm argmin}\limits_{R\in {\mathcal R} } \left\| R - V^T (Y^{k}+\frac{1}{\beta}Z^{k} )V \right\|_F^2 \\
& = & {\mathcal P}_{{\mathcal R} } \left( V^T \left( Y^{k} + \frac{1}{\beta}Z^{k} \right) V \right), \\
& = & U \ \Diag\left( {\mathcal P}_{\Delta_{p+1}} (d) \right) \ U^T ,
\end{array}
\]
where the second equality holds by completing the square and the third
equality holds due to $V^TV = I$; $d$ is the vector of eigenvalues of
$V^T \left( Y^{k} + \frac{1}{\beta}Z^{k} \right) V$ with eigenvector
matrix $U$; \textdef{$\Diag$}
forms the diagonal matrix;
and $ {\mathcal P}_{\Delta_{p+1}} $ is the projection operator onto the simplex $\Delta_{p+1} = \{ x\in \mathbb{R}^{n_0+1-p} : \bar{e}^Tx = 1+p \}$.

\subsubsection{$Y$-Update}
\label{sec:Y_update}
The update rule for $Y$ is as follows.
\begin{equation} \label{Y_update_complete_sqaure}
\begin{array}{rcl}
Y^{k+1}
& = & \mathop{\rm argmin}\limits_{Y\in {\mathcal Y} } {\mathcal L} _A(R^{k+1},Y,Z^{k+\frac{1}{2}} ) \\
& = & \mathop{\rm argmin}\limits_{Y\in {\mathcal Y} } \left\| Y - \left( V R^{k+1} V^T - \frac{1}{\beta} (\hat{E}+Z^{k+\frac{1}{2}} ) \right)  \right\|_F^2 \\
& = & {\mathcal P}_{\text{box}} \left( G_{\hat{{\mathcal J} }^c } \left( V R^{k+1} V^T - \frac{1}{\beta} (\hat{E}+Z^{k+\frac{1}{2}} ) \right)  \right),
\end{array}
\end{equation}
where ${\mathcal P}_{\text{box}}$ is the projection onto the polyhedral set $\{Y
\in \mathbb{S}^{n_0+1} : 0\le Y \le 1\}$. This is done after fixing the
gangster positions.
\index{${\mathcal P}_{\text{box}}$}

\subsection{Bounding}
\label{sec:Bounding}

In this section we present some strategies for obtaining lower and upper bounds to (\textbf{IQP}) in~\eqref{problem:IQP}.

\subsubsection{Lower Bounds from Lagrange Relaxation}
\label{sec:Lower_Bound}

We follow the approaches in \cite{GHILW:20,OliveiraWolkXu:15} and
obtain lower bounds via the Lagrangian dual to our \textbf{DNN}\,
relaxation in \eqref{eq:SDP_FR};
see also~\cite{Eckstein:2020}.
Define the dual functional $g:\mathbb{S}^{n_0+1} \to \mathbb{R}$ by
\index{$g$, dual functional}
\index{dual functional, $g$}
\[
g(Z) := \min\limits_{R\in {\mathcal R} , Y\in {\mathcal Y}  }\<\hat{E},Y\rangle + \< Z,Y-VRV^T\> .
\]
By the linearity of the constraint $Y = VRV^T$, and by the compactness
and convexity of the set constraints ${\mathcal R} \cap {\mathcal Y} $, we note that
\[
\begin{array}{rl}
p^*_{\textbf{DNN}}=\min\limits_{R\in {\mathcal R} , Y\in {\mathcal Y}  } \<\hat{E},Y\rangle + \< Z,Y-VRV^T\>
& = \min\limits_{Y\in \mathcal{Y} } \<\hat{E}+Z,Y\>  + \min\limits_{R\in {\mathcal R} } \<-V^TZV,R\>  \\
& = \min\limits_{Y\in \mathcal{Y} } \<\hat{E}+Z,Y\> - (p+1) \lambda_{\max}(V^TZV),
\end{array}
\]
where \textdef{$\lambda_{\max}$} is the maximum eigenvalue function.
Hence we have a valid lower bound to the optimal value $p^*_{\textbf{DNN}\,}$ of the model \eqref{eq:SDP_FR}:
\begin{equation*}
\label{eq:lbd_computation}
p^*_{\textbf{DNN}\,} = \max\limits_Z  g(Z)
\ge g(Z)
= \min\limits_{Y\in \mathcal{Y} } \<\hat{E}+Z,Y\> - (p+1) \lambda_{\max}(V^TZV),
\end{equation*}
where the first equality holds since the constraint qualification holds for the model \eqref{eq:SDP_FR}.
We note that the computation for $\min\limits_{Y\in \mathcal{Y} } \<\hat{E}+Z,Y\rangle $ can be efficiently performed.

\subsubsection{Upper Bounds from Nearest Binary Feasible Solutions}
\label{sec:Upper_Bound}

In this section we discuss two strategies for obtaining upper bounds to
the \textbf{SCP}\, problem.  These strategies are derived from those presented
in~\cite{ForbesVrisWolk:11} and we include them here for completeness.
We obtain upper bounds by finding feasible solutions to the original
integer model in \eqref{problem:IQP}.
Let $(R^{\text{out}},Y^{\text{out}},Z^{\text{out}})$ be the output of the algorithm.
\begin{enumerate}
\item
Let $x^{\text{approx}} \in \mathbb{R}^{n_0}$ be the second through to the last
elements of the first column of $Y^{\text{out}}$.
Note that $0\leq x^{\text{approx}} \leq 1$.
Then the nearest feasible solution to (\textbf{IQP}) from $x^{\text{approx}}$ can be found by solving the following projection:
\begin{equation}
\label{eq:projection_onto_fea}
\min_x \left\{ \|x - x^{\text{approx}}\|^2 \ : \ Ax = \bar{e}_p, \ x\in \{0,1\}^{n_0}  \right\}.
\end{equation}
It is shown in~\cite{ForbesVrisWolk:11} that
solving \eqref{eq:projection_onto_fea} is equivalent to solving the
following \emph{linear program}:
\begin{equation}
\label{eq:nearest_sol_LP}
\min_x \left\{ \<x , x^{\text{approx}} \> \ : \ Ax = \bar{e}_p, \ x \ge 0  \right\}.
\end{equation}
\item
We now let $x^{\text{approx}}$ be the second through to the last elements of the most dominant eigenvector of $Y^{\text{out}}$.
Note that we again have $0\leq x^{\text{approx}} \leq 1$, by the Perron-Frobenius
theorem. We again obtain the nearest feasible solution to $x^{\text{approx}}$ by
solving \eqref{eq:nearest_sol_LP}.
\end{enumerate}

\begin{remark}
In fact, solving \eqref{eq:nearest_sol_LP} does not require using any LP software; we can obtain the optimal solution for \eqref{eq:nearest_sol_LP} as follows.
We partition $x^{\text{approx}}$ into $p$ subvectors of sizes $m_i=|{\mathcal V} _i|$, for $i=1,\ldots,p$.
Let $x^{i}\in \mathbb{R}^{m_i}$ be the subvector of $x^{\text{approx}}$ associated with $i$-th rotamer set ${\mathcal V} _i$, i.e.,
$x^{\text{approx}} = [ x^{1};x^{2}; \ldots ;x^{p}]$.
We define $\hat{x}^i \in \mathbb{R}^{m_i}$ as follows:
\[
\hat{x}^i_j
= \left\{ \begin{array}{ll}
1 , & \text{if } x^i_j = \max\limits_{\ell \in [m_i] } \left\{x^i_\ell \right\}   \\
0, & \text{otherwise}.\\
\end{array} \right.
\]
If there is subvector $ \hat{x}^{i}$ with more than one $1$ in its components, we pick only one $1$ and set the remaining to be $0$.
We then form $\hat{x} = [\hat{x}^1;\hat{x}^2; \ldots;\hat{x}^p ] \in \mathbb{R}^{n_0}$.
It is clear that $\hat{x}$ is feasible for \eqref{prob:original form}.
We use $\hat{x}^T E \hat{x}$ as an upper bound to the \textbf{SCP}\, problem.
\end{remark}

\section{Numerical Experiments}
\label{sec:Numerical_Experiments}

We present the numerical experiments for \Cref{algo:PRSM_algorithm}.
We use the bounding strategies presented in \Cref{sec:Bounding}
to prove optimality.
We denote the entire process \textbf{rPRSM } where its name follows from \cite{GHILW:20}.
This section is organized as follows.
In \Cref{sec:stopping_Criteria} we present the parameter settings and stopping criteria.
In \Cref{sec:Data} we explain how we process the data from the Protein Data Bank (PDB) to obtain the energy matrix $E$.
In \Cref{sec:Experiments} we finally present the numerical results using
\textbf{rPRSM } and show that we provably solve all the given instances to optimality.
\index{Protein Data Bank, PDB}
\index{PDB, Protein Data Bank}

\subsection{Stopping Criteria and Parameter Settings}
\label{sec:stopping_Criteria}

\paragraph{Stopping Criteria}

We terminate \textbf{rPRSM } when either of the following conditions is satisfied.
\begin{enumerate}
\item
Maximum number of iterations, denoted by ``maxiter'' is achieved.

\item
For given tolerance $\epsilon$, the following bound on the primal and
dual residuals holds for $t$ sequential times:
\[
\max \left\{ \frac{\|Y^k- \widehat{V}R^k \widehat{V}^T\|_F}{\|Y^k\|_F}
	, \beta \|Y^{k}-Y^{k-1}\|_F \right\} <\epsilon.
\]

\item
Let $\{ l_1,\ldots, l_k \}$ and  $\{ u_1,\ldots,u_k \}$
be sequences of lower and upper bounds discussed in \Cref{sec:Lower_Bound} and \Cref{sec:Upper_Bound}, respectively.
Any of the lower bounds achieve the best upper bound, i.e.,
\[
\min \{ l_1,\ldots, l_k \} \ge \max \{ u_1,\ldots,u_k \} .
\]

\end{enumerate}

\paragraph{Parameter Settings}

We use the following parameters related to the implementation of \Cref{algo:PRSM_algorithm}:
\[
\beta = \max\{ \lfloor 0.5*n_0/p \rfloor,1\},  \ \
\gamma = 0.9 .
\]
The parameters related to stopping criteria are:
\[
\text{maxiter} = p(n_0+1) + 10^4 , \ \
\epsilon = 10^{-10}, \ \
t = 100.
\]
For the initial iterates for \textbf{rPRSM}, we use
\[
Y^0= 0 , \ Z^0 = {\mathcal P}_{{\mathcal Z} _A} (Y^0) .
\]

\subsection{Energy Matrix Computation}
\label{sec:Data}

In this section we give a brief description for acquiring the energy matrix $E$.
Our implementation relies on the usage of a Python script executing as an extension of the UCSF Chimera\footnote{The UCSF Chimera software can be found in \url{https://www.cgl.ucsf.edu/chimera/download.html}.} application. A detailed implementation can be found in \cite[Chapter 7]{BurkowskiForbesJ2015CaVT}.
We used protein data files from the PDB to obtain the coordinates of all atoms in the protein.  To get the energy values required by the  algorithm, the native side chain conformations were replaced by rotamers extracted from a rotamer library provided by the Dunbrack Laboratory \cite{Dunbrack341254}.

Some approaches use an energy evaluation based on a piece-wise linear approximation of the Lennard-Jones potential formula
(e.g., \cite{XuJinbo2006Faaa,CanutescuAdrianA2003Agaf}).
Here, we used the Lennard-Jones potential formula, which provides a more accurate energy value computation.
In brief, the Lennard-Jones potential formula engages the Euclidean distance between a pair of atoms with some parameters dependant on the type of amino acids.
A more detailed explanation of these energy computations can be found in \cite[Chapter 6-7]{BurkowskiForbesJ2015CaVT}.
\index{dead end elimination}
We finally used a strategy (known as `dead end elimination') to reduce the size of the rotamer sets associated with each amino acid.
The basic idea behind this strategy is that a rotamer can be removed from its rotamer set if there is another rotamer in that set that gives a better energy value regardless of the rotamer selections for the neighbouring amino acids.
Among various approaches for the dead end elimination, we followed the Goldstein's criteria \cite{Goldstein:94}.

Let $U$ be a side-chain conformation of a protein.
The energy of the conformation $U$ is
\[
E(U) = \sum_{i=1}^{n_0} E_{\text{self}} (u_i) + \sum_{i=1}^{n_0-1}\sum_{j=i+1}^{n_0} E_{\text{pair}} (u_i,u_j) ,
\]
where $u_i$ is a side-chain conformation of an amino acid,
$E_{\text{self}} (u_i)$ is the energy corresponding to $u_i$ and the backbone,
and $E_{\text{pair}} (u_i,u_j)$ is the energy formed by $u_i$ and  $u_j$, a rotamer associated with a neighbouring amino acid.
In our formulation, we placed $E_{\text{self}} (u_i)$ along the diagonal of $E$ and $E_{\text{pair}} (u_i,u_j)$ on the appropriate off-diagonal positions of $E$ as shown in \Cref{sec:Problem_Formulation}.

\subsection{Experiments with Real World Data}
\label{sec:Experiments}

In this section we demonstrate the strength of our approach by presenting the numerical experiment using the 35 instances from the PDB.
All instances in \Cref{table:PDBinstance1} are tested using MATLAB
2018b on a Dell PowerEdge M630 with two Intel Xeon E5-2637v3 4-core 3.5 GHz (Haswell) with 64 Gigabyte memory.
The following list defines the column headers used in Table
\ref{table:PDBinstance1}.
\begin{enumerate}[itemsep=-1mm]
\item {\bf problem}: instance name;
\item $\bm{p}$: the number of amino acids;
\item $\bm{n_0}$: the total number of rotamers;
\item {\bf lbd}: the lower bound obtained by running \textbf{rPRSM};
\item {\bf ubd}: the upper bound obtained by running \textbf{rPRSM};
\item {\bf rel.gap}: relative gap of each instance using \textbf{rPRSM},
where
\begin{equation*}
\label{def:rel_gap}
\text{relative gap} := 2\ \frac{ | \text{best feasible upper
bound}- \text{best lower bound} | }{\ | \text{best feasible upper
bound}+\text{best lower bound}+1 | } ;
\end{equation*}
\item {\bf iter}: number of iterations used by \textbf{rPRSM } with tolerance
	$\epsilon = 10^{-10}$;
	\item {\bf time(sec)}: CPU time (in seconds) used by \textbf{rPRSM}.
\end{enumerate}
\begin{table}[h!]
\small
\centering
\caption{The Numerical Result of the 35 PDB Instances}
\label{table:PDBinstance1}
\vspace{0.15cm}
\begin{tabular}{|l|cc|ccccc|}
\hline \rule{-0pt}{8pt}
\textbf{problem} & $\bm{p}$ & $\bm{n_0}$ & \textbf{lbd} & \textbf{ubd}   & \textbf{rel.gap} & \textbf{iter} & \textbf{time(sec)} \\
\hline
1CRN   & 37 & 130 & -48.46 & -48.46 & 5.55e-13 & 1579 & 8.37 \\
1AAC  & 85 & 441 & -318.97 & -318.97 & 3.57e-15 & 400 & 10.02 \\
1TIM & 192 & 822 & 571.82 & 571.82 & 1.65e-14 & 1800 & 126.00 \\
1AHO & 54 & 140 & -17.06 & -17.06 & 1.50e-12 & 1428 & 8.05 \\
1CC7  & 66 & 498 & -140.03 & -140.03 & 6.11e-16 & 400 & 11.70 \\
1CEX  & 146 & 415 & 200.75 & 200.75 & 8.26e-14 & 13000 & 222.40 \\
1IGD  & 50 & 116 & -114.38 & -114.38 & 1.97e-13 & 1083 & 4.80 \\
1CTJ  & 61 & 242 & -168.87 & -168.87 & 2.29e-13 & 3958 & 40.83 \\
1CZ9  & 111 & 688 & -292.55 & -292.55 & 1.69e-14 & 2800 & 150.14 \\
1CZP  & 83 & 626 & -134.89 & -134.89 & 5.39e-10 & 4200 & 179.98 \\ 1ITM  & 119 & 736 & 131.21 & 131.21 & 2.40e-14 & 7700 & 472.91 \\
1MFM & 118 & 456 & -272.98 & -272.98 & 5.80e-14 & 3200 & 68.52 \\
1PLC  & 82 & 224 & -162.54 & -162.54 & 7.88e-12 & 1200 & 9.58 \\
1QJ4  & 221 & 1121 & -596.60 & -596.60 & 2.52e-13 & 2600 & 307.76\\
1QQ4  & 143 & 851 & -235.68 & -235.68 & 7.77e-13 & 2231 & 224.81 \\
1VFY  & 63 & 173 & -166.32 & -166.32 & 1.63e-11 & 15400 & 138.33 \\
4RXN  & 48 & 161 & -128.13 & -128.13 & 6.46e-15 & 2700 & 16.69 \\
1ARE  & 26 & 260 & -32.93 & -32.93 & 1.12e-14 & 400 & 4.64 \\
1BXL  & 163 & 680 & 175.56 & 175.56 & 7.64e-13 & 3029 & 150.06 \\
1CB3  & 8 & 16 & -7.81 & -7.81 & 2.53e-10 & 800 & 0.13 \\
1EDN & 21 & 25 & 31.79 & 31.79 & 3.32e-14 & 100 & 0.02 \\
1EDP & 17 & 20 & -14.03 & -14.03 & 1.41e-14 & 100 & 0.02 \\
1TZ4 & 30 & 55 & -36.41 & -36.41 & 1.92e-12 & 5100 & 5.14 \\
1TZ5 & 31 & 57 & -25.55 & -25.55 & 1.56e-14 & 300 & 0.23 \\
2KAP  & 53 & 169 & -64.28 & -64.28 & 2.37e-13 & 1677 & 11.50 \\
2KGU  & 29 & 82 & 409.29 & 409.29 & 1.53e-15 & 2541 & 5.74 \\
2ZNF  & 14 & 16 & -28.73 & -28.73 & 2.52e-16 & 100 & 0.02 \\
6PAX  & 114 & 310 & -336.12 & -336.12 & 3.79e-14 & 2200 & 27.41 \\
5P21 & 144 & 892 & -359.18 & -359.18 & 4.20e-14 & 9100 & 781.38 \\
1PCH & 69 & 92 & -120.82 & -120.82 & 2.39e-13 & 1806 & 5.39 \\
2EWH & 63 & 215 & -171.91 & -171.91 & 1.34e-10 & 1700 & 15.70 \\
1AGT & 33 & 50 & 29.41 & 29.41 & 3.11e-14 & 200 & 0.07 \\
1RTE & 100 & 500 & -331.02 & -331.02 & 9.46e-15 & 1900 & 51.57 \\
4PFK & 247 & 975 & -622.01 & -622.01 & 4.13e-13 & 2500 & 235.62 \\
1MBO & 125 & 707 & -216.37 & -216.37 & 1.32e-09 & 1900 & 106.66 \\
\hline
\end{tabular}
\end{table}

\paragraph{Discussion}

We observe from the last two columns of Table \ref{table:PDBinstance1} that all 35 instances are solved within a reasonable amount of time.
The column \textbf{rel.gap} illustrates the relative gaps of the  instances and we observe that these gaps are essentially 0.
We recall from \Cref{sec:Upper_Bound} that we obtain the upper bounds via finding feasible solutions to (\textbf{IQP}).
That we have the relative gap essentially $0$ grants us the attainment of the \emph{globally optimal} solutions to the \textbf{SCP}\, problem.
Approaches involving heuristic algorithms do not provide a natural means of certifying optimality, relying solely on a comparison of the rotameric solution with native $\chi_1$ and $\chi_2$ angles from the PDB while ignoring optimality of the discretized solution.
Here, we highlight that we provide not only the globally optimal solutions but also a way to certify their optimality.

During the course of the energy matrix computation presented in \Cref{sec:Data}, we typically observe some very large elements in $E$.
These large values (generally more than 10 digits) stem from collisions of pairs of distinct rotamers.
Intuitively, a reasonable conformation should not include rotamers that yield very large energy values.
In general, having very large values in data is prone to numerical instabilities.
However, this ill-posed data does not take place as a problem in our implementation.
Recall that we solve the $Y$-subproblem \eqref{Y_update_complete_sqaure} as follows:
\[
Y^{k+1} = {\mathcal P}_{\text{box}} \left( G_{{\mathcal J} ^c } \left( V R^{k+1} V^T - \frac{1}{\beta} (\hat{E}+Z^{k+\frac{1}{2}} ) \right)  \right).
\]
If there is a very large element $(i,j)$ in $\hat{E} = \blkdiag(0,E)$, the operator ${\mathcal P}_{\text{box}} $ sets the element $(i,j)$ of $Y^{k+1}$ to $0$.
Hence, for those positions $(i,j)$ with very large values, the constraint $Y_{i,j}=0$ is implicitly imposed.
We can interpret this as having implicit gangster constraints on these entries.

\section{Conclusions}

We presented a simple way of formulating the relaxation of
the \textbf{SCP}\, problem.
We began by formulating the \textbf{SCP}\, problem into an \textbf{IQP } and derived a
facially reduced \textbf{SDP}\, relaxation.
We then identified some redundant constraints to the \textbf{IQP } to complete a \textbf{DNN}\, relaxation.
\textbf{FR}\, provided a natural splitting of the variables.
Hence we adopted the \textbf{rPRSM } to solve the \textbf{DNN}\, relaxation of the \textbf{SCP}\, problem.
Our approach \emph{solved all $35$ PDB chosen instances to optimality}
thus illustrating the efficiency of this approach for this NP-hard
problem.

\begin{comment}

\section{Recover these after the first submission}

\subsection{Removing Collisions}
Collisions between rotamers are indicated by \emph{huge} energy values
$E_{ij} >> 0$. We can take advantage of these to eliminate edges in
the graph and increase the size of the gangster set.
\index{collisions}
\begin{lemma}[collision elimination]
Suppose that $x$ is feasible for the \textbf{IQP } in \eqref{problem:IQP}, and let
$u = x^TEx$ be its objective value, and thus $p^*\leq u$.
Let $N_E = \sum_{\{ij : E_{ij}<0\}}E_{ij}$ and suppose that
\begin{noteH}
?????we can do better than $N_E$ using only one negative from each set?
or from each node?????
\end{noteH}
\[
2E_{i_0j_0} > u-N_E, \text{  for some   }  i_0,j_0
\]
holds. Then for any optimal solution $x^*$ to \textbf{IQP}, we have $x^*_{i_0}x^*_{j_0} = 0$.
\end{lemma}
\begin{proof}
\index{$U$, set of selected rotamers}
\index{set of selected rotamers, $U$}
Let $x^*$ be an optimal solution to \textbf{IQP}.
Let $U$ be the set of selected rotamers found in optimal solution $x^*$.
Suppose to the contrary that $x^*$ holds $x^*_{i_0}x^*_{j_0} = 1$, i.e., $x^*_{i_0}=x^*_{j_0} = 1$.
Then we reach the following contradiction:
\[
p^*_{\textbf{IQP}} = \<x^*,Ex^*\>  =  2E_{i_0j_0}  + \sum\limits_{i,j\in U\setminus \{ i_0,j_0 \}} E_{i,j} \ \ge \ 2E_{i_0j_0}  +N_E > u .
\]
\end{proof}

\begin{remark}
\begin{noteH}
First we note about choosing $Y$ for the first term, we use ????? see
the code???? roundoff error for $E+Z$???
Can there be roundoff error here from HUGE numbers?

In fact do we need any of the remark?????
\end{noteH}

Second, the calculation $V^TZV$ can yield roundoff error for the
eigenvalue calculation, especially when amplified by $p$. So we would
like to do iterative refinement. Suppose we have
\[
V^TZV u - \lambda u = \delta \neq 0, \|u\|=1,
\]
where $\lambda < 0$ is the approximate \emph{largest
eigenvalue}, i.e.,~$1/|\lambda|$ is the approximation of the
spectral radius of the $(V^TZV)^{-1}$. Suppose that
$Z=-R^TR$ is the Cholesky factorization. Then we need to solve
\[
-V^TR^TRV u^\prime =  u.
\]

The inverse power method can be applied here to improve the accuracy and
reduce $\|\delta\|$. We repeat
\[
u^\prime = V^T(Z(V u)), u^\prime = u^\prime/\|u^\prime\|,
\lambda^\prime = (u^\prime)^TV^T(Z(V u^\prime)). ???? change to inverse
power method???
\]

\url{https://link.springer.com/article/10.1007/s13160-018-0310-3}

????find blas and/or lapack for the multiplication $V'ZV$???????
can we avoid the roundoff error from the symmetrization after the mult??

let $W=V^TZV$ and we find the eigenvalues eigenvectors, $U,D$ and we note that
the first $k$ largest are nearly equal.
Let $W = U(:,1:k)'*W*U(:,1:k)$.

\end{remark}

\end{comment}

\cleardoublepage
\phantomsection
\addcontentsline{toc}{section}{Index}
\printindex
\label{ind:index}

\cleardoublepage
\phantomsection
\addcontentsline{toc}{section}{Bibliography}