#  A feasible method for solving an SDP relaxation of the quadratic knapsack
problem

Tianyun Tang , Kim-Chuan Toh Department of Mathematics, National University of
Singapore, Singapore 119076 (ttang@u.nus.edu). Department of Mathematics, and
Institute of Operations Research and Analytics, National University of
Singapore, Singapore 119076 (mattohkc@nus.edu.sg). The research of this author
is supported by the Ministry of Education, Singapore, under its Academic
Research Fund Tier 3 grant call (MOE-2019-T3-1-010).

###### Abstract

In this paper, we consider an SDP relaxation of the quadratic knapsack problem
(QKP). After applying low rank factorization, we get a non-convex problem,
whose feasible region is an algebraic variety with certain good geometric
properties which we analyse. We derive a rank condition under which these two
formulations are equivalent. This rank condition is much weaker than the
classical rank condition if the coefficient matrix has certain special
structures. We also prove that under an appropriate rank condition, the non-
convex problem has no spurious local minima without assuming linearly
independent constraint qualification. We design a feasible method that can
escape from non-optimal non-regular points. Numerical experiments are
conducted to verify the high efficiency and robustness of our algorithm as
compared to other solvers. In particular, our algorithm is able to solve a
one-million dimensional sparse SDP problem accurately in about 20 minutes on a
modest computer.

##  1 Introduction.

###  1.1 An SDP relaxation of quadratic knapsack problem.

The binary quadratic knapsack problem (QKP) was introduced in [20] by Gallo et
al. as follows:

${\rm(QKP)}\ \max\left\\{x^{\top}Cx:\ a^{\top}x\leq\tau,\
x\in\\{0,1\\}^{n}\right\\},$

where $C\in\mathbb{S}^{n}$ is the profit matrix, $a\in\mathbb{R}^{n}$ is the
positive weight vector and $\tau>0$ is the knapsack capacity. We consider the
nontrivial case where $n>1.$ We also make the following natural assumption on
the data $(a,\tau)$.

######  Assumption 1.

The weight vector $a\in\mathbb{R}^{n}$ and capacity $\tau\in\mathbb{R}$
satisfy the condition that for any $i\in[n]$, $0<a_{i}<\tau$ and
$e^{\top}a>\tau,$ where $e$ is the vector of all ones.

The problem (QKP) has many applications such as very large-scale integration
(VLSI) and compiler design [16, 26]. When $C$ is a diagonal matrix, (QKP)
reduces to a binary linear knapsack problem. (QKP) is NP-hard because it can
be reduced to the clique problem (which is closely related to the densest
$k-$subgraph problem) by choosing $\tau=k\in\mathbb{N}^{+},$ $a$ to be the all
ones vector and $C$ to be the adjacency matrix of an undirected graph. Like
many combinatorial optimization problems, (QKP) is extremely difficult to be
solved exactly. Many SDP relaxations are provided in the literature to
calculate an upper bound of (QKP). In [32], Pisinger discussed 4 types of SDP
relaxations of (QKP) whose optimal values are ordered as $U^{1}_{\rm HRW}\geq
U^{2}_{\rm HRW}\geq U^{3}_{\rm HRW}\geq U^{4}_{\rm HRW}$ respectively. These
SDP models were originally provided by Helmberg et al. in [22]. However, the
authors only mentioned the SDP model with optimal value $U^{3}_{\rm HRW}$ for
a special case of (QKP) and it was later generalized by Pisinger in [32].
Among them, the first three models have the same number of constraints and the
last one has multiple inequality constraints. However, it has been shown in
the numerical experiments of [22] that even though the first three SDP models
are simpler than the last one, $U^{2}_{\rm HRW}$ is already quite close to
$U^{4}_{\rm HRW}.$ Therefore, in this paper, we only consider the third SDP
model with both tightness and simplicity. We write it as follows:

${\rm(SQK)}\ \min\left\\{-\left\langle C,X\right\rangle:\ a^{\top}Xa-\tau
a^{\top}x\leq 0,\ {\rm diag}\left(X\right)=x,\ Y:=\begin{pmatrix}1&x^{\top}\\\
x&X\end{pmatrix}\in\mathbb{S}^{n+1}_{+}\right\\}.$

While the original problem is a maximization problem, we add a minus sign to
the objective function to transform it into a minimization problem. Note that
(SQK) is an SDP problem with the matrix variable $Y$ of size
$(n+1)\times(n+1).$ The SDP constraint $Y\in\mathbb{S}^{n+1}_{+}$ is
equivalent to $X-xx^{\top}\in\mathbb{S}^{n}_{+}$ by the Schur complement
lemma. For notational convenience, we will sometimes use $(1,x^{\top};x,X)$ to
denote the matrix $\begin{pmatrix}1&x^{\top}\\\ x&X\end{pmatrix}.$

###  1.2 Our contributions.

In this paper, instead of solving (SQK), we consider the following problem
that transfers the inequality constraint in (SQK) into an affine constraint:

${\rm(SQKE)}\ \min\left\\{-\left\langle C,X\right\rangle:\ a^{\top}Xa-\tau
a^{\top}x=0,\ {\rm diag}\left(X\right)=x,\ Y:=\begin{pmatrix}1&x^{\top}\\\
x&X\end{pmatrix}\in\mathbb{S}^{n+1}_{+}\right\\}.$

Although in general, (SQKE) is not equivalent to (SQK), we can easily
construct an optimal solution of (SQK) from the solution of (SQKE) in section
2. The feasible region of (SQKE) is nonempty because it has a feasible
solution $Y:=(1,0_{1\times n}\,;\,0_{n\times 1},0_{n\times n})$. Note that
from the last two constraints in (SQKE), we have that for any $i\in[n],$
$(1,x_{i}\,;\,x_{i},x_{i})\succeq 0.$ This implies that $0\leq x_{i}\leq 1,$
which in addition implies that every entry of $X$ is inside the interval
$[-1,1].$ Therefore, the feasible region of (SQKE) is compact. (SQKE) is a
typical linear SDP problem, which can be solved with guaranteed convergence by
various well developed solvers such as SDPT3 [35, 36], MOSEK [28], and SDPNAL
[40, 43]. Those solvers are efficient and accurate for problems of moderate
size. However, if $n$ is large, (SQKE) will be challenging to solve by the
aforementioned solvers due to excessive computing cost and memory demand. In
order to solve (SQKE) efficiently, we have conducted both theoretical and
computational studies, which are listed as follows.

####  1.2.1 A modified BM factorization.

In order to overcome the large dimensionality, we consider the following low
rank version of (SQKE),

${\rm(SQKELR)}\ \min\left\\{-\big{\langle}C,\,RR^{\top}\big{\rangle}:\
\|a^{\top}R\|^{2}-\tau a^{\top}Re_{1}=0,\ {\rm
diag}\big{(}RR^{\top}\big{)}=Re_{1},\ R\in\mathbb{R}^{n\times r}\right\\},$

where $e_{1}\in\mathbb{R}^{r}$ is such that its first entry is 1 and all other
entries are zero. The above factorization model comes from the fact that for
any $Y:=(1,x^{\top}\,;\,x,X)\in\mathbb{S}^{n+1}_{+}$ of rank $r$, it can be
factorized as
$\begin{pmatrix}e_{1}&R^{\top}\end{pmatrix}^{\top}\begin{pmatrix}e_{1}&R^{\top}\end{pmatrix}$
for some $R\in\mathbb{R}^{n\times r}.$ Note that our factorization is
different from the classical Burer-Monteiro factorization (see [12, 13])
because we fix the first row of the low rank matrix to be $e_{1}^{\top}.$ By
doing so, we can avoid dealing with the constraint $Y_{11}=1$ because it is
automatically satisfied. Another advantage of this factorization is that we
can prove that there are only finitely many non-regular points in the feasible
region of (SQKELR). However, this is not the case for the usual Burer-Monteiro
factorization. The reason is that for any matrix $R\in\mathbb{R}^{n\times r}$
and orthogonal matrix $Q\in\mathbb{R}^{r\times r},$
$RR^{\top}=RQQ^{\top}R^{\top}.$ This implies that $R$ is equivalent to $RQ$
and one non-regular point will implies infinite many non-regular points.

####  1.2.2 An algebraic variety with good geometric properties.

There are many algorithms and approaches for solving the low rank formulation
of linear SDP problems. A general solver is SDPLR by Burer and Monteiro (see
[12, 13]), which use an augmented Lagrangian method to solve the factorized
SDP problems. Another approach is based on a feasible method, which strictly
preserves the constraints. A feasible method has the advantage that it can
terminate in advance when the primal feasibility is more important than the
optimality of a solution. A famous application of feasible methods is the low
rank SDP relaxation of max-cut problem (see [9, 11, 12, 38]). Because a
feasible method usually relies on the manifold structure of the feasible
region, in section 3, we will analyze the geometric properties of the feasible
region of (SQKELR) i.e.,

$\mathcal{K}_{n,r}:=\left\\{R\in\mathbb{R}^{n\times r}:\
\|a^{\top}R\|^{2}-\tau a^{\top}Re_{1}=0,\ {\rm
diag}\big{(}RR^{\top}\big{)}=Re_{1}\right\\}.$ (1)

We will show that for a generic data
$\left(a,\tau\right)\in\mathbb{R}^{n}\times\mathbb{R}$ that satisfies
Assumption 1, $\mathcal{K}_{n,r}$ is smooth everywhere except for the trivial
point $R=0$. We will also study the local geometric properties of
$\mathcal{K}_{n,r}$ at non-regular points, which is useful for our theoretical
analysis and algorithmic design.

####  1.2.3 A new rank bound for (SQKE) to be equivalent to (SQKELR).

An interesting problem to investigate is when is (SQKE) equivalent to
(SQKELR). Note that there are only $m:=n+2$ affine constraints and one of the
affine constraints of (SQKE) is $Y_{11}=1.$ The compactness of the nonempty
feasible region of (SQKE) implies that it has an optimal solution of rank
$\leq\lceil\sqrt{2(n+2)}\rceil$ (see, e.g., [2, 18, 29]). Thus, when
$r\geq\lceil\sqrt{2(n+2)}\rceil,$ (SQKELR) and (SQKE) have the same optimal
function value. Note that $\lceil\sqrt{2m}\rceil$ is the theoretical rank-
bound, which only relies on the number of constraints of an SDP problem while
ignoring any structural information of the coefficient data. In subsection
4.1, we will derive a new rank-bound that is related to the structure of the
profit matrix $C.$ This new bound is much stronger than
$\lceil\sqrt{2(n+2)}\rceil$ in the case of a binary linear knapsack problem.

####  1.2.4 (SQKELR) doesn’t have spurious local minima.

Another interesting problem to investigate is when does (SQKELR) have no
spurious local minima. Boumal et al. did a series of works (see [4, 9, 10]) to
study when does the Burer-Monteiro formulation of a linear SDP problem have no
spurious local minima. Their analysis is based on the assumption that the
feasible region of the low rank SDP problem is a smooth manifold. In
subsection 4.2, we will show that any second order stationary point of
(SQKELR) is a global optimal solution under a suitable rank condition that is
similar to the results in [4, 9, 10]. Compared with their works, our result
has the advantage that we don’t need any regularity assumption.

####  1.2.5 A Riemannian optimization approach that can handle non-regular
points.

In section 5, we will discuss the details of our algorithmic implementation.
We will use the theoretical results in previous sections to design an
algorithm that combines a Riemannian optimization method with a verification
and escaping strategy at non-regular points. Our algorithm is similar to
algorithm 1 in [34] and both of them can stop after finite number of outer
iterations. We also apply rounding to the output of our algorithm to get an
integer solution that is feasible for (QKP). In section 6, we conduct
numerical experiments to verify the efficiency of our algorithm compared with
the solver SDPLR and another recent solver SketchyCGAL developed by Yurtsever
et al. in [41]. Note that all these three algorithms make use of the low rank
property of (SQKE) but our algorithm is the only feasible method that exploit
the geometric properties of (SQKELR). Our algorithm can be more than 100 times
faster than the other two solvers for some high dimensional problems with
$n\geq 5000$. For a large sparse profit matrix $C$ of size $n=10^{6}$, our
algorithm can find an accurate optimal solution of (SQKE) within half an hour.
We also test the behaviour of our algorithm with a rounding procedure for
solving (QKP). We show that our algorithm can find a nearly optimal solution
for (QKP) very efficiently.

###  1.3 Notations.

In this paper, we often omit stating the dimension of a vector or matrix if it
is already clear from the context. Without mentioning explicitly, $\|.\|$
always denote the matrix Frobenius norm. We denote the trace of a square
matrix $X$ as $\mathrm{Tr}(X)$. For two matrices with the same size,
$\left\langle A,B\right\rangle:=\mathrm{Tr}\left(AB^{\top}\right)$ is the
matrix inner product; $e$ is the vectors of all ones; $e_{1}$ is the vector
such that its first entry is 1 and all the other entries are zero; ${\rm
S}^{r-1}:=\left\\{x\in\mathbb{R}^{r}:\ \|x\|=1\right\\};$ $\mathbb{S}^{n}$ and
$\mathbb{S}^{n}_{+}$ denote the space of $n\times n$ symmetric matrices and
its subset of symmetric positive semidefinite matrices, respectively. For
notational simplicity, we use $X\succeq 0$ to denote $X\in\mathbb{S}^{n}_{+}$.
For $x\in\mathbb{R}^{n},$ ${\rm diag}\left(x\right)\in\mathbb{S}^{n}$ is the
diagonal matrix whose diagonal is given by $x.$ For $X\in\mathbb{R}^{n\times
n}$, ${\rm diag}(X)\in\mathbb{R}^{n}$ is the diagonal vector of $X$, and ${\rm
Diag}(X)={\rm diag}({\rm diag}(X)).$ For $A,B\in\mathbb{R}^{m\times n}$,
$A\circ B\in\mathbb{R}^{m\times n}$ is defined by $\left(A\circ
B\right)_{ij}=a_{ij}b_{ij}.$

###  1.4 Organization of the rest of the paper.

In section 2, we discuss when is (SQK) equivalent to (SQKE) and how to solve
(SQK) via solving (SQKE). In section 3, we study the geometric properties of
the algebraic variety $\mathcal{K}_{n,r}.$ In section 4, we discuss the
equivalence between (SQKE) and (SQKELR), which includes a new rank-bound for
(SQKELR). In section 5, we design an algorithm to solve (SQKELR). In section
6, we conduct numerical experiments to compare our algorithm with SDPLR and
SketchyCGAL. In section 7, we give a brief conclusion.

##  2 Relation between (SQK) and (SQKE).

In this section, we consider how to construct an optimal solution of (SQK) by
solving (SQKE). We have the following proposition.

######  Proposition 2.1.

Suppose $(a,\tau)$ satisfy Assumption 1. If every entry of $C$ is nonnegative,
then (SQK) and (SQKE) have the same optimal function value.

###### Proof.

Proof. Note that Assumption 1 implies that $\tau>0.$ We only have to prove
that there exists an optimal solution $Y:=(1,x^{\top}\,;\,x,X)$ of (SQK) such
that $a^{\top}Xa-\tau a^{\top}x=0.$ In order to prove this, assume that
$a^{\top}Xa-\tau a^{\top}x<0$ for some optimal solution $Y$ of (SQK) because
otherwise there is nothing to prove. Consider another matrix
$\widehat{Y}:=(1,\hat{x}^{\top}\,;\,\hat{x},\widehat{X})=ee^{\top}.$ From
Assumption 1, we have that

$a^{\top}\widehat{X}a-\tau a^{\top}\hat{x}=\big{(}a^{\top}e\big{)}^{2}-\tau
a^{\top}e>0.$ (2)

Also, because every entry of $C$ is nonnegative, and every entry of $Y$ is
inside the interval $[-1,1]$ as shown in the introduction. We have that

$-\big{\langle}C,\,\widehat{X}\big{\rangle}\leq-\big{\langle}C,\,X\big{\rangle}.$
(3)

Now, consider $Y_{t}:=t\widehat{Y}+(1-t)Y.$ From (2), (3) and the convexity of
(SQK), we have that there exists $t\in(0,1)$ such that $Y_{t}$ is inside the
feasible region of (SQKE) with its function value less than or equal to the
function value at $Y.$ From the optimality of $Y,$ we have that $Y_{t}$ is
also an optimal solution of (SQK). Thus, we have found an optimal solution of
(SQK), which is also an optimal solution of (SQKE).

Proposition 2.1 requires that the profit matrix $C$ is nonnegative. However,
the following proposition tells us that if $C$ is not nonnegative, we can
still construct an optimal solution of (SQK) by solving (SQKE) and the
following simpler problem.

${\rm(SQKS)}\ \min\left\\{-\left\langle C,X\right\rangle:\ {\rm
diag}\left(X\right)=x,\ \begin{pmatrix}1&x^{\top}\\\
x&X\end{pmatrix}\in\mathbb{S}^{n+1}_{+}\right\\}.$

######  Proposition 2.2.

Suppose $(a,\tau)$ satisfy Assumption 1. If (SQK) and (SQKE) don’t have the
same optimal function value, then any optimal solution of (SQKS) is also an
optimal solution of (SQK).

###### Proof.

Proof. Let $V,V_{E},V_{S}$ be the optimal function values of (SQK), (SQKE) and
(SQKS) respectively. It is easy to see that $V_{E}\geq V\geq V_{S}.$ Because
(SQK) and (SQKE) have different optimal function values, we have that
$V_{E}>V\geq V_{S}.$ Now, assume by contradiction that there exists an optimal
solution $Y_{S}:=(1,x^{\top}\,;\,x,X)$ of (SQKS) such that $a^{\top}Xa-\tau
a^{\top}x>0.$ Because $V_{E}>V,$ we have that there exists an optimal solution
$Y:=(1,\hat{x}^{\top}\,;\,\hat{x},\widehat{X})$ of (SQK) such that
$a^{\top}\widehat{X}a-\tau a^{\top}\hat{x}<0.$ Now, define
$Y_{E}^{t}=tY+(1-t)Y_{S}.$ We have that there exists $t\in(0,1)$ such that
$Y_{E}^{t}$ is feasible for (SQKE). Also, from $V_{E}>V\geq V_{S},$ we have
that the function value of $Y_{E}^{t}$ is strictly smaller than $V_{E}.$ This
contradicts to the fact that $V_{E}$ is the optimal function value of (SQKE).

In order to solve (SQK), we can first solve (SQKS) to get an optimal solution
$Y.$ If $Y$ is feasible for (SQK), then $Y$ is also an optimal solution of
(SQK). Otherwise if $Y$ is not feasible for (SQK), then by Proposition 2.2,
(SQK) and (SQKE) must have the same optimal function value. In this case, we
only have to solve (SQKE) to get an optimal solution of (SQK). Note that
(SQKS) is easy to solve because its low rank factorization only has the
nonlinear constraint ${\rm diag}\left(RR^{\top}\right)=Re_{1},$ which is
equivalent to that ${\rm
diag}\left(\left(2R-ee_{1}^{\top}\right)\left(2R-ee_{1}^{\top}\right)^{\top}\right)=e.$
Thus, solving the low rank version of (SQKS) is equivalent to minimizing a
quadratic function on the following oblique manifold

${\rm OB}_{n,r}:=\left\\{R\in\mathbb{R}^{n\times r}:\ {\rm
diag}\left(RR^{\top}\right)=e\right\\}.$

This problem can be solved efficiently by a widely used toolbox called manopt
(see [8]). Since (SQK) can always be reduced to (SQKE), in the sequel, we will
focus on analysing the property of (SQKE) instead of (SQK).

##  3 Geometric properties of $\mathcal{K}_{n,r}$.

In this section, we focus on the geometric properties of the feasible region
of (SQKELR), $\mathcal{K}_{n,r}$ defined in (1). Note that $\mathcal{K}_{n,r}$
is not empty because $0_{n\times r}\in\mathcal{K}_{n,r}.$ We call this point a
trivial point because it is of little use in practice. Actually, this point
corresponds to the case that we don’t select any item in (QKP).

###  3.1 Regularity of $\mathcal{K}_{n,r}$.

From now on, we say that a point $R\in\mathcal{K}_{n,r}$ is regular if the
linear independent constraint qualification (LICQ) holds at $R.$ Otherwise we
say that $R$ is non-regular. Note that LICQ is a sufficiently condition for a
point to be smooth on an algebraic variety. If $R\in\mathcal{K}_{n,r}$ is
regular, then $\mathcal{K}_{n,r}$ is a smooth manifold locally around $R$ and
its tangent cone is a linear space, which is also called the tangent space. It
is easy to check that for any $(a,\tau)$ satisfying Assumption 1, $0_{n\times
r}$ is a non-regular point. The following proposition characterizes all the
non-regular points on $\mathcal{K}_{n,r}.$

######  Proposition 3.1.

Suppose $(a,\tau)$ satisfy Assumption 1. The non-regular points of
$\mathcal{K}_{n,r}$ are given by

$\mathcal{N}_{n,r}:=\\{0_{n\times r}\\}\sqcup\left\\{ve_{1}^{\top}:\
v\in\\{0,1\\}^{n},\ a^{\top}v=\tau\right\\}.$ (4)

###### Proof.

Proof. Since the trivial point $0_{n\times r}$ is non-regular, we only have to
characterize the set of non-trivial and non-regular points of
$\mathcal{K}_{n,r}.$ Assume that $R\in\mathcal{K}_{n,r}$ is non-trivial and
non-regular. Because the LICQ doesn’t hold at $R,$ we have that there exists
nonzero $\left(\mu,\lambda\right)\in\mathbb{R}^{n}\times\mathbb{R}$ such that

${\rm diag}(\mu)\left(2R-ee_{1}^{\top}\right)+\lambda a\left(2a^{\top}R-\tau
e_{1}^{\top}\right)=0.$ (5)

Because $R\in\mathcal{K}_{n,r},$ the constraint ${\rm
diag}\left(RR^{\top}\right)=Re_{1}$ implies that $2R-ee_{1}^{\top}\in{\rm
OB}_{n,r}.$ That is, every row of $2R-ee_{1}^{\top}$ is on the unit sphere of
$\mathbb{R}^{1\times r}.$ If there exists $i\in[n]$ such that $\mu_{i}=0,$
then from (5) and the fact that every entry of $a$ is positive, we have that
$\lambda\left(2a^{\top}R-\tau e_{1}^{\top}\right)=0.$ This implies that ${\rm
diag}(\mu)\left(2R-ee_{1}^{\top}\right)=0$ and so $\mu=0.$ Because
$\left(\mu,\lambda\right)$ is nonzero, we have that $\lambda\neq 0.$ Thus,
$\left(2a^{\top}R-\tau e_{1}^{\top}\right)=0$ and so
$a^{\top}R=\frac{1}{2}\tau e_{1}^{\top}.$ It follows that
$\|a^{\top}R\|^{2}-\tau a^{\top}Re_{1}=-\frac{\tau^{2}}{4}<0,$ which
contradicts to the fact that $R\in\mathcal{K}_{n,r}.$ Thus we conclude the
fact that for any $i\in[n],$ $\mu_{i}\neq 0.$

Now from (5), we have that

$2R-ee_{1}^{\top}+{\rm diag}\left(\mu\right)^{-1}\lambda
a\big{(}2a^{\top}R-\tau e_{1}^{\top}\big{)}=0.$ (6)

Because $2R-ee_{1}^{\top}\in{\rm OB}_{n,r}$ and ${\rm
diag}\left(\mu\right)^{-1}\lambda a\left(2a^{\top}R-\tau e_{1}^{\top}\right)$
is rank 1, we have that $2R-ee_{1}^{\top}=ub^{\top}$ for some
$u\in\\{-1,1\\}^{n}$ and $b\in{\rm S}^{r-1}$ such that the first nonzero entry
of $b$ is positive. From (6), we have that

$b^{\top}=\alpha\left(a^{\top}ub^{\top}+(a^{\top}e-\tau)e_{1}^{\top}\right),$
(7)

for some nonzero normalization factor $\alpha\in\mathbb{R}.$ If $b$ is not
parallel to $e_{1}$, then from (7), we have that $a^{\top}e-\tau=0,$ which
contradicts the fact that $a^{\top}e>\tau.$ Thus $b$ is parallel to $e_{1}$,
and hence $b=e_{1},$ since $b\in{\rm S}^{r-1}$ and its first nonzero entry is
positive. Therefore, we have that $R=\left(\frac{e+u}{2}\right)e_{1}^{\top}.$
Define $v:=\frac{e+u}{2},$ we have that $v\in\\{0,1\\}^{n}$ and
$R=ve_{1}^{\top}.$ Because $R\in\mathcal{K}_{n,r},$ we have that

$0=\|a^{\top}R\|^{2}-\tau a^{\top}Re_{1}=\big{(}a^{\top}v\big{)}^{2}-\tau
a^{\top}v.$ (8)

Because $R$ is non-trivial and every entry of $a$ is positive, we have that
$a^{\top}v>0.$ This together with (8) implies that $a^{\top}v=\tau.$
Therefore, $R\in\mathcal{N}_{n,r}.$

For any nontrivial point $R=ve_{1}^{\top}\in\mathcal{N}_{n,r}.$ Because
$v\in\\{0,1\\}^{n}$ and $a^{\top}v=\tau,$ we have that
$R=ve_{1}^{\top}\in\mathcal{K}_{n,r}.$ Thus, we only have to check that
$ve_{1}^{\top}$ is a non-regular point of $\mathcal{K}_{n,r}.$ This can be
easily verified by choosing $\left(\mu,\lambda\right)=\left({\rm
diag}(a)\left(2v-e\right),\frac{-1}{\tau}\right)$ in (5).

From Proposition 3.1, one can see that the non-trivial non-regular points of
$\mathcal{K}_{n,r}$ corresponds to the integer feasible solutions of (QKP)
such that the knapsack upper bound are attained. Proposition 3.1 has the
following corollary, which says that $\mathcal{K}_{n,r}$ is likely to be
smooth everywhere except for the trivial point.

######  Corollary 3.2.

For a generic $(a,\tau)\in\mathbb{R}^{n}\times\mathbb{R}$ that satisfies
Assumption 1, $\mathcal{N}_{n,r}=\\{0_{n\times r}\\}$.

###### Proof.

Proof. First, note that the set of vectors $(a,\tau)$ satisfying Assumption 1
is an open set in $\mathbb{R}^{n}\times\mathbb{R}$ with positive measure. From
(4), it is easy to see that the set of $(a,\tau)$’s that satisfies Assumption
1 and $\mathcal{N}_{n,r}\neq\\{0_{n\times r}\\}$ is the following set

$\mathcal{W}_{n,\tau}:=\cup_{v\in\\{0,1\\}^{n}}\left\\{(a,\tau)\in\mathbb{R}^{n}\times\mathbb{R}:\
a^{\top}v=\tau,\forall\ i\in[n]\ 0<a_{i}<\tau,\ e^{\top}a>\tau\right\\}.$ (9)

Note that the affine constraint $a^{\top}v=\tau$ implies that
$\mathcal{W}_{n,\tau}$ is the union of finitely many sets of measure zero in
$\mathbb{R}^{n}\times\mathbb{R}$. Therefore, we have that for a generic
$(a,\tau)\in\mathbb{R}^{n}\times\mathbb{R}$ that satisfies Assumption 1,
$\mathcal{N}_{n,r}=\\{0_{n\times r}\\}.$

###  3.2 Local geometric properties at non-regular points of
$\mathcal{K}_{n,r}$.

Although Corollary 3.2 implies that there are no non-trivial non-regular
points of $\mathcal{K}_{n,r}$ with probability one with respect to the data
$\left(a,\tau\right)$ satisfying Assumption 1, for theoretical completeness,
we will discuss the local geometric properties at non-regular points of
$\mathcal{K}_{n,r}$ in this subsection. Before presenting our result, we need
the definition of the tangent cone and second order tangent set. For any
closed set $S\subset\mathbb{R}^{n}$ and $x\in S,$ the inner tangent cone (see
[6] section 2) $\mathcal{T}^{i}_{S}(x)$ is defined as follows:

$\mathcal{T}^{i}_{S}(x):=\left\\{h\in\mathbb{R}^{n}:\ {\rm
dist}\left(x+th,S\right)=o(t),\ t\geq 0\right\\},$ (10)

where ${\rm dist}(x,S):=\inf_{y\in S}\|x-y\|.$ For any vector
$h\in\mathcal{T}^{i}_{S}(x),$ the inner second order tangent set is defined as
follows (see [6] definition 3.28):

$\mathcal{T}^{i,2}_{S}(x,h):=\left\\{w\in\mathbb{R}^{n}:\ {{\rm
dist}}\Big{(}x+th+\frac{1}{2}t^{2}w,S\Big{)}=o(t^{2}),t\geq 0\right\\}.$ (11)

Note that there are also other kinds of tangent cones and second order tangent
sets (see section 2 and section 3 in [6] for more details). However, in this
paper, we only use the definition of inner tangent cone and inner second order
tangent set. Thus, for convenience, we use $\mathcal{T}_{S}(x)$ and
$\mathcal{T}^{2}_{S}(x,h)$ to denote $\mathcal{T}^{i}_{S}(x)$ and
$\mathcal{T}^{i,2}_{S}(x,h)$ respectively.

######  Proposition 3.3.

Suppose $(a,\tau)$ satisfies Assumption 1 and $r\geq 2.$ Let
$R=ve_{1}^{\top}\in\mathcal{N}_{n,r},$ where $v\in\\{0,1\\}^{n}.$ Then the
tangent cone of $\mathcal{K}_{n,r}$ at $R$ satisfies

$\widetilde{\mathcal{T}}:=\left\\{[0,H]:\left\langle
aa^{\top}-\sigma_{v}\tau\cdot{\rm diag}\left(a\circ
d\right),HH^{\top}\right\rangle=0,\
H\in\mathbb{R}^{n\times(r-1)}\right\\}\subset\mathcal{T}_{\mathcal{K}_{n,r}}\left(R\right),$
(12)

where $d=2v-e$ and $\sigma_{v}:=1$ if $v\neq 0,$ and $\sigma_{v}:=-1$ if
$v=0.$ In addition, for any $h=[0,H]\in\widetilde{\mathcal{T}},$ $w=[-2{\rm
diag}\left(HH^{\top}\right)\circ
d,\,0]\in\mathcal{T}^{2}_{\mathcal{K}_{n,r}}\left(R,h\right).$

###### Proof.

Proof. For any $[0,H]\in\widetilde{\mathcal{T}},$ define the matrix valued
function
$R(t,W):(-\delta,\delta)\times\mathbb{R}^{n\times(r-1)}\rightarrow\mathbb{R}^{n\times
r}$ such that

$R(t,W):=\left[\left(v-\frac{e}{2}\right)\circ\sqrt{e-4{\rm
diag}\big{(}\left(tH+t^{3}W\right)\left(tH+t^{3}W\right)^{\top}\big{)}}+\frac{e}{2},\
tH+t^{3}W\right],$ (13)

where the square-root is taken componentwise and $\delta>0$ is sufficiently
small for the square root to be well-defined. Note that the first column of
$R(t,W)$ is derived from solving the equation ${\rm
diag}\left(RR^{\top}\right)-Re_{1}=0.$ It is easy to check that ${\rm
diag}\left(R(t,W)R(t,W)^{\top}\right)-R(t,W)e_{1}=0$. By using the Taylor
expansion of the square root term of $R(t,W),$ we get

$\displaystyle R(t,W)=$ $\displaystyle\Big{[}v-\left(t^{2}{\rm
diag}\big{(}HH^{\top}\big{)}+2t^{4}{\rm diag}\big{(}HW^{\top}\big{)}+t^{4}{\rm
diag}^{2}\big{(}HH^{\top}\big{)}\right)\circ d+t^{6}g(t,W),$
$\displaystyle\;\;tH+t^{3}W\Big{]},$ (14)

where $g(t,W)$ is a smooth function. In order to get
$R(t,W)\in\mathcal{K}_{n,r},$ we only have to guarantee that

$4\|a^{\top}R(t,W)\|^{2}-4\tau a^{\top}R(t,W)e_{1}=\|2a^{\top}R(t,W)-\tau
e_{1}^{\top}\|^{2}-\tau^{2}=0.$ (15)

Note that from the definition of $\mathcal{N}_{n,r}$ and $\sigma_{v}$, we have
that $2a^{\top}v-\tau=\sigma_{v}\tau.$ By using (3.2), we have that

$\displaystyle\|2a^{\top}R(t,W)-\tau e_{1}^{\top}\|^{2}-\tau^{2}$
$\displaystyle=-4t^{2}\sigma_{v}\tau\cdot a^{\top}\big{(}{\rm
diag}\big{(}HH^{\top}\big{)}\circ d\big{)}+4t^{4}\left(a^{\top}\big{(}{\rm
diag}\big{(}HH^{\top}\big{)}\circ
d\big{)}\right)^{2}-8t^{4}\sigma_{v}\tau\cdot a^{\top}\big{(}{\rm
diag}\big{(}HW^{\top}\big{)}\circ d\big{)}$
$\displaystyle\quad-4t^{4}\sigma_{v}\tau\cdot a^{\top}\big{(}{\rm
diag}^{2}\big{(}HH^{\top}\big{)}\circ
d\big{)}+8t^{4}\big{\langle}a^{\top}H,\,a^{\top}W\big{\rangle}+4t^{2}\|a^{\top}H\|^{2}+t^{6}k\left(t,W\right)$
$\displaystyle=8t^{4}\left\langle aa^{\top}H-\sigma_{v}\tau\cdot{\rm
diag}\left(a\circ d\right)H,W\right\rangle$
$\displaystyle\quad+4t^{4}\left(\big{(}a^{\top}\big{(}{\rm
diag}\big{(}HH^{\top}\big{)}\circ d\big{)}\big{)}^{2}-\sigma_{v}\tau\cdot
a^{\top}\big{(}{\rm diag}^{2}\big{(}HH^{\top}\big{)}\circ
d\big{)}\right)+t^{6}k(t,W),$ (16)

where the second equality comes from that $\left\langle
aa^{\top}-\sigma_{v}\tau\cdot{\rm diag}\left(a\circ
d\right),HH^{\top}\right\rangle=0.$ Also, $k(t,W)$ is a smooth function.

######  Claim 3.4.

If $H\neq 0,$ then $aa^{\top}H-\sigma_{v}\tau\cdot{\rm diag}\left(a\circ
d\right)H\neq 0.$

###### Proof.

Proof of Claim 3.4. Assume by contradiction that
$aa^{\top}H-\sigma_{v}\tau\cdot{\rm diag}\left(a\circ d\right)H=0.$ Because
every entry of $a$ is positive, we have that
$H=\frac{da^{\top}H}{\tau\sigma_{v}}.$ If $a^{\top}H=0,$ then
$H=\frac{da^{\top}H}{\tau\sigma_{v}}=0,$ which contradicts to the fact that
$H\neq 0.$ Thus, $a^{\top}H\neq 0.$ Because
$H=\frac{da^{\top}H}{\tau\sigma_{v}},$ we have that
$a^{\top}H=a^{\top}\frac{da^{\top}H}{\tau\sigma_{v}},$ which implies that
$a^{\top}d=\tau\sigma_{v}.$ If $v\neq 0,$ then $\sigma_{v}=1$, and from (4),
$a^{\top}v=\tau.$ This together with $a^{\top}d=\tau$ and $d=2v-e$ implies
that $a^{\top}e=\tau,$ which contradicts to the fact that $a^{\top}e>\tau.$ If
$v=0,$ then $\sigma_{v}=-1$ and $d=2v-e=-e.$ It follows that
$a^{\top}e=-a^{\top}d=-\tau\sigma_{v}=\tau.$ Similarly, we get a
contradiction.

Now we define the following function
$F:(-\delta,\delta)\times\mathbb{R}^{n\times(r-1)}\rightarrow\mathbb{R}$ such
that

$\displaystyle F(t,W):=8\left\langle aa^{\top}H-\sigma_{v}\tau\cdot{\rm
diag}\left(a\circ d\right)H,W\right\rangle$
$\displaystyle\qquad\qquad+4\Big{(}\big{(}a^{\top}\big{(}{\rm
diag}(HH^{\top})\circ d\big{)}\big{)}^{2}-\sigma_{v}\tau\cdot
a^{\top}\big{(}{\rm diag}^{2}(HH^{\top})\circ d\big{)}\Big{)}+t^{2}k(t,W).$
(17)

If $H\neq 0,$ then from Claim 3.4, there exists
$W_{0}\in\mathbb{R}^{n\times(r-1)}$ such that $F(0,W_{0})=0.$ Moreover, since
$\nabla_{W}F(0,W_{0})=8\left(aa^{\top}H-\sigma_{v}\tau\cdot{\rm
diag}\left(a\circ d\right)H\right)\neq 0.$ From Lemma A.1, we have that there
exists a continuously differentiable mapping
$\widehat{W}:(-\epsilon,\epsilon)\rightarrow\mathbb{R}^{n\times(r-1)}$ such
that $\widehat{W}(0)=W_{0}$ and $F(t,\widehat{W}(t))=0.$ This implies that
$R(t,\widehat{W}(t))\in\mathcal{K}_{n,r}$ for any $t\in(-\epsilon,\epsilon).$
Because ${\rm dist}\big{(}R(t,\widehat{W}(t)),\left[v-t^{2}{\rm
diag}\left(HH^{\top}\right)\circ d,tH\right]\big{)}=O(t^{3}),$ we have that
$h=[0,H]\in\mathcal{T}_{\mathcal{K}_{n,r}}\left(R\right)$ and $w=[-2{\rm
diag}\left(HH^{\top}\right)\circ
d,0]\in\mathcal{T}_{\mathcal{K}_{n,r}}^{2}\left(R,h\right).$ If $H=0,$ then it
is obvious that $0=[0,0]\in\mathcal{T}_{\mathcal{K}_{n,r}}\left(R\right)$ and
$w=[-2{\rm diag}\left(HH^{\top}\right)\circ
d,0]=0\in\mathcal{T}_{\mathcal{K}_{n,r}}^{2}\left(R,0\right).$

######  Remark 3.5.

In order to save space, we don’t provide the complete local geometric
properties at non-regular points of $\mathcal{K}_{n,r}$ in Proposition 3.3.
The converse of (12) is also correct and it can be proved similarly based on
the Taylor expansion used the in the proof of Proposition 3.3. However, the
above properties are already enough for our later theoretical analysis and
algorithmic design.

##  4 Equivalence between (SQKE) and (SQKELR).

In this section, we will study the equivalence between (SQKE) and (SQKELR). We
will first derive a new rank condition under which (SQKE) and (SQKELR) have
the same global optimal value. After that, we will discuss when a second order
stationary point of (SQKLRE) is also a global optimal solution.

###  4.1 A new rank-bound for (SQKELR).

In this subsection, we discuss the equivalence between (SQKE) and (SQKELR). As
we mentioned in the introduction, (SQKE) and (SQKELR) are equivalent when
$r\geq\sqrt{2(n+2)}.$ Here, we will provide a new rank-bound which relies on
the structure of the matrix $C.$ Before showing our main result, we present
some background knowledge on the minimum rank problem of a graph. Consider a
simple undirected graph $G$ with vertices set $V(G)=[n]$ and edges set
$E(G)\subset\binom{[n]}{2}.$ We define the following set

$\mathcal{Q}(G):=\left\\{A\in\mathbb{S}^{n}:\ \forall i\neq j,\ A_{ij}\neq
0\Leftrightarrow ij\in E(G)\right\\},$ (18)

which is the set of symmetric matrices whose sparsity patterns are specified
by the adjacency matrix of $G.$ The $minimum\ rank$ of $G$ (see [24]) is
defined as follows:

${\rm mr}(G):=\min\left\\{{\rm rank}\left(A\right):\
A\in\mathcal{Q}(G)\right\\}.$ (19)

For a given graph $G$, determining ${\rm mr}(G)$ is a difficult problem that
has received considerable attention (see [15] for a survey). In the next
theorem, we give a new rank-bound of (SQKELR) using ${\rm mr}(G).$

######  Theorem 4.1.

Consider the problem (SQKE), and $(a,\tau)$ satisfies Assumption 1. Suppose
$C\in\mathcal{Q}(G)$ for some graph $G$ and there exists another graph $G_{1}$
such that $V(G_{1})=[n]$ and $E(G)\subset E(G_{1}).$ Then there is an optimal
solution $Y^{*}=(1,{x^{*}}^{\top}\,;\,x^{*},X^{*})$ of (SQKE) such that ${\rm
rank}\left(Y^{*}\right)\leq n+2-{\rm mr}(G_{1}).$

###### Proof.

Proof. Because $C\in\mathcal{Q}(G)$ and $E(G)\subset E(G_{1}),$ there exists a
matrix sequence $\\{C_{k}\\}\subset\mathbb{S}^{n}$ such that for any
$k\in\mathbb{N}^{+}$, $C_{k}\in\mathcal{Q}(G_{1})$ and
$\lim_{k\rightarrow\infty}\|C_{k}-C\|=0.$ Let (SQKE)k be the optimization
problem obtained by replacing $C$ with $C_{k}$ in (SQKE). Suppose
$Y_{k}=(1,{x_{k}^{*}}^{\top}\,;\,x_{k}^{*},X_{k}^{*})$ is an optimal solution
of problem (SQKE)${}_{k}.$ From the compactness of the feasible region of
(SQKE), and $\lim_{k\rightarrow\infty}\|C_{k}-C\|=0,$ we have that any
accumulation point of $\\{Y_{k}\\}$ is an optimal solution of (SQKE). Without
loss of generality, we may assume that $Y_{i_{k}}\rightarrow Y^{*}$ for some
convergent sequence $\\{Y_{i_{k}}\\}$. Thus,
$Y^{*}=(1,{x^{*}}^{\top}\,;\,x^{*},X^{*})$ is an optimal solution of (SQKE).

######  Claim 4.2.

For any $k\in\mathbb{N}^{+},$ ${\rm rank}\left(Y_{i_{k}}\right)\leq n+2-{\rm
mr}(G_{1}).$

###### Proof.

Proof of Claim 4.2. Note that ${\rm mr}(G_{1})\leq|G_{1}|-1=n-1$ (see
observation 1.2 of [15]). Assume by contradiction that there exists
$k\in\mathbb{N}^{+}$ such that ${\rm rank}\left(Y_{i_{k}}\right)>n+2-{\rm
mr}(G_{1}).$ Then we have that ${\rm rank}\left(Y_{i_{k}}\right)\geq 4.$ Let
(SQKELR)${}_{i_{k}}$ be the low rank formulation of (SQKE)${}_{i_{k}}$ with
$r={\rm rank}\left(Y_{i_{k}}\right).$ Because $Y_{i_{k}}$ the an optimal
solution of (SQKE)${}_{i_{k}},$ we have that (SQKELR)${}_{i_{k}}$ is
equivalent to (SQKE)${}_{i_{k}}.$ Moreover, consider the low rank
factorization
$Y_{i_{k}}=\begin{pmatrix}e_{1}&R_{i_{k}}^{\top}\end{pmatrix}^{\top}\begin{pmatrix}e_{1}&R_{i_{k}}^{\top}\end{pmatrix}$
in (SQKELR)${}_{i_{k}}.$ Then $R_{i_{k}}$ is an optimal solution of
(SQKELR)${}_{i_{k}}$ such that ${\rm rank}\left(R_{i_{k}}\right)\geq{\rm
rank}\left(Y_{i_{k}}\right)-1\geq 3.$ Because ${\rm
rank}\left(R_{i_{k}}\right)\geq 3$, by Proposition 3.1, it satisfies the LICQ
of (SQKELR)${}_{i_{k}}$. Thus $R_{i_{k}}$ satisfies the following first order
KKT condition of (SQKELR)${}_{i_{k}}$

$-2\left(C_{i_{k}}+{\rm diag}\left(\mu\right)\right)R_{i_{k}}+\mu
e_{1}^{\top}-\lambda a\big{(}2a^{\top}R_{i_{k}}-\tau e_{1}^{\top}\big{)}=0,$
(20)

for some $\left(\mu,\lambda\right)\in\mathbb{R}^{n}\times\mathbb{R}.$
Rearranging (20), we get

$\big{(}\mu,-2\left(C_{i_{k}}+{\rm
diag}\left(\mu\right)\right)\big{)}\begin{pmatrix}e_{1}^{\top}\\\
R_{i_{k}}\end{pmatrix}=\lambda a\big{(}2a^{\top}R_{i_{k}}-\tau
e_{1}^{\top}\big{)}.$ (21)

Because $C_{i_{k}}\in\mathcal{Q}(G_{1}),$ we have that

${\rm rank}\Big{(}\big{(}\mu,-2\left(C_{i_{k}}+{\rm
diag}\left(\mu\right)\right)\big{)}\Big{)}\geq{\rm
rank}\Big{(}-2\left(C_{i_{k}}+{\rm diag}\left(\mu\right)\right)\Big{)}\geq{\rm
mr}\left(G_{1}\right).$ (22)

Applying (22) in (21), we have that

$\displaystyle 1\geq{\rm rank}\left(\big{(}\mu,-2\left(C_{i_{k}}+{\rm
diag}\left(\mu\right)\right)\big{)}\begin{pmatrix}e_{1}^{\top}\\\
R_{i_{k}}\end{pmatrix}\right)$ $\displaystyle\geq{\rm
rank}\Big{(}\big{(}\mu,-2\left(C_{i_{k}}+{\rm
diag}\left(\mu\right)\right)\big{)}\Big{)}+{\rm
rank}\left(\begin{pmatrix}e_{1}^{\top}\\\
R_{i_{k}}\end{pmatrix}\right)-(n+1),$ $\displaystyle\geq{\rm mr}(G_{1})+{\rm
rank}\left(\begin{pmatrix}e_{1}^{\top}\\\ R_{i_{k}}\end{pmatrix}\right)-(n+1)$
(23)

where the second inequality follows from Sylvester’s rank inequality (see 10.5
of [30]). Therefore, ${\rm rank}\left(Y_{i_{k}}\right)={\rm
rank}\big{(}\begin{pmatrix}e_{1}\
R_{i_{k}}^{\top}\end{pmatrix}^{\top}\big{)}\leq n+2-{\rm mr}(G_{1}),$ which
contradicts to the assumption that ${\rm rank}\left(Y_{i_{k}}\right)>n+2-{\rm
mr}(G_{1}).$

From Claim 4.2, since $Y_{i_{k}}\rightarrow Y^{*},$ we have that ${\rm
rank}\left(Y^{*}\right)\leq n+2-{\rm mr}(G_{1}).$

The rank-bound of $n+2-{\rm rank}(G_{1})$ in Theorem 4.1 is different from the
rank-bound of $\sqrt{2(n+2)}$ because it is related to the structure of the
underlying graph $G$ of matrix $C$. Many existing results in the literature
for the minimum rank problem can be applied to estimate the number $n+2-{\rm
rank}(G_{1})$ (see [15, 24, 23]). A simple application of Theorem 4.1 is on
(SQKE) for a binary linear knapsack problem. In this case, $C$ is a diagonal
matrix and its underlying graph $G$ is an empty graph. We have that
$E(G)\subset E(P_{n}),$ where $P_{n}$ is a path with $n$ vertices. Because
${\rm mr}(P_{n})=n-1$ (see observation 1.2 of [15]), the rank-bound becomes
$n+2-{\rm rank}(P_{n})=3,$ which is much better than $\sqrt{2(n+2)}$ when $n$
is large. The SDP relaxation of a linear knapsack problem has been
demonstrated to be tighter than its linear programming relaxation (see lemma
2.8 in [22]). Our results additionally implies that we can use (SQKE) to get a
better bound than linear programming relaxation with almost the same memory
usage. Another interesting case is when the matrix $C$ is a banded matrix such
that $C_{ij}=0$ if $|i-j|>k.$ In that case, Theorem 4.1 implies that (SQKE)
has an optimal solution of rank $\leq 2+k.$

###  4.2 Second order stationarity of a non-regular point of (SQKELR).

In this subsection, we discuss when is a second order stationary point of
(SQKELR) also a global optimal solution. Traditional optimality conditions are
applicable only to regular points in the feasible region of an optimization
problem (see Section 3 of [6]). Thus, when $R$ is a regular point of (SQKELR),
we will call $R$ a second order stationary point if it satisfies the first
order and second order KKT conditions (see theorem 12.1 and 12.5 in [39]).
However, since $\mathcal{K}_{n,r}$ may have non-regular points, we need to
define the second order stationarity for non-regular points of (SQKELR). From
Proposition 3.3 on the local geometric properties of a non-regular point in
$\mathcal{K}_{n,r}$, we have the following lemma.

######  Lemma 4.3.

Suppose $(a,\tau)$ satisfies Assumption 1 and $r\geq 2.$ Let
$R=ve_{1}^{\top}\in\mathcal{N}_{n,r},$ where $v\in\\{0,1\\}^{n}$ and
$d\in\\{-1,1\\}^{n},\sigma_{v}\in\\{-1,1\\}$ are defined as in Proposition
3.3. If $R$ is a local minimum of (SQKELR), then the optimal value of the
following problem is nonnegative.

$\min\left\\{\big{\langle}\left(2{\rm diag}\left(\left(Cv\right)\circ
d\right)-C\right)H,\,H\big{\rangle}:\ \begin{array}[]{l}\left\langle
aa^{\top}-\sigma_{v}\tau\cdot{\rm diag}\left(a\circ
d\right),HH^{\top}\right\rangle=0,\\\\[5.0pt] \|H\|^{2}=1,\
H\in\mathbb{R}^{n\times(r-1)}\end{array}\right\\}.$ (24)

###### Proof.

Proof. Assume by contradiction that there exists a feasible solution $H$ of
(24) such that

$\left\langle\left(2{\rm diag}\left(\left(Cv\right)\circ
d\right)-C\right)H,H\right\rangle<0.$ (25)

From Proposition 3.3, we have that
$h:=[0,H]\in\mathcal{T}_{\mathcal{K}_{n,r}}\left(R\right)$ and
$w:=\left[-2{\rm diag}\left(HH^{\top}\right)\circ
d,0\right]\in\mathcal{T}^{2}_{\mathcal{K}_{n,r}}(R,h).$ Define

$\widehat{R}\left(t,H\right)=\left[v-t^{2}\big{(}{\rm diag}(HH^{\top})\circ
d\big{)},tH\right]\in\mathbb{R}^{n\times r}.$ (26)

We have that ${\rm
dist}\big{(}\widehat{R}\left(t,H\right),\mathcal{K}_{n,r}\big{)}=o(t^{2}).$
Define the function $F:\mathbb{R}^{n\times r}\rightarrow\mathbb{R}$ such that
$F(R)=-\left\langle C,RR^{\top}\right\rangle$ is the objective function of
(SQKELR). We have that

$\displaystyle F\left({\rm
Proj}_{\mathcal{K}_{n,r}}\big{(}\widehat{R}\left(t,H\right)\big{)}\right)=F\left(\widehat{R}\left(t,H\right)\right)+o(t^{2})$
$\displaystyle=F(R)+t^{2}\big{\langle}\left(2{\rm
diag}\left(\left(Cv\right)\circ d\right)-C\right)H,\,H\big{\rangle}+o(t^{2}),$
(27)

where the second equality is obtained by expanding
$\widehat{R}\left(t,H\right)\widehat{R}\left(t,H\right)^{\top}$ and some
algebraic manipulations. From (25) and (4.2), we have that for $t>0$
sufficiently small,

$F\left({\rm
Proj}_{\mathcal{K}_{n,r}}\left(\widehat{R}\left(t,H\right)\right)\right)<F(R).$
(28)

This contradicts to the assumption that $R$ is a local minimum of (SQKELR).
Hence the lemma is proved.

Lemma 4.3 provides us a necessary condition for a non-regular point of
(SQKELR) to be a local minimum. The equality (4.2) also provides us a way to
escape from any non-regular point that does not satisfy the necessary
condition. Next, we give the following definition.

######  Definition 4.4.

Suppose $(a,\tau)$ satisfies Assumption 1 and $r\geq 2.$ Let
$R=ve_{1}^{\top}\in\mathcal{N}_{n,r},$ where $v\in\\{0,1\\}^{n}$ and
$d\in\\{-1,1\\}^{n},\sigma_{v}\in\\{-1,1\\}$ are defined as in Proposition
3.3. We say that $H\in\mathbb{R}^{n\times(r-1)}$ is an escaping direction of
$R$ if $H$ is a feasible solution of (24) with a negative function value.

Obviously, we can escape from a non-regular point of (SQKELR) by moving along
an escaping direction of it. However, if there is no escaping direction, what
can we tell about a non-regular point $R=ve_{1}^{\top}\in\mathcal{N}_{n,r}$?
To answer the question, we first consider the following SDP relaxation of
problem (24)

$\min\left\\{\left\langle 2{\rm diag}\left(\left(Cv\right)\circ
d\right)-C,X\right\rangle:\ \begin{array}[]{l}\left\langle
aa^{\top}-\sigma_{v}\tau\cdot{\rm diag}\left(a\circ
d\right),X\right\rangle=0,\\\\[5.0pt] \mathrm{Tr}(X)=1,\
X\in\mathbb{S}^{n}_{+}\end{array}\right\\}.$ (29)

The dual problem of (29) is as follows:

$\max\left\\{\beta:\ 2{\rm diag}\left(\left(Cv\right)\circ
d\right)-C-\alpha\left(aa^{\top}-\sigma_{v}\tau\cdot{\rm diag}\left(a\circ
d\right)\right)-\beta
I\in\mathbb{S}^{n}_{+},\alpha,\beta\in\mathbb{R}\right\\}.$ (30)

We show in Lemma A.2 that the duality gap between (29) and (30) is zero and
the optimal values of both (29) and (30) are attainable, i.e., there exists a
KKT solution of (29).

Now, we can define the second order stationarity of a non-regular point as
follows.

######  Definition 4.5.

Suppose $(a,\tau)$ satisfies Assumption 1 and $r\geq 3.$ Let
$R=ve_{1}^{\top}\in\mathcal{N}_{n,r},$ where $v\in\\{0,1\\}^{n}$ and
$d\in\\{-1,1\\}^{n},\sigma_{v}\in\\{-1,1\\}$ are defined as in Proposition
3.3. We say that $R$ is a second order stationary point of (SQKELR) if the
optimal function value of the SDP (29) is nonnegative.

######  Remark 4.6.

Verifying the second order stationarity of a non-regular point of (SQKELR) is
equivalent to solving the SDP problem (29), which can be done in polynomial
time. Since there are only two affine constraints in (29), there exists an
optimal solution $X^{*}$ such that ${\rm rank}\left(X^{*}\right)\leq 2.$
Observe that if $\big{\langle}\left(2{\rm diag}\left(\left(Cv\right)\circ
d\right)-C\right),\,X^{*}\big{\rangle}<0$ and $X^{*}=HH^{\top},$ where
$H\in\mathbb{R}^{n\times(r-1)},$ then from Definition 4.4, we have that $H$ is
an escaping direction and hence $R=ve_{1}^{\top}$ is not a local minimum for
(SQKELR).

With Definition 4.5, we have the following simple corollary.

######  Corollary 4.7.

Suppose $a,\tau$ satisfy Assumption 1 and $r\geq 3.$ Let
$R=ve_{1}^{\top}\in\mathcal{N}_{n,r},$ where $v\in\\{0,1\\}^{n}$ and
$d\in\\{-1,1\\}^{n},\sigma_{v}\in\\{-1,1\\}$ are defined as in Proposition
3.3. If $R$ is a local minimum of (SQKELR), then $R$ is a second order
stationary point.

###  4.3 When is a second order stationary point of (SQKELR) global optimal.

In this subsection, with the definition of second order stationarity, we move
on to study when is a second order stationary point of (SQKELR) also a global
optimal solution. The following lemma gives a sufficient condition for a
second order stationary point of (SQKELR) to be a global optimal solution.

######  Lemma 4.8.

Suppose $(a,\tau)$ satisfies Assumption 1 and $r\geq 3.$ If
$R\in\mathcal{K}_{n,r}$ is a second order stationary point of (SQKELR) such
that ${\rm rank}\left(R\right)\leq r-2$, then
$Y:=\begin{pmatrix}e_{1}&R^{\top}\end{pmatrix}^{\top}\begin{pmatrix}e_{1}&R^{\top}\end{pmatrix}$
is a global optimal solution of (SQKE) and so $R$ is a global optimal solution
of (SQKELR).

###### Proof.

Proof. Case 1. $R$ is regular.
Since $R\in\mathcal{K}_{n,r}$ is a regular point, the LICQ holds at $R.$ Thus,
it satisfies the following first order and second order KKT conditions for
(SQKELR), i.e., there exists
$\left(\mu,\lambda\right)\in\mathbb{R}^{n}\times\mathbb{R}$ such that

$\displaystyle\|a^{\top}R\|^{2}-\tau a^{\top}Re_{1}=0,\ {\rm
diag}\big{(}RR^{\top}\big{)}-Re_{1}=0,$ (31) $\displaystyle-2CR-2{\rm
diag}(\mu)R+\mu e_{1}^{\top}-2\lambda aa^{\top}R+\lambda\tau ae_{1}^{\top}=0,$
(32) $\displaystyle\big{\langle}-2C-2{\rm diag}\left(\mu\right)-2\lambda
aa^{\top},\,HH^{\top}\big{\rangle}\geq 0\ \forall
H\in\mathcal{T}_{\mathcal{K}_{n,r}}\left(R\right).$ (33)

Because LICQ is satisfied at $R,$
$\mathcal{T}_{\mathcal{K}_{n,r}}\left(R\right)$ corresponds to the following
tangent space of $\mathcal{K}_{n,r}$ at $R$

$T_{R}\mathcal{K}_{n,r}:=\left\\{h\in\mathbb{R}^{n\times r}:\
2\big{\langle}aa^{\top}R,\,H\big{\rangle}-\tau a^{\top}He_{1}=0,\ 2{\rm
diag}\big{(}RH^{\top}\big{)}-He_{1}=0\right\\}.$ (34)

Since ${\rm rank}\left(R\right)\leq r-2,$ there exists $b\in\mathbb{R}^{r}$
such that $b\neq 0,$ $b_{1}=0$ and $Rb=0.$ It is easy to see that for any
$d\in\mathbb{R}^{n},$ $db^{\top}\in T_{R}\mathcal{K}_{n,r}.$ Substituting
$db^{\top}$ into (33), we get that

$-2C-2{\rm diag}\left(\mu\right)-2\lambda aa^{\top}\in\mathbb{S}^{n}_{+}.$
(35)

Now, we will move on to use (31), (32) and (35) to prove the global optimality
of $R.$ Note that we only have to prove that $Y=\left(e_{1}\
R^{\top}\right)^{\top}\left(e_{1}\ R^{\top}\right)$ is a global optimal
solution of (SQKE) through the KKT conditions of (SQKE). First rearrange (32)
to get the following equation

$\left(\frac{\mu+\lambda\tau a}{2},-C-{\rm diag}\left(\mu\right)-\lambda
aa^{\top}\right)\begin{pmatrix}e_{1}^{\top}\\\ R\end{pmatrix}=0.$ (36)

After multiplying $R^{\top}$ to the left-hand side and $e_{1}$ to the right-
hand side of (32), we get the following equation

$R^{\top}\left(-C-{\rm diag}\left(\mu\right)-\lambda
aa^{\top}\right)Re_{1}+\frac{1}{2}R^{\top}\mu+\frac{1}{2}\lambda\tau
R^{\top}a=0.$ (37)

Dividing equation (32) by 2 and taking its transpose, we get the following
equation

$R^{\top}\left(-C-{\rm diag}\left(\mu\right)-\lambda
aa^{\top}\right)=-\frac{e_{1}\mu^{\top}+\lambda\tau e_{1}a^{\top}}{2}.$ (38)

Substituting (38) into (37) and taking its transpose, we get

$\frac{-\left(\mu^{\top}+\lambda\tau
a^{\top}\right)Re_{1}}{2}e_{1}^{\top}+\frac{\mu^{\top}+\lambda\tau
a^{\top}}{2}R=0.$ (39)

Combining (36) and (39), we get

$\begin{pmatrix}\frac{-\left(\mu^{\top}+\lambda\tau
a^{\top}\right)Re_{1}}{2}&\frac{\mu^{\top}+\lambda\tau a^{\top}}{2}\\\\[5.0pt]
\frac{\mu+\lambda\tau a}{2}&-C-{\rm diag}\left(\mu\right)-\lambda
aa^{\top}\end{pmatrix}\begin{pmatrix}e_{1}^{\top}\\\ R\end{pmatrix}=0.$ (40)

Note that matrix
$Y:=\begin{pmatrix}e_{1}&R^{\top}\end{pmatrix}^{\top}\begin{pmatrix}e_{1}&R^{\top}\end{pmatrix}$
is inside the feasible region of (SQKE). Now consider the problem (SQKE). If
we consider $(0,0\,;\,0,-C)$ as its profit matrix, $\mu$ as the dual variable
of the constraint ${\rm diag}\left(X\right)=x,$ $\lambda$ as the dual variable
of the constraint $a^{\top}Xa-\tau a^{\top}x=0$ and
$\frac{1}{2}\left(\mu^{\top}+\lambda\tau a^{\top}\right)Re_{1}$ as the dual
variable for $Y_{11}=1,$ then (43) implies that the complementarity condition
of (SQKE) hold for the primal variable $Y$ and the dual variables just
mentioned. In order to prove that $Y$ is an optimal solution of (SQKE), we
only have to check its dual feasibility i.e., whether the dual slack ${\bf S}$
of (SQKE) is positive semidefinite, where

${\bf S}:=\begin{pmatrix}\frac{-\left(\mu^{\top}+\lambda\tau
a^{\top}\right)Re_{1}}{2}&\frac{\mu^{\top}+\lambda\tau a^{\top}}{2}\\\\[5.0pt]
\frac{\mu+\lambda\tau a}{2}&-C-{\rm diag}\left(\mu\right)-\lambda
aa^{\top}\end{pmatrix}.$ (41)

The positive semidefiniteness of ${\bf S}$ comes from (35) and the following
equation

${\bf S}=\begin{pmatrix}-\left(Re_{1}\right)^{\top}\\\
I\end{pmatrix}\left(-C-{\rm diag}\left(\mu\right)-\lambda
aa^{\top}\right)\begin{pmatrix}-Re_{1}&I\end{pmatrix},$ (42)

which can be verified by using (38) several times.
Case 2. $R$ is nonregular.
In this case, suppose $R=ve_{1}^{\top}\in\mathcal{N}_{n,r},$ where
$v\in\\{0,1\\}^{n}$ and $d\in\\{-1,1\\}^{n},\sigma_{v}\in\\{-1,1\\}$ are
defined as in Proposition 3.3. From Definition 4.5, we have that the optimal
value of (29) is nonnegative. From Lemma A.2, there exists a KKT solution
$\left(X,\alpha,\gamma\right)$ of (29) that satisfies the following KKT
condition

$\displaystyle\left\langle aa^{\top}-\sigma_{v}\tau\cdot{\rm diag}\left(a\circ
d\right),X\right\rangle=0,\ \mathrm{Tr}\left(X\right)=1,\
X\in\mathbb{S}^{n}_{+}$ (43) $\displaystyle 2{\rm
diag}\left(\left(Cv\right)\circ
d\right)-C-\alpha\left(aa^{\top}-\sigma_{v}\tau\cdot{\rm diag}\left(a\circ
d\right)\right)-\gamma I\in\mathbb{S}^{n}_{+}$ (44) $\displaystyle\left\langle
2{\rm diag}\left(\left(Cv\right)\circ
d\right)-C-\alpha\left(aa^{\top}-\sigma_{v}\tau\cdot{\rm diag}\left(a\circ
d\right)\right)-\gamma I,X\right\rangle=0.$ (45)

Also, since the optimal function value of (29) is nonnegative, we have that

$\left\langle 2{\rm diag}\left(\left(Cv\right)\circ
d\right)-C,X\right\rangle\geq 0.$ (46)

Now, combining (43), (45) and (46), we have that $\gamma\geq 0.$ Then, from
(44), we have that

$2{\rm diag}\left(\left(Cv\right)\circ
d\right)-C-\alpha\left(aa^{\top}-\sigma_{v}\tau\cdot{\rm diag}\left(a\circ
d\right)\right)\in\mathbb{S}^{n}_{+}.$ (47)

Comparing (47) and (35), we may choose the dual variables of (SQKELR) as
follows

$\mu:=-2\left(Cv\right)\circ d-\alpha\sigma_{v}\tau\cdot a\circ
d,\quad\lambda:=\alpha.$ (48)

Then from (47), (35) holds. Now, in order to prove
$Y:=\begin{pmatrix}e_{1}&R^{\top}\end{pmatrix}^{\top}\begin{pmatrix}e_{1}&R^{\top}\end{pmatrix}$
is the global optimal solution of (SQKE), we only have to verify whether (32)
holds for $R,$ $\mu$ and $\lambda.$ Because once (32) holds, the remaining
proof will be exactly the same as in Case 1. We have that

$\displaystyle-2CR-2{\rm diag}\left(\mu\right)R+\mu e_{1}^{\top}-2\lambda
aa^{\top}R+\lambda\tau ae_{1}^{\top}$ $\displaystyle=\left(-2Cv-2\mu\circ
v+\mu-2\lambda aa^{\top}v+\lambda\tau a\right)e_{1}^{\top}$
$\displaystyle=\left(-2Cv-2\left(-2\left(Cv\right)\circ
d-\alpha\sigma_{v}\tau\cdot a\circ d\right)\circ v-2\left(Cv\right)\circ
d-\alpha\sigma_{v}\tau\cdot a\circ d-2\alpha aa^{\top}v+\alpha\tau
a\right)e_{1}^{\top}$ $\displaystyle=\Big{(}-2Cv+4(Cv)\circ
v-2\left(Cv\right)\circ
d\Big{)}e_{1}^{\top}+\tau\alpha\Big{(}2\sigma_{v}a\circ v-\sigma_{v}a\circ
d-\frac{2a^{\top}v}{\tau}a+a\Big{)}e_{1}^{\top}$
$\displaystyle=\tau\alpha\Big{(}\sigma_{v}a-\frac{2a^{\top}v}{\tau}a+a\Big{)}e_{1}^{\top},$
(49)

where the last equality comes from that $d=2v-e.$ If $v\not=0$, then
$\sigma_{v}=1,$ together with the fact that $v^{\top}a=\tau$, we get
$\sigma_{v}a-\frac{2a^{\top}v}{\tau}a+a=0.$ Also, if $v=0$, then
$\sigma_{v}=-1,$ and $\sigma_{v}a-\frac{2a^{\top}v}{\tau}a+a=0.$ Thus, we have
verified (32).

######  Remark 4.9.

Compared with the results in [9, 10], our rank deficient condition is slightly
stronger, i.e., there exists at least two zero singular values of a certain
second order stationary point. One reason is that the first row of the
factorization of $Y$ is fixed to be $e_{1}^{\top}$. Also, note that the proof
of Lemma 4.8 provides a way to recover the dual variable of (SQKELR), which is
also the dual variable for (SQKE). When $R$ is regular, we can solve (32) to
get a least square solution $(\mu,\lambda).$ The corresponding LICQ property
implies that there exists a unique least square solution of (32), which is
exactly the dual variable. When $R$ is non-regular, then (48) is an explicit
formula for the dual variable of (SQKELR).

Now, we move on to study when will the rank deficient condition in Lemma 4.8
hold. The following result is inspired by Lemma 3 in [4] and Theorem 2 in [9].
However, our rank-bound is also related to the minimum rank number of the
underlying graph of $C$ in (SQKELR).

######  Lemma 4.10.

Suppose $(a,\tau)$ satisfies Assumption 1, $r\geq 3$ and $R$ is a regular
first order stationary point of (SQKELR). The following statements hold.

  * (i)

If $\frac{(r-2)(r-1)}{2}>n+1,$ then for a generic $C\in\mathbb{S}^{n},$ ${\rm
rank}\left(R\right)\leq r-2.$

  * (ii)

If $C\in\mathcal{Q}(G)$ for some graph $G$ and $r\geq n+4-{\rm mr}(G),$ then
${\rm rank}\left(R\right)\leq r-2.$

###### Proof.

Proof. Because $R$ is a regular first order stationary point of (SQKELR), we
have that the KKT condition (32) and its variant (36) hold at $R.$ After
applying Sylvester’s rank inequality (see 10.5 of [30]) to (36), we have that

${\rm rank}\left(-C-{\rm diag}\left(\mu\right)-\lambda aa^{\top}\right)\leq
n+1-{\rm rank}\left(R\right).$ (50)

Now, suppose ${\rm rank}\left(R\right)\geq r-1$ and
$\frac{(r-2)(r-1)}{2}>n+1,$ then from (50), we have that

${\rm rank}\left(-C-{\rm diag}\left(\mu\right)-\lambda aa^{\top}\right)\leq
n+2-r,$ (51)

which implies that $-C-{\rm diag}\left(\mu\right)-\lambda
aa^{\top}\in\mathcal{M}_{\leq n+2-r},$ where

$\mathcal{M}_{\leq n+2-r}:=\left\\{A\in\mathbb{S}^{n}:\ {\rm
rank}\left(A\right)\leq n+2-r\right\\}.$

Then we have that $-C\in\mathcal{M}_{\leq n+2-r}+\left\\{{\rm
diag}\left(\mu\right)+\lambda aa^{\top}:\ \mu\in\mathbb{R}^{n},\
\lambda\in\mathbb{R}\right\\}$ whose dimension is bounded by
$\frac{n(n+1)}{2}-\frac{(r-2)(r-1)}{2}+n+1<\frac{n(n+1)}{2}=\dim\left(\mathbb{S}^{n}\right).$
This means that $C\in\mathbb{S}^{n}$ is contained in a set of measure zero.
Thus, (i) is proved.

If $r\geq n+4-{\rm mr}(G)$, then from (50), we have that

$\displaystyle{\rm rank}\left(R\right)\leq n+1-{\rm rank}\left(-C-{\rm
diag}\left(\mu\right)-\lambda aa^{\top}\right)$ $\displaystyle\leq n+2-{\rm
rank}\left(-C-{\rm diag}\left(\mu\right)\right)\leq n+2-{\rm mr}(G)\leq r-2.$
(52)

Thus, (ii) is verified.

Note that we did not consider non-regular points of $\mathcal{K}_{n,r}$ in
Lemma 4.10. This is because from Proposition 3.1, the rank of any non-regular
point $R$ of $\mathcal{K}_{n,r}$ already satisfies the condition that ${\rm
rank}\left(R\right)\leq 1\leq r-2$ for any $r\geq 3.$ Now, we are able to
state our main result.

######  Theorem 4.11.

Suppose $(a,\tau)$ satisfies Assumption 1 and $r\geq 3.$ If $R$ is a second
order stationary point of (SQKELR) and at least one of the following two
conditions are satisfied

  * 1.

$\frac{(r-2)(r-1)}{2}>n+1$ and $C\in\mathbb{S}^{n}$ is generic;

  * 2.

$C\in\mathcal{Q}(G)$ for some graph $G$ and $r\geq n+4-{\rm mr}(G),$

then
$Y:=\begin{pmatrix}e_{1}&R^{\top}\end{pmatrix}^{\top}\begin{pmatrix}e_{1}&R^{\top}\end{pmatrix}$
is a global optimal solution of (SQKE) and so $R$ is a global optimal solution
of (SQKELR).

###### Proof.

Proof. Theorem 4.11 follows directly from the combination of Lemma 4.8 and
Lemma 4.10.

##  5 Algorithms and implementation details.

In this section, we consider how to design a feasible algorithm to solve
(SQKELR) by using the results established in the last few sections. We discuss
several techniques involved in our algorithm.

###  5.1 Projection and retraction.

As mentioned in Corollary 3.2, $\mathcal{K}_{n,r}$ is a smooth manifold for a
generic input data $(a,\tau)$, except for a trivial point. Also, even if
$a\in\mathcal{W}_{n,\tau},$ which is defined in (9), from (4), there are only
finitely many non-regular points in $\mathcal{K}_{n,r}.$ Therefore, we can use
a Riemannian optimization method to solve (SQKELR) (see [7]). Two important
operations in a Riemannian optimization method are projection and retraction.
Consider a regular point $R\in\mathcal{K}_{n,r}.$ Because LICQ holds at $R,$
$\mathcal{K}_{n,r}$ is locally a smooth manifold around $R.$ The projection
mapping of $\mathcal{K}_{n,r}$ at $R$ is simply the orthogonal projection
operator of its tangent space $T_{R}\mathcal{K}_{n,r}.$ However, the
retraction mapping (see definition 3.41 and definition 5.41 in [7] for the
definition of retraction and second order retraction)
$\mathcal{R}_{R}:\mathcal{T}_{R}\mathcal{K}_{n,r}\rightarrow\mathcal{K}_{n,r}$
is not unique. For simple manifolds like the unit sphere ${\rm S}^{r-1}$, its
retraction can be defined as a normalisation mapping
$x\rightarrow\frac{x}{\|x\|_{2}}$. This retraction has an explicit formula and
is also the closest point retraction. For more complicated manifold, it may be
difficult for us to find a retraction mapping with a simple closed form. In
this case, we have to compute the retraction iteratively (see [1, 42]). In our
algorithm, we use the Newton retraction introduced in [42]. This retraction is
equivalent to applying the Gauss-Newton method to the nonlinear system of the
manifold and it is initialized by the point at which we want to retract (see
algorithm 1 of [42]). It has been proved that this is a second order
retraction in [42].

###  5.2 A Riemannian optimization approach that can escape from non-optimal
singular points.

In this subsection, we suppose $(a,\tau)$ satisfies Assumption 1 and $r\geq
3.$ When we meet a non-regular point $R=ve_{1}^{\top}$ of $\mathcal{K}_{n,r},$
we may first solve (29). If the optimal value of (29) is negative, then from
Remark 4.6, we can factorize the optimal solution to get an escaping
direction. After that we can move along this direction to escape from this
non-optimal non-regular point while reducing the function value of (SQKELR).
If the optimal value of (29) is nonnegative, i.e., $R$ is a non-regular second
order stationary point, then from Lemma 4.8, we have that
$Y:=\begin{pmatrix}e_{1}&R^{\top}\end{pmatrix}^{\top}\begin{pmatrix}e_{1}&R^{\top}\end{pmatrix}$
is an optimal solution of (SQKE) and so $R$ is a global optimal solution of
(SQKELR). In practice, it may happen that the iterations of a Riemannian
optimization method slowly approach a non-regular point of
$\mathcal{K}_{n,r}$. Thus, to mitigate the slow convergence towards a non-
regular point, we have to apply rounding to every nearly non-regular point to
a non-regular point. Due to the constraint ${\rm
diag}\left(RR^{\top}\right)=Re_{1},$ we have that $Re_{1}\in[0,1]^{n}.$ Define
the rounding function
$\texttt{Round}:\mathcal{K}_{n,r}\rightarrow\\{0,1\\}^{n}e_{1}^{\top}$ such
that for any $R\in\mathcal{K}_{n,r},$ $\texttt{Round}(R)=ve_{1}^{\top}$ and

$v_{i}=\begin{cases}1&{\rm if}\ \left(Re_{1}\right)_{i}\geq 0.5\\\ 0&{\rm
otherwise}\end{cases}.$ (53)

With the function $\texttt{Round}(\cdot),$ we may define the following set for
any $\delta>0.$

$\mathcal{K}_{n,r}\left(\delta\right):=\left\\{R\in\mathcal{K}_{n,r}:\
\|R-\texttt{Round}(R)\|<\delta\right\\}.$ (54)

It is easy to see that $\mathcal{N}_{n,r}\in\mathcal{K}_{n,r}(\delta)$ for any
$\delta>0.$

######  Lemma 5.1.

Suppose $(a,\tau)$ satisfies Assumption 1 and $r\geq 3.$ Then there exists
$\delta_{0}>0$ such that for any $\delta\in(0,\delta_{0})$ and any
$R\in\mathcal{K}_{n,r}\left(\delta\right),$
$\texttt{Round}(R)\in\mathcal{N}_{n,r}$.

###### Proof.

Proof. Assume on the contrary that there exists $\delta_{k}\downarrow 0$ such
that for any $k\in\mathbb{N}^{+},$ there exists
$R_{k}\in\mathcal{K}_{n,r}\left(\delta_{k}\right)$ such that
$\texttt{Round}\left(R_{k}\right)=v_{k}e_{1}^{\top}\notin\mathcal{N}_{n,r}.$
This implies that

$\|a^{\top}v_{k}\|^{2}-\tau a^{\top}v_{k}\neq 0,\ \forall k\in\mathbb{N}^{+}.$
(55)

Since there are only finitely many elements in $\\{\|a^{\top}v\|^{2}-\tau
a^{\top}v:\ v\in\\{0,1\\}^{n}\\},$ we have that there exists $\epsilon>0$ such
that

$\left|\|a^{\top}v_{k}\|^{2}-\tau a^{\top}v_{k}\right|\geq\epsilon,\ \forall
k\in\mathbb{N}^{+}.$ (56)

However, from $R_{k}\in\mathcal{K}_{n,r}(\delta_{k}),$ we have that
$\|a^{\top}R_{k}\|^{2}-\tau a^{\top}R_{k}e_{1}=0$ and
$\lim\limits_{k\rightarrow\infty}\|R_{k}-v_{k}e_{1}^{\top}\|=0.$ Thus, we have
that

$\displaystyle\lim\sup\limits_{k\rightarrow\infty}\left|\|a^{\top}v_{k}\|^{2}-\tau
a^{\top}v_{k}\right|=\lim\sup\limits_{k\rightarrow\infty}\left|\|a^{\top}v_{k}e_{1}^{\top}\|^{2}-\tau
a^{\top}v_{k}e_{1}^{\top}e_{1}\right|$
$\displaystyle=\lim\sup\limits_{k\rightarrow\infty}\left|\|a^{\top}R_{k}\|^{2}-\tau
a^{\top}R_{k}e_{1}\right|=0,$ (57)

which contradicts to (56).

With Lemma 5.1, we can design the following algorithm, which can be seen as a
modified version of algorithm 1 in [34].

Algorithm 1:

Define $f\left(R\right):=\left\langle-C,RR^{\top}\right\rangle.$ Choose
$R_{0}\in\mathcal{K}_{n,r}$ and $\delta_{0}\in(0,1).$ Set $k=0.$

Step 1.

Set $R_{k}$ as the initial point, $f_{k}:=f(R_{k}).$ Use any Riemannian
optimization method to solve (SQKELR) such that all the generated iteration
point $R^{(j)}\in\mathcal{K}_{n,r}$ ($\forall j\in\\{0,1,2,\ldots\\}$, with
$R^{(0)}=R_{k}$) satisfies $f(R^{(j)})\leq f_{k}$. If the KKT residue
corresponding to $R^{(j)}$ as defined in (6) for the convex SDP (SQKE)
satisfies the required accuracy tolerance, terminate the Riemannian
optimization method and stop. If at some iteration of the Riemannian
optimization method, we encounter an iterate denoted by $\widehat{R}_{k}$ that
satisfies $\widehat{R}_{k}\in\mathcal{K}_{n,r}(\delta_{k}),$ goto Step 2.

Step 2.

Case 1. $\texttt{Round}\big{(}\widehat{R}_{k}\big{)}\notin\mathcal{K}_{n,r}$:
set $R=\widehat{R}_{k},$ goto step 2’.
Case 2. $\texttt{Round}\big{(}\widehat{R}_{k}\big{)}\in\mathcal{K}_{n,r}$ and
it is an optimal solution: Stop.
Case 3. $\texttt{Round}\big{(}\widehat{R}_{k}\big{)}\in\mathcal{K}_{n,r}$ and
is non-optimal: Use the escaping strategy described above to find another
regular point (regularity comes from the finiteness of non-regular points of
$\mathcal{K}_{n,r}$) $R^{+}\in\mathcal{K}_{n,r}$ such that
$f(R^{+})<f(\texttt{Round}(\widehat{R}_{k}))-\epsilon\big{(}\texttt{Round}\big{(}\widehat{R}_{k}\big{)}\big{)},$
where $\epsilon\big{(}\texttt{Round}\big{(}\widehat{R}_{k}\big{)}\big{)}>0$ is
only related to $\texttt{Round}\big{(}\widehat{R}_{k}\big{)}.$ If
$f(R^{+})<f(\widehat{R}_{k})$ set $R\leftarrow R^{+};$ Otherwise, set
$R\leftarrow\widehat{R}_{k}.$ (In the latter case, $R$ is also a regular
point. This is because otherwise
$f(R^{+})<f(\texttt{Round}(\widehat{R}_{k}))=f(\widehat{R}_{k})$ and this case
won’t happen.) Goto Step 2’.

Step 2’.

Set $\delta_{k+1}\leftarrow\delta_{k}/2$, $R_{k+1}\leftarrow R$, $k\leftarrow
k+1,$ goto Step 1.

######  Remark 5.2.

Similar to Theorem 4.1 of [34], we can show that for any initialization,
Algorithm 1 will terminate after finitely many outer iterations. Here we omit
the proof because of their similarity. In other words, after finitely many
escaping steps, the iterations of Algorithm 1 will either terminate at some
optimal non-regular point or enter a subset
$\mathcal{K}_{n,r}\setminus\mathcal{K}_{n,r}\left(\delta_{k}\right)$ of
$\mathcal{K}_{n,r}$ that is bounded away from the non-regular points of
$\mathcal{K}_{n,r}.$ In this case, Algorithm 1 is reduced to optimization on a
smooth manifold and any results in Riemannian optimization can be applied in
our algorithm.

In Step 1 of Algorithm 1, we choose a Riemannian gradient descent method with
Barzilai-Borwein (BB) step and non-monotone linesearch (see [3, 19, 25, 33]).
We omit the details of the algorithm because it is frequently used in the
literature. The parameters we choose for BB step is the same as the Riemannian
BB algorithm in [34].

###  5.3 A rounding procedure to the output of (SQKELR).

In practice, we can apply rounding procedures to the output of the SDP
relaxations of integer programming problems to get an approximate solution. A
famous example is the SDP relaxation of a max-cut problem. In [21], Goemans
and Williamson designed a randomized rounding procedure for the output of a
max-cut SDP whose expected value is at least 0.878 of the optimal value of the
max-cut problem. In this subsection, we describe a simple rounding procedure
for the output of the problem (SQKELR). Suppose $R\in\mathbb{R}^{n\times r}$
is the optimal solution of (SQKELR) and $d=(d_{1},d_{2},\ldots,d_{n})^{\top}$
is the first column of $R.$ We first sort the entries of $d$ in descending
order i.e., $d_{i_{1}}\geq d_{i_{2}}\geq\ldots\geq d_{i_{n}}.$ Let $k\in[n]$
to be the maximum integer such that $\sum_{\ell=1}^{k}a_{i_{\ell}}\leq\tau.$
Then, define the output $x\in\\{0,1\\}^{n}$ of our rounding procedure as
follows: $x_{i_{\ell}}=1$ if $\ell\leq k$ and $x_{i_{\ell}}=0$ if $\ell>k.$ It
is easy to see that $x$ is a feasible solution for (QKP). Also, the complexity
of our rounding procedure is $O(n\log n)$ because of the sorting step. The
(SQK) model together with the rounding procedure can be considered as a
heuristic to find a lower bound of (QKP). We can measure the accuracy of the
lower bound by comparing it with the upper bound provided by the SDP
relaxation (SQK). (Here We remove the minus sign of the objective function of
(SQK) to get the upper bound of (QKP).) We will test the accuracy and
efficiency of this method later in section 6.3.

##  6 Numerical experiments.

In this section, we conduct numerical experiments to verify the efficiency of
our feasible method described in the previous section. Our dataset satisfies
Assumption 1 and the profit matrix $C$ is always nonnegative. Then, from
Proposition 2.1, (SQK) is equivalent to (SQKE) so we directly test the problem
(SQKE). We use the KKT residue (given by $\max\\{{\rm Rp},{\rm Rd},{\rm
pdgap}\\}$) for (SQKE) to measure the accuracy of a solution
$Y=\left(\alpha,x^{\top};x,X\right)\in\mathbb{S}^{n+1}_{+}$, where

$\displaystyle{\rm primal\ residue:}\ {\rm Rp}:=\frac{\sqrt{\|{\rm
diag}\left(X\right)-x\|^{2}+(\alpha-1)^{2}+\left(a^{\top}Xa-\tau
a^{\top}x\right)^{2}}}{2},$ $\displaystyle{\rm dual\ residue:}\ {\rm
Rd}:=\frac{\|\Pi_{\mathbb{S}^{n+1}_{-}}\left({\bf S}\right)\|}{1+\|{\bf
S}\|},\ {\rm duality\ gap:}\ {\rm
pdgap:}=\frac{|\left\langle-C,X\right\rangle-y|}{1+|\left\langle-C,X\right\rangle|+|y|},$
(58)

where ${\bf S}$ is the dual slack of (SQKE) defined in (41), and the dual
variable $\left(\mu,\lambda\right)$ can be constructed from Remark 4.9; $y$ is
the dual variable for the constraint $Y_{11}=1.$ For our method,
$y=\frac{1}{2}\left(\mu^{\top}+\lambda\tau a^{\top}\right)Re_{1}.$ Note that
Algorithm 1 uses factorization so we do not have to check whether $Y$ is
positive semidefinite.

The algorithms we compare are SDPLR by Burer and Monteiro and SketchyCGAL by
Yurtsever et al. in [41]. Their source codes can be downloaded from
https://sburer.github.io/projects.html and
https://github.com/alpyurtsever/SketchyCGAL, respectively. Note that all three
algorithms make use of the low rank property of the linear SDP to save
storage. Therefore, our comparison will demonstrate the advantage of
exploiting the geometric property of a problem. We will stop the Riemannian
optimization method in Step 1 of Algorithm 1 when the KKT residue $\max\\{{\rm
Rp},\ {\rm Rd},\ {\rm pdgap}\\}<10^{-6}.$ Because it is expensive to compute
the Rd too frequently, we only check the KKT residue when the normalized norm
of Riemannian gradient is small enough, i.e.,

$\frac{\|{\rm
Proj}_{T_{R}\mathcal{K}_{n,r}}\left(-2CR\right)\|}{\max\\{1,\|R\|\\}}<{\rm
tolg}.$ (59)

We first set tolg $:=10^{-6}.$ When (59) is satisfied but $\max\\{{\rm Rp},\
{\rm Rd},\ {\rm pdgap}\\}\geq 10^{-6},$ we will update the tolerance as ${\rm
tolg}\leftarrow{\rm tolg}/10$ and continue running our algorithm. Note that we
will always check the KKT residue for a non-regular point because its
Riemannian gradient may not be well-defined.

We also set the stopping tolerance for SDPLR to be $10^{-6}.$ For SketchyCGAL,
we find that it cannot find an accurate solution, so we set the tolerance to
be $10^{-1}.$ The parameter $\delta_{0}$ in Algorithm 1 is chosen to be $0.1.$
For Algorithm 1, we use random initialization with normal distribution. For
SDPLR and SketchyCGAL, we use their defaulted initialization. From our
experiment, SDPLR and SketchCGAL usually fail to converge without a pre-scaled
input data. Therefore, in order to improve their performance, we scale the
capacity $\tau$ and weight vector $a$ as $a\leftarrow a/\tau$ and
$\tau\leftarrow 1.$ For all three algorithms, we set the maximum running time
to be 3600 seconds. All the experiments are run using Matlab R2021b on a
Workstation with a Intel(R) Xeon(R) CPU E5-2680 v3 @ 2.50GHz Processor and
128GB RAM.

###  6.1 The SDP relaxation of Binary linear knapsack problem.

In this subsection, we consider binary linear knapsack problems. In this case,
$C$ is a diagonal matrix. We use the randomly generated data in [31]. The
dataset can be downloaded from http://artemisa.unicauca.edu.co/
johnyortega/instances_01_KP/. In order to verify the new rank-bound from
Theorem 4.1, we choose $r=3$ in Algorithm 1. Note that smaller $r$ does not
always lead to faster convergence because it may introduce more spurious local
minima. From our numerical tests, we find that the behaviours of SDPLR and
SketchyCGAL will always be better if we choose $r$ to be slightly larger than
$3.$ Therefore, we choose their rank to be $r=10$ to give them some advantage
over Algorithm 1 while the problem size is still $O(n)$. Note that we have
also tested the setting $r=20$, but the performance of SDPLR and SketchyCGAL
are rather similar to that for the setting $r=10$.

Table 1: Comparison of Algorithm1, SDPLR and SketchyCGAL for binary linear knapsack SDP. knapPI_1, knapPI_2 and knapPI_3 are uncorrelated instances, weakly correlated instances and strongly correlated instances, respectively (see [31]).  problem | algorithm | obj | Rp | Rd | pdgap | time
---|---|---|---|---|---|---
knapPI_1 | Algorithm1 | -5.4538020e+04 | 1.64e-10 | 5.38e-09 | 1.46e-08 | 2.89e-01
n = 1000 | SDPLR | -5.4537907e+04 | 9.83e-07 | 2.48e-06 | 5.59e-05 | 9.47e+00
opt = 54503 | SketchyCGAL | -5.1024775e+04 | 6.27e-02 | 1.45e-08 | 4.77e-02 | 3.19e+02
knapPI_2 | Algorithm1 | -9.0573608e+03 | 6.46e-16 | 9.93e-10 | 5.12e-08 | 5.88e-01
n = 1000 | SDPLR | -9.0545910e+03 | 8.19e-07 | 3.24e-05 | 4.05e-04 | 1.05e+01
opt = 9052 | SketchyCGAL | -8.5724562e+03 | 9.68e-02 | 0.00e+00 | 5.15e-02 | 2.08e+03
knapPI_3 | Algorithm1 | -1.4406317e+04 | 2.23e-13 | 1.13e-07 | 7.86e-09 | 3.98e-01
n = 1000 | SDPLR | -1.4406040e+04 | 3.37e-07 | 4.73e-05 | 2.72e-04 | 2.00e+01
opt = 14390 | SketchyCGAL | -1.4364806e+04 | 1.41e-01 | 0.00e+00 | 4.96e-02 | 5.98e+02
knapPI_1 | Algorithm1 | -1.1064590e+05 | 6.64e-12 | 1.57e-08 | 2.07e-09 | 1.02e+00
n = 2000 | SDPLR | -1.1064600e+05 | 8.91e-07 | 1.27e-06 | 6.86e-05 | 8.81e+01
opt = 110625 | SketchyCGAL | -1.0143376e+05 | 1.06e-01 | 6.86e-11 | 4.71e-02 | 7.29e+02
knapPI_2 | Algorithm1 | -1.8052806e+04 | 2.13e-09 | 6.93e-07 | 3.79e-07 | 6.24e-01
n = 2000 | SDPLR | -1.8051935e+04 | 9.85e-07 | 7.43e-05 | 1.11e-05 | 1.16e+02
opt = 18051 | SketchyCGAL | -1.7201318e+04 | 1.97e-01 | 0.00e+00 | 2.02e-01 | 3.60e+03
knapPI_3 | Algorithm1 | -2.9012870e+04 | 9.18e-11 | 4.15e-10 | 5.72e-09 | 1.44e+00
n = 2000 | SDPLR | -2.9010766e+04 | 9.19e-07 | 8.41e-05 | 9.98e-05 | 1.36e+02
opt = 28919 | SketchyCGAL | -2.9631203e+04 | 1.77e-01 | 0.00e+00 | 3.78e-02 | 8.14e+02
knapPI_1 | Algorithm1 | -2.7645834e+05 | 1.71e-12 | 5.99e-07 | 2.39e-08 | 1.24e+01
n = 5000 | SDPLR | -2.7645644e+05 | 9.09e-07 | 2.85e-05 | 4.20e-05 | 1.84e+03
opt = 276457 | SketchyCGAL | -2.5881374e+05 | 2.28e-01 | 4.30e-08 | 3.35e-02 | 2.33e+03
knapPI_2 | Algorithm1 | -4.4357542e+04 | 2.45e-12 | 9.60e-08 | 3.05e-08 | 1.24e+01
n = 5000 | SDPLR | -4.4354071e+04 | 9.90e-07 | 7.06e-05 | 4.67e-04 | 1.42e+03
opt = 44356 | SketchyCGAL | -2.5213507e+04 | 6.57e-01 | 4.97e-08 | 6.07e-01 | 3.60e+03
knapPI_3 | Algorithm1 | -7.2563369e+04 | 3.51e-09 | 8.99e-10 | 9.00e-09 | 6.67e+00
n = 5000 | SDPLR | -7.2518590e+04 | 1.69e-04 | 4.47e-04 | 2.14e-05 | 2.42e+03
opt = 72505 | SketchyCGAL | -6.8704702e+04 | 5.81e-01 | 5.52e-08 | 6.02e-02 | 3.60e+03
knapPI_1 | Algorithm1 | -5.6364672e+05 | 9.87e-10 | 6.70e-07 | 2.30e-07 | 3.13e+01
n = 10000 | SDPLR | - | - | - | - | -
opt = 563647 | SketchyCGAL | - | - | - | - | -
knapPI_2 | Algorithm1 | -9.0204226e+04 | 5.55e-11 | 7.45e-07 | 1.23e-07 | 4.53e+01
n = 10000 | SDPLR | - | - | - | - | -
opt = 90204 | SketchyCGAL | - | - | - | - | -
knapPI_3 | Algorithm1 | -1.4694938e+05 | 2.50e-09 | 5.12e-08 | 1.61e-09 | 3.95e+01
n = 10000 | SDPLR | - | - | - | - | -
opt = 146919 | SketchyCGAL | - | - | - | - | -

From the results in Table LABEL:LKP, we can see that Algorithm 1 can solve all
the problems accurately. This verifies the validity of the new rank-bound for
binary linear knapsack problems. SDPLR and SketchyCGAL can only solve problems
of size up to $5000$ to a moderate and low accuracy, respectively. For
problems of size $10000,$ both of them run out of time and return a solution
that is far away from being feasible, so we do not show their output in the
table. In all these examples, Algorithm 1 is much more efficient than SDPLR
and SketchyCGAL. For some instances with $n=5000$, Algorithm 1 is more than
100 times faster than them.

###  6.2 The SDP relaxation of Quadratic knapsack problem.

In this subsection, we consider the quadratic knapsack problem (QKP) with a
quadratic profit matrix $C.$ We randomly generate the profit matrix and weight
vector following the procedure proposed by Gallo et al. in [20]. The entries
of the profit matrix $C_{ij}=C_{ji}$ are zero with probability $(1-p)$ and
otherwise integers randomly generated uniformly in the range $[1,100].$ We set
the problem density as $p=0.25.$ The coefficients of the weight vector $a_{i}$
are integers uniformly distributed in the range $[1,50].$ While the procedure
in [20] also randomly generates the knapsack capacity $\tau$ in the range
$[1,e^{\top}a],$ for convenience, we directly set the knapsack capacity to be
$\beta\cdot e^{\top}a,$ where $\beta$ is chosen to be $0.1,0.5$ and $0.9$.
This data-generating procedure has been widely accepted in the literature (see
[5, 14, 32, 37] to just name a few). For Algorithm 1, we choose
$r=\lceil\sqrt{2(n+1)}\rceil+2,$ which satisfies the rank condition in (i) of
Lemma 4.10. For SDPLR, we use their default rank. For SketchyCGAL, we choose
the rank to be $20$ because we find that the rank of the output is always
smaller than the classical theoretical upper bound of $\sqrt{2(n+2)}.$
Actually, when we set $r=20$ for Algorithm 1, it can find an accurate optimal
solution of (SQKE) even faster than the case for
$r=\lceil\sqrt{2(n+1)}\rceil+2$ (see the running time for Algorithm 1 in Table
1.1 in the next subsection). We choose $r=\lceil\sqrt{2(n+1)}\rceil+2$ for
Algorithm 1 to be consistent to Lemma 4.10. Note that we have also tested the
setting $r=50$ for SketchyCGAL, but its performance is slightly slower than
that for the setting $r=20$.

Table 2: Comparison of Algorithm1, SDPLR and SketchCGAL for quadratic knapsack SDP.  problem | algorithm | obj | Rp | Rd | pdgap | time
---|---|---|---|---|---|---
$n$ = 1000 | Algorithm1 | -1.1464123e+06 | 2.13e-15 | 4.61e-08 | 1.61e-10 | 1.39e+00
$\beta$ = 1.00e-01 | SDPLR | -1.1464103e+06 | 7.93e-07 | 1.03e-07 | 7.44e-06 | 2.37e+01
$p$ = 2.50e-01 | SketchyCGAL | -1.2127791e+06 | 7.41e-02 | 1.83e-07 | 1.81e-02 | 7.18e+02
$n$ = 1000 | Algorithm1 | -3.6163898e+06 | 7.07e-11 | 3.68e-08 | 6.05e-09 | 7.31e-01
$\beta$ = 3.00e-01 | SDPLR | -3.6163884e+06 | 9.96e-07 | 1.05e-06 | 2.23e-05 | 5.13e+01
$p$ = 2.50e-01 | SketchyCGAL | -7.3592186e+06 | 3.96e-01 | 4.70e-04 | 3.52e-01 | 3.60e+03
$n$ = 1000 | Algorithm1 | -6.2329083e+06 | 3.05e-15 | 3.47e-08 | 1.01e-09 | 4.96e-01
$\beta$ = 5.00e-01 | SDPLR | -6.2329001e+06 | 9.70e-07 | 8.31e-07 | 2.55e-07 | 7.55e+01
$p$ = 2.50e-01 | SketchyCGAL | -9.2273123e+06 | 2.06e-01 | 3.34e-02 | 5.88e-01 | 3.60e+03
$n$ = 1000 | Algorithm1 | -8.8340915e+06 | 2.14e-10 | 2.01e-08 | 1.19e-10 | 6.26e-01
$\beta$ = 7.00e-01 | SDPLR | -8.8341094e+06 | 9.89e-07 | 3.92e-07 | 4.59e-06 | 1.17e+02
$p$ = 2.50e-01 | SketchyCGAL | -1.2078637e+07 | 2.11e-01 | 4.75e-01 | 9.46e-01 | 3.60e+03
$n$ = 2000 | Algorithm1 | -4.8166859e+06 | 2.83e-10 | 2.57e-09 | 1.46e-09 | 4.30e+00
$\beta$ = 1.00e-01 | SDPLR | -4.8166801e+06 | 9.53e-07 | 6.07e-06 | 2.73e-06 | 1.24e+02
$p$ = 2.50e-01 | SketchyCGAL | - | - | - | - | -
$n$ = 2000 | Algorithm1 | -1.4947423e+07 | 4.66e-15 | 4.03e-09 | 1.83e-09 | 2.62e+00
$\beta$ = 3.00e-01 | SDPLR | -1.4947429e+07 | 9.97e-07 | 7.55e-07 | 9.12e-06 | 1.14e+03
$p$ = 2.50e-01 | SketchyCGAL | - | - | - | - | -
$n$ = 2000 | Algorithm1 | -2.5092783e+07 | 5.67e-10 | 1.03e-08 | 5.42e-11 | 2.53e+00
$\beta$ = 5.00e-01 | SDPLR | -2.5092737e+07 | 8.64e-07 | 9.21e-07 | 3.69e-06 | 1.29e+03
$p$ = 2.50e-01 | SketchyCGAL | - | - | - | - | -
$n$ = 2000 | Algorithm1 | -3.5316003e+07 | 3.98e-10 | 1.23e-09 | 1.45e-10 | 2.34e+00
$\beta$ = 7.00e-01 | SDPLR | -3.5316024e+07 | 9.76e-07 | 1.55e-07 | 1.85e-06 | 1.31e+03
$p$ = 2.50e-01 | SketchyCGAL | -5.0470000e+07 | 3.06e-01 | 5.87e-01 | 9.96e-01 | 3.60e+03
$n$ = 5000 | Algorithm1 | -3.1439664e+07 | 1.02e-11 | 1.07e-09 | 6.85e-09 | 1.28e+01
$\beta$ = 1.00e-01 | SDPLR | -3.1439583e+07 | 9.38e-07 | 4.14e-08 | 9.06e-06 | 1.68e+03
$p$ = 2.50e-01 | SketchyCGAL | - | - | - | - | -
$n$ = 5000 | Algorithm1 | -9.5441832e+07 | 9.24e-15 | 1.01e-08 | 2.13e-10 | 1.64e+01
$\beta$ = 3.00e-01 | SDPLR | - | - | - | - | -
$p$ = 2.50e-01 | SketchyCGAL | - | - | - | - | -
$n$ = 5000 | Algorithm1 | -1.5831412e+08 | 1.03e-14 | 8.48e-09 | 7.72e-11 | 1.88e+01
$\beta$ = 5.00e-01 | SDPLR | - | - | - | - | -
$p$ = 2.50e-01 | SketchyCGAL | - | - | - | - | -
$n$ = 5000 | Algorithm1 | -2.2086974e+08 | 1.12e-14 | 1.58e-08 | 1.05e-11 | 1.64e+01
$\beta$ = 7.00e-01 | SDPLR | - | - | - | - | -
$p$ = 2.50e-01 | SketchyCGAL | -3.1589361e+08 | 4.02e-01 | 6.44e-01 | 9.99e-01 | 3.60e+03
$n$ = 10000 | Algorithm1 | -1.2146659e+08 | 1.09e-10 | 1.47e-08 | 6.19e-08 | 5.78e+01
$\beta$ = 1.00e-01 | SDPLR | - | - | - | - | -
$p$ = 2.50e-01 | SketchyCGAL | - | - | - | - | -
$n$ = 10000 | Algorithm1 | -3.7203424e+08 | 1.51e-14 | 1.15e-08 | 6.97e-11 | 7.75e+01
$\beta$ = 3.00e-01 | SDPLR | - | - | - | - | -
$p$ = 2.50e-01 | SketchyCGAL | - | - | - | - | -
$n$ = 10000 | Algorithm1 | -6.2572380e+08 | 1.73e-14 | 1.15e-08 | 6.77e-12 | 6.70e+01
$\beta$ = 5.00e-01 | SDPLR | - | - | - | - | -
$p$ = 2.50e-01 | SketchyCGAL | - | - | - | - | -
$n$ = 10000 | Algorithm1 | -8.8036592e+08 | 1.87e-14 | 2.73e-09 | 6.84e-11 | 7.19e+01
$\beta$ = 7.00e-01 | SDPLR | - | - | - | - | -
$p$ = 2.50e-01 | SketchyCGAL | - | - | - | - | -

From the results in Table LABEL:SQKR, we can see that Algorithm 1 is much more
efficient than the other two algorithms in all the examples. Algorithms 1 can
solve all the instances accurately in a short time. In particular, it can
solve each of the SDP problems (SQKE) with $n=10,000$ very accurately in about
a minute. SDPLR can only solve the instances of moderate size. As for
SketchyCGAL, it reaches the maximum running time in all the instances except
for the first one. Also, the output of Algorithm 1 is always more accurate
than its competitors. Moreover, when $\beta$ is increased, the running time of
SDPLR increases dramatically, while the running time of Algorithm 1 becomes
even better. This means that Algorithm 1 is more robust than SDPLR. One
possible reason is that for SDPLR, the algorithm is based on ALM. For
infeasible methods like ALM and ADMM, their speed is very sensitive to the
penalty parameter. Since the best penalty parameter is related to the input
data, a different input data may result in significant difference in the
running time. However, Algorithm 1 is a feasible method which does not need
any penalty parameter updating. Thus, Algorithm 1 is more robust than SDPLR
and this also shows the advantage of using a feasible method.

###  6.3 Quadratic knapsack problem.

In this section, we solve (QKP) approximately by first solving (SQKELR) and
then followed by applying our rounding procedure described in subsection 5.3.
In order to test the accuracy of a feasible solution $d$ of (QKP), we compute
the relative gap between $\left\langle C,X\right\rangle$ for some optimal
solution $Y=\begin{pmatrix}1&x^{\top}\\\ x&X\end{pmatrix}$ of (SQKE) and
$d^{\top}Cd$ i.e.,

${\rm relgap}=\frac{|\left\langle
C,X\right\rangle-d^{\top}Cd|}{1+|d^{\top}Cd|}.$ (60)

We will compare our algorithm with a method called DP Heuristic proposed by
Fomeni and Letchford in [17], which is a modified dynamic programming
algorithm. The source codes can be downloaded from
https://sites.google.com/view/franklindjeumoufomeni/research/quadratic-
knapsack-problem. It has been shown in the numerical experiments of [17] that
their algorithm can find a feasible solution for (QKP) that is quite close to
the exact solution. We generate the dataset in the same way as mentioned in
the previous subsection. However, since their codes require that the knapsack
capacity $\tau$ to be an integer, we choose $\tau$ to be $\lceil\beta\cdot
e^{\top}a\rceil,$ where $\beta\in\\{0.1,0.5\\}.$ We also choose the density
parameter $p\in\\{0.1,0.5\\}.$ Although the theoretical rank bound for
(SQKELR) is $\lceil\sqrt{2(n+1)}\rceil+2$, we find that the rank of an optimal
solution of (SQKELR) is usually less than 20. Therefore, we choose the rank
parameter $r:=\min\left\\{20,\lceil\sqrt{2(n+1)}\rceil+2\right\\}$ to further
accelerate Algorithm 1. We set the maximum running time to be 3600 seconds.

Table 3: Comparison between Algorithm1 with rounding procedure and DP Heuristic. $n\ |$ $\beta\ |$ $p\ |$ $f_{\rm SDP}$ | algorithm | obj | relgap | time
---|---|---|---|---
$100|\ 0.1|\ 0.1$ | Algorithm1 | 5524 | 1.49e-01 | 1.25e-01
6.3497166e+03 | DP Heuristic | 5736 | 1.07e-01 | 1.88e-01
$100|\ 0.1|\ 0.5$ | Algorithm1 | 29283 | 1.45e-02 | 5.56e-02
2.9708357e+04 | DP Heuristic | 29216 | 1.69e-02 | 1.91e-01
$100|\ 0.5|\ 0.1$ | Algorithm1 | 28535 | 1.48e-02 | 3.67e-02
2.8957497e+04 | DP Heuristic | 28535 | 1.48e-02 | 3.94e-01
$100|\ 0.5|\ 0.5$ | Algorithm1 | 125252 | 1.06e-03 | 3.98e-02
1.2538509e+05 | DP Heuristic | 125252 | 1.06e-03 | 3.35e-01
$500|\ 0.1|\ 0.1$ | Algorithm1 | 120312 | 2.41e-02 | 5.09e-01
1.2321499e+05 | DP Heuristic | 120554 | 2.21e-02 | 7.45e+01
$500|\ 0.1|\ 0.5$ | Algorithm1 | 598541 | 2.38e-03 | 5.67e-01
5.9996349e+05 | DP Heuristic | 599355 | 1.02e-03 | 7.88e+01
$500|\ 0.5|\ 0.1$ | Algorithm1 | 648454 | 4.31e-03 | 3.16e-01
6.5124711e+05 | DP Heuristic | 648915 | 3.59e-03 | 9.35e+01
$500|\ 0.5|\ 0.5$ | Algorithm1 | 3165783 | 1.91e-04 | 4.47e-01
3.1663868e+06 | DP Heuristic | 3166162 | 7.10e-05 | 9.08e+01
$1000|\ 0.1|\ 0.1$ | Algorithm1 | 500226 | 9.65e-03 | 7.21e-01
5.0505314e+05 | DP Heuristic | 501407 | 7.27e-03 | 1.17e+03
$1000|\ 0.1|\ 0.5$ | Algorithm1 | 2414103 | 4.22e-03 | 5.01e-01
2.4243016e+06 | DP Heuristic | 2420049 | 1.76e-03 | 1.19e+03
$1000|\ 0.5|\ 0.1$ | Algorithm1 | 2574765 | 1.64e-03 | 4.91e-01
2.5789825e+06 | DP Heuristic | 2575821 | 1.23e-03 | 1.34e+03
$1000|\ 0.5|\ 0.5$ | Algorithm1 | 12428992 | 2.10e-03 | 7.54e-01
1.2455036e+07 | DP Heuristic | 12453855 | 9.48e-05 | 1.34e+03
$5000|\ 0.1|\ 0.1$ | Algorithm1 | 12644311 | 1.07e-03 | 6.92e+00
1.2657897e+07 | DP Heuristic | - | - | -
$5000|\ 0.1|\ 0.5$ | Algorithm1 | 60961427 | 1.77e-03 | 5.51e+00
6.1069318e+07 | DP Heuristic | - | - | -
$5000|\ 0.5|\ 0.1$ | Algorithm1 | 62978380 | 5.24e-04 | 6.12e+00
6.3011370e+07 | DP Heuristic | - | - | -
$5000|\ 0.5|\ 0.5$ | Algorithm1 | 315829856 | 1.99e-05 | 6.95e+00
3.1583613e+08 | DP Heuristic | - | - | -
$10000|\ 0.1|\ 0.1$ | Algorithm1 | 48436428 | 8.39e-04 | 2.44e+01
4.8477060e+07 | DP Heuristic | - | - | -
$10000|\ 0.1|\ 0.5$ | Algorithm1 | 243618278 | 8.80e-04 | 2.13e+01
2.4383258e+08 | DP Heuristic | - | - | -
$10000|\ 0.5|\ 0.1$ | Algorithm1 | 251900539 | 1.04e-04 | 2.68e+01
2.5192683e+08 | DP Heuristic | - | - | -
$10000|\ 0.5|\ 0.5$ | Algorithm1 | 1251839038 | 1.82e-04 | 2.68e+01
1.2520674e+09 | DP Heuristic | - | - | -

From the results in Table 1.1, we can see that the output of DP Heuristic is
slightly more accurate than the SDP model with rounding procedure. This is
reasonable because DP Heuristic is a dynamic programming algorithm with a lot
of delicated enhancements. However, the overall accuracy of the SDP model and
DP Heuristic is at the same level. Moreover, Algorithm 1 is much more
efficient than DP Heuristic. It has been mentioned in [17] that the
computational complexity of DP Heuristic is $O(n^{2}\tau).$ Since we choose
$\tau=\lceil\beta\cdot e^{\top}a\rceil=\Theta(n),$ the complexity becomes
$O(n^{3}),$ which becomes extremely slow when $n$ is large. This can also be
seen from Table 1.1. When $n=100,$ both DP Heuristic and Algorithm 1 can
terminate within a second. However, when $n=1000,$ DP Heuristic needs more
than 1000 seconds to return a solution while the running time of our algorithm
is still less than a second. In this case, our algorithm is more than 1000
times faster than DP Heuristic. When $n=5000$ or $10000,$ DP Heuristic reaches
the maximum running time but our algorithm can terminate within 30 seconds.

###  6.4 Testing Algorithm 1 for problems with non-regular points.

Although Algorithm 1 is able to verify the global optimality of a non-regular
point and escape from it if it is not optimal, we have not met any non-regular
points in the numerical experiments of the previous subsections. One reason is
that for a generic input data
$\left(a,\tau\right)\in\mathbb{R}^{n}\times\mathbb{R}$ that satisfies
Assumption 1, $\mathcal{K}_{n,r}$ is regular everywhere except for a trivial
point (see Corollary 3.2). Also, a trivial point implies that we do not select
any item in the knapsack problem, which is unlikely to happen in practice. In
order to test the ability of Algorithm 1 to deal with non-regular points, we
generate a special dataset
$\left(C^{\prime},a^{\prime},\tau^{\prime}\right)\in\mathbb{S}^{n}\times\mathbb{R}^{n}\times\mathbb{R}$
such that the optimal solution of (SQKE) contains a non-regular point, i.e.,
an integer solution. For an even number $n\in\mathbb{N}^{+},$ we first
generate a random data
$\left(C,a,\tau\right)\in\mathbb{S}^{n}\times\mathbb{R}^{n}\times\mathbb{R}$
using the same method as in subsection 6.2. We then modify the data
$\left(C,a,\tau\right)$ to get a new data
$\left(C^{\prime},a^{\prime},\tau^{\prime}\right)\in\mathbb{S}^{n}\times\mathbb{R}^{n}\times\mathbb{R}$
as follows: For any $i,j\in[n]$,

$C^{\prime}_{ij}:=\begin{cases}C_{ij}&{i\ {\rm and}\ j\ {\rm are\ even}}\\\
0&{\rm otherwise}\end{cases},\quad a^{\prime}_{i}:=a_{2\lceil
i/2\rceil},\quad\tau^{\prime}=\frac{1}{2}e^{\top}a^{\prime}.$ (61)

Consider problem (SQKE) and (SQKELR) with input data
$\left(C^{\prime},a^{\prime},\tau^{\prime}\right).$ Define two vectors
$v_{1},v_{2}\in\\{0,1\\}^{n}$ such that the index of nonzero entries of
$v_{1}$ and $v_{2}$ are odd numbers and even numbers respectively. From (61),
we have that $a^{\prime\top}v_{1}=a^{\prime\top}v_{2}=\tau^{\prime}.$ Thus,
both $v_{1}e_{1}^{\top}$ and $v_{2}e_{1}^{\top}$ are inside
$\mathcal{N}_{n,r}$ and so they are non-regular points of $\mathcal{K}_{n,r}.$
Also, it is easy to verify that $R=v_{2}e_{1}^{\top}$ is a global optimal
solution of (SQKELR) and $Y=\left(1\ v_{2}^{\top}\right)^{\top}\left(1\
v_{2}^{\top}\right)$ is a global optimal solution of (SQKE). Actually, if for
any $i=2k\in[n],$ there exists $j=2\ell\in[n]$ such that $C^{\prime}_{ij}>0,$
then $R=v_{2}e_{1}^{\top}$ and $Y=\left(1\ v_{2}^{\top}\right)^{\top}\left(1\
v_{2}^{\top}\right)$ are the only global optimal solution of (SQKELR) and
(SQKE) respectively. In this numerical experiment, we choose
$n\in\\{1000,2000,5000,7000,10000\\},$ and $p=0.25.$ We choose $r=3$ because
all the instances have a rank-one optimal solution. We generate the initial
point $R_{0}$ by adding a small perturbation to the non-regular non-optimal
point $v_{1}e_{1}^{\top}.$ The Matlab code for the initialization is: R0 =
v1*[1,0,0]+rand(n,3)/(1000*n); In this case, the escaping procedure of
Algorithm 1 is likely to be triggered at the beginning of the algorithm. To
demonstrate the usefulness of the rounding procedure in Algorithm 1, we also
tested Algorithm 1 without the rounding procedure, i.e., the Riemannian
gradient descent method with BB step and non-monotone line search mentioned in
subsection 5.2. We call this method RieBB.

Table 4: Comparison between Algorithm 1 and the Riemannian optimization method, RieBB without the rounding procedure, for quadratic knapsack SDPs with non-regular optimal solutions. problem | Algorithm | obj | Rp | Rd | pdgap | time
---|---|---|---|---|---|---
$n$ = 1000 | Algorithm1 | -3.1407370e+06 | 0.00e+00 | 0.00e+00 | 1.31e-16 | 3.42e-01
$p$ = 0.25 | RieBB | -3.1407370e+06 | 2.99e-10 | 0.00e+00 | 5.19e-12 | 4.78e+00
$n$ = 2000 | Algorithm1 | -1.2725874e+07 | 0.00e+00 | 0.00e+00 | 1.30e-16 | 7.57e-01
$p$ = 0.25 | RieBB | -1.2725874e+07 | 9.64e-11 | 0.00e+00 | 3.70e-14 | 1.03e+00
$n$ = 5000 | Algorithm1 | -7.8879089e+07 | 2.22e-16 | 0.00e+00 | 2.09e-16 | 4.75e+00
$p$ = 0.25 | RieBB | -7.8879089e+07 | 3.03e-10 | 3.08e-15 | 1.90e-12 | 6.43e+00
$n$ = 7000 | Algorithm1 | -1.5456007e+08 | 1.78e-15 | 3.86e-18 | 1.57e-15 | 9.95e+00
$p$ = 0.25 | RieBB | -1.5456007e+08 | 1.35e-12 | 0.00e+00 | 3.86e-16 | 1.12e+02
$n$ = 10000 | Algorithm1 | -3.1552221e+08 | 8.88e-15 | 0.00e+00 | 1.26e-15 | 2.00e+01
$p$ = 0.25 | RieBB | -3.1552221e+08 | 7.03e-12 | 0.00e+00 | 1.13e-15 | 4.51e+01

From the results in Table LABEL:SQKSIG, we can see that Algorithm 1 can solve
all of the instances efficiently. Also, note that the KKT residues of the
solutions are close to the rounding error. This is because our test problems
has an integer solution $v_{2}e_{1}^{\top},$ which can be obtained exactly by
the rounding step in Algorithm 1. The performance of RieBB is unstable. For
the instance with $n=7000$, it is more than 10 times slower than Algorithm 1.
This is because the Riemannian optimization method, RieBB, may suffer from
degeneracy issue when the iterations approach a non-regular point, and that
may slow down the algorithm. However, this will not happen in Algorithm 1
because the rounding procedure allows the algorithm to detect a non-regular
point in advance and escape from it if it is non-optimal.

###  6.5 Testing Algorithm 1 for large sparse problems.

It has been verified in the previous subsections that Algorithm 1 is more
efficient than other SDP solvers that can also utilise the low rank property
of a linear SDP problem. In this subsection, we move on to test some large
instances of (SQKE) where $n\geq 10000$ and $C$ is a sparse random matrix. The
matrix $C$ is generated in the same way as in subsection 6.2 with density
$p=\log(n)/n.$ We choose $\beta=0.5$ and $n\in\\{10^{4},2\cdot 10^{4},5\cdot
10^{4},10^{5},2\cdot 10^{5},5\cdot 10^{5},10^{6}\\}.$ We also set $r=20$ as in
subsection 6.3. Note that when $n$ is very large, computing the dual residue
Rd (see (6)) will be extremely slow because it requires the full eigenvalue
decomposition of a large dense matrix ${\bf S},$ which is the dual slack
defined in (41). However, due to the special structure of ${\bf S},$ i.e., it
can be written as the sum of a sparse matrix and a low rank matrix, we can use
the solver lobpcg developed by Knyazev in [27] to compute the smallest
eigenvalue of ${\bf S}$ efficiently. The source codes can be downloaded from
https://github.com/lobpcg/blopex. We redefine the dual residue ${\rm
Rd}:=\frac{\max\\{0,-\lambda_{\min}\left({\bf S}\right)\\}}{1+\|{\bf
S}\|_{F}}$, so that only the smallest eigenvalue of ${\bf S}$ needs to be
computed. It is easy to see that ${\rm Rd}=0$ if and only if ${\bf
S}\in\mathbb{S}^{n+1}_{+}.$

Table 5: Testing Algorithm 1 for quadratic knapsack SDPs with large sparse profit matrix $C.$ problem | obj | relgap | Rp | Rd | pdgap | time
---|---|---|---|---|---|---
$n$ = 10000 | -2.7522454e+06 | 1.86e-02 | 6.76e-14 | 5.97e-08 | 7.11e-08 | 7.29e+00
$n$ = 20000 | -5.8646881e+06 | 1.88e-02 | 4.32e-11 | 9.77e-08 | 5.59e-09 | 1.25e+01
$n$ = 50000 | -1.5887212e+07 | 1.74e-02 | 2.37e-12 | 5.98e-08 | 3.97e-07 | 2.68e+01
$n$ = 100000 | -3.3509978e+07 | 1.81e-02 | 3.72e-12 | 3.52e-08 | 2.81e-08 | 4.16e+01
$n$ = 200000 | -7.1116522e+07 | 2.41e-02 | 2.44e-11 | 1.54e-09 | 4.50e-08 | 7.98e+01
$n$ = 500000 | -1.8974901e+08 | 2.64e-02 | 7.80e-09 | 4.90e-09 | 1.57e-07 | 3.39e+02
$n$ = 1000000 | -3.9832359e+08 | 2.79e-02 | 2.47e-10 | 1.51e-09 | 2.40e-08 | 1.26e+03

From the results in Table LABEL:SQKLARGE, we can see that Algorithm 1 can
solve all the problems successfully. For the instance of size $n=10^{6},$
Algorithm 1 can compute a solution with high accuracy in about 20 minutes.
Also, it should be noted that for a sparse profit matrix $C,$ the relative gap
between the SDP bound and the objective value of the rounding solution is
close to $2\%,$ which is considered to be very good for an NP-hard integer
programming problem.

##  7 Conclusion.

In this paper, we study one of the SDP relaxations of the quadratic knapsack
problem. We consider the low rank formulation of this SDP problem to reduce
the problem’s dimensionality. In order to solve the low rank problem
efficiently, we explore the geometric properties of its feasible region
$\mathcal{K}_{n,r}$, which is an algebraic variety. We prove that this
algebraic variety is smooth everywhere for a generic input data except for a
trivial point. We also study the local geometric properties of
$\mathcal{K}_{n,r}$ at non-regular points on this algebraic variety. In order
to find the global optimal solution of this non-convex problem, we derive a
new-rank bound for these two problems to be equivalent. We also prove that
there is no spurious local minima under some rank condition without using any
regularity condition. All these good properties allow us to design a feasible
method based on Riemannian optimization, which can also handle non-regular
points. We conduct numerical experiments to verify the efficiency and
robustness of our feasible method. We also apply a rounding procedure to the
optimal solution of the SDP relaxation to get a feasible solution of the
quadratic knapsack problem. Numerical experiments show that our feasible
method with rounding strategy is much more efficient than another heuristic
based on dynamic programming. Traditional feasible methods focus on simple
smooth manifolds like oblique manifold and Stiefel manifold and thus their
application range is quite limited. However, this work demonstrates the
possibility of extending a Riemannian optimization framework to algebraic
varieties. It would be interesting to generalize this idea for solving more
general SDP problems with more complicated constraints. We leave this for
future work.

## References

  * [1] Absil, P. A. and Malick, J. (2012). Projection-like retractions on matrix manifolds. SIAM Journal on Optimization, 22(1), 135-158.
  * [2] Barvinok, A. I. (1995). Problems of distance geometry and convex properties of quadratic maps. Discrete & Computational Geometry, 13(2), 189-202.
  * [3] Barzilai, J. and Borwein, J. M. (1988). Two-point step size gradient methods. IMA Journal of Numerical Analysis, 8(1), 141-148.
  * [4] Bhojanapalli, S., Boumal, N., Jain, P. and Netrapalli, P. (2018, July). Smoothed analysis for low-rank solutions to semidefinite programs in quadratic penalty form. In Conference On Learning Theory (pp. 3243-3270). PMLR.
  * [5] Billionnet, A. and Soutif, É. (2004). An exact method based on Lagrangian decomposition for the 0–1 quadratic knapsack problem. European Journal of Operational Research, 157(3), 565-575.
  * [6] Bonnans, J. F. and Shapiro, A. (2013). Perturbation analysis of optimization problems. Springer Science & Business Media.
  * [7] Boumal, N. (2020). An introduction to optimization on smooth manifolds. Available online.
  * [8] Boumal, N., Mishra, B., Absil, P. A. and Sepulchre, R. (2014). Manopt, a Matlab toolbox for optimization on manifolds. Journal of Machine Learning Research, 15(1), 1455-1459.
  * [9] Boumal, N., Voroninski, V. and Bandeira, A. (2016). The non-convex Burer-Monteiro approach works on smooth semidefinite programs. Advances in Neural Information Processing Systems, 29.
  * [10] Boumal, N., Voroninski, V. and Bandeira, A. S. (2020). Deterministic guarantees for Burer‐Monteiro factorizations of smooth semidefinite programs. Communications on Pure and Applied Mathematics, 73(3), 581-608.
  * [11] Burer, S. and Monteiro, R. D. (2001). A projected gradient algorithm for solving the maxcut SDP relaxation. Optimization Methods and Software, 15(3-4), 175-200.
  * [12] Burer, S. and Monteiro, R. D. (2003). A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization. Mathematical Programming, 95(2), 329-357.
  * [13] Burer, S. and Monteiro, R. D. (2005). Local minima and convergence in low-rank semidefinite programming. Mathematical Programming, 103(3), 427-444.
  * [14] Caprara, A., Pisinger, D. and Toth, P. (1999). Exact solution of the quadratic knapsack problem. INFORMS Journal on Computing, 11(2), 125-137.
  * [15] Fallat, S. M. and Hogben, L. (2007). The minimum rank of symmetric matrices described by a graph: a survey. Linear Algebra and its Applications, 426(2-3), 558-582.
  * [16] Ferreira, C. E., Martin, A., de Souza, C. C., Weismantel, R. and Wolsey, L. A. (1996). Formulations and valid inequalities for the node capacitated graph partitioning problem. Mathematical Programming, 74(3), 247-266.
  * [17] Fomeni, F. D. and Letchford, A. N. (2014). A dynamic programming heuristic for the quadratic knapsack problem. INFORMS Journal on Computing, 26(1), 173-182.
  * [18] Friedland, S. and Loewy, R. (1976). Subspaces of symmetric matrices containing matrices with a multiple first eigenvalue. Pacific Journal of Mathematics, 62(2), 389-399.
  * [19] Gao, B. and Absil, P. A. (2022). A Riemannian rank-adaptive method for low-rank matrix completion. Computational Optimization and Applications, 81(1), 67-90.
  * [20] Gallo, G., Hammer, P. L. and Simeone, B. (1980). Quadratic knapsack problems. In Combinatorial optimization (pp. 132-149). Springer, Berlin, Heidelberg.
  * [21] Goemans, M. X. and Williamson, D. P. (1995). Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming. Journal of the ACM (JACM), 42(6), 1115-1145.
  * [22] Helmberg, C., Rendl, F. and Weismantel, R. (2000). A semidefinite programming approach to the quadratic knapsack problem. Journal of Combinatorial Optimization, 4(2), 197-215.
  * [23] Hogben, L. and Van Der Holst, H. (2007). Forbidden minors for the class of graphs G with $\xi(G)\leq 2$. Linear Algebra and its Applications, 423(1), 42-52.
  * [24] Hogben, L. (2010). Minimum rank problems. Linear Algebra and its Applications, 432(8), 1961-1974.
  * [25] Iannazzo, B. and Porcelli, M. (2018). The Riemannian Barzilai–Borwein method with nonmonotone line search and the matrix geometric mean computation. IMA Journal of Numerical Analysis, 38(1), 495-517.
  * [26] Johnson, E. L., Mehrotra, A. and Nemhauser, G. L. (1993). Min-cut clustering. Mathematical Programming, 62(1), 133-151.
  * [27] Knyazev, A. V. (2001). Toward the optimal preconditioned eigensolver: Locally optimal block preconditioned conjugate gradient method. SIAM Journal on Scientific Computing, 23(2), 517-541.
  * [28] https://www.mosek.com/.
  * [29] Pataki, G. (1998). On the rank of extreme matrices in semidefinite programs and the multiplicity of optimal eigenvalues. Mathematics of Operations Research, 23(2), 339-358.
  * [30] Petersen, K. B. and Pedersen, M. S. (2008). The matrix cookbook. Technical University of Denmark, 7(15), 510.
  * [31] Pisinger, D. (2005). Where are the hard knapsack problems? Computers & Operations Research, 32(9), 2271-2284.
  * [32] Pisinger, D. (2007). The quadratic knapsack problem—a survey. Discrete Applied Mathematics, 155(5), 623-648.
  * [33] Raydan, M. (1997). The Barzilai and Borwein gradient method for the large scale unconstrained minimization problem. SIAM Journal on Optimization, 7(1), 26-33.
  * [34] Tang, T. and Toh, K. C. (2021). Solving graph equipartition SDPs on an algebraic variety. arXiv preprint arXiv:2112.04256.
  * [35] Toh, K. C., Todd, M. J. and Tütüncü, R. H. (1999). SDPT3—a Matlab software package for semidefinite programming, version 1.3. Optimization Methods and Software, 11(1-4), 545-581.
  * [36] Tütüncü, R. H., Toh, K. C. and Todd, M. J. (2003). Solving semidefinite-quadratic-linear programs using SDPT3. Mathematical Programming, 95(2), 189-217.
  * [37] Wang, H., Kochenberger, G. and Glover, F. (2012). A computational study on the quadratic knapsack problem with multiple constraints. Computers & Operations Research, 39(1), 3-11.
  * [38] Wen, Z. and Yin, W. (2013). A feasible method for optimization with orthogonality constraints. Mathematical Programming, 142(1), 397-434.
  * [39] Wright, S. and Nocedal, J. (1999). Numerical optimization. Springer Science, 35(67-68), 7.
  * [40] Yang, L., Sun, D. and Toh, K. C. (2015). SDPNAL+: a majorized semismooth Newton-CG augmented Lagrangian method for semidefinite programming with nonnegative constraints. Mathematical Programming Computation, 7(3), 331-366.
  * [41] Yurtsever, A., Tropp, J. A., Fercoq, O., Udell, M. and Cevher, V. (2021). Scalable semidefinite programming. SIAM Journal on Mathematics of Data Science, 3(1), 171-200.
  * [42] Zhang, R. (2020). Newton retraction as approximate geodesics on submanifolds. arXiv preprint arXiv:2006.14751.
  * [43] Zhao, X. Y., Sun, D. and Toh, K. C. (2010). A Newton-CG augmented Lagrangian method for semidefinite programming. SIAM Journal on Optimization, 20(4), 1737-1765.

##  Appendix A Useful auxiliary results

###  A.1 A corollary of the implicit function theorem

######  Lemma A.1.

Let $F:\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{k}$ be a
continuously differentiable mapping in a neighbourhood of
$\left(x_{0},y_{0}\right)\in\mathbb{R}^{n}\times\mathbb{R}^{m}.$ Suppose
$F(x_{0},y_{0})=0$ and ${\rm
D}_{y}F(x_{0},y_{0}):\mathbb{R}^{m}\rightarrow\mathbb{R}^{k}$ is a surjective
linear mapping. Then there exists $\delta>0$ and a continuously differentiable
mapping $y:B_{\delta}(x_{0})\rightarrow\mathbb{R}^{k}$ such that
$y(x_{0})=y_{0}$ and $F(x,y(x))=0.$ Here,
$B_{\delta}(x_{0}):=\left\\{x\in\mathbb{R}^{n}:\ {\rm
dist}(x,x_{0})<\delta\right\\}.$

###### Proof.

Proof. Since ${\rm D}_{y}F(x_{0},y_{0})$ is surjective, we have that $m\geq
k.$ If $m=k,$ then Lemma A.1 directly follows from the implicit function
theorem and the mapping $y$ is unique. If $m>k$, we may find a linear mapping
$\mathcal{A}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m-k}$ such that the linear
operator $\triangle y\rightarrow\left({\rm D}_{y}F(x_{0},y_{0})[\triangle
y],\mathcal{A}(\triangle y)\right)$ is a bijective mapping from
$\mathbb{R}^{m}$ to $\mathbb{R}^{m}.$ Thus, we can define a new mapping
$\widehat{F}:\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ such
that $\widehat{F}(x,y):=\left(F(x,y),\mathcal{A}(y)-\mathcal{A}(y_{0})\right)$
and apply the implicit function theorem to $\widehat{F}.$

###  A.2 (29) and (30) are strictly feasible

######  Lemma A.2.

Suppose $(a,\tau)$ satisfies Assumption 1. Then both (29) and (30) have a
strictly feasible solution.

###### Proof.

Proof. It is easy to see that the feasibility of (30) is strictly satisfied by
choosing $\beta$ to be sufficiently negative. We only have to find a strictly
feasible solution for (29).

Case 1. $v\neq 0.$
In this case, $\sigma_{v}=1.$ Let $X_{1}:=\frac{ee^{\top}}{n}.$ We have that
$X_{1}\succeq 0,$ $\mathrm{Tr}(X_{1})=1.$ Also,

$\left\langle aa^{\top}-\sigma_{v}\tau\cdot{\rm diag}\left(a\circ
d\right),X_{1}\right\rangle=\frac{\left(e^{\top}a\right)^{2}-\tau
d^{\top}a}{n}>\frac{\left(e^{\top}a\right)^{2}-\tau e^{\top}a}{n}>0,$ (62)

where the last inequality comes from the fact that $e^{\top}a>\tau.$ Because
$v\neq 0,$ there exists $i\in[n]$ such that $d_{i}=1.$ For $t\in[0,1],$ define
$Y_{t}=(1-t){\rm diag}\left(e_{i}\right)+\frac{t}{n}I.$ We have that
$Y_{t}\succeq 0$ and $\mathrm{Tr}(Y_{t})=1.$ Also,

$\left\langle aa^{\top}-\sigma_{v}\tau\cdot{\rm diag}\left(a\circ
d\right),Y_{0}\right\rangle=a_{i}^{2}-\tau a_{i}<0,$ (63)

where the last inequality comes from Assumption 1. From continuity, there
exists $t\in(0,1)$ such that $\left\langle aa^{\top}-\sigma_{v}\tau\cdot{\rm
diag}\left(a\circ d\right),Y_{t}\right\rangle<0$ and $Y_{t}\succ 0.$ From
(62), there exists $\beta\in\left(0,1\right)$ such that $\beta
X_{1}+(1-\beta)Y_{t}\succ 0,$ $\left\langle I,\beta
X_{1}+(1-\beta)Y_{t}\right\rangle=1$ and

$\left\langle aa^{\top}-\sigma_{v}\tau\cdot{\rm diag}\left(a\circ
d\right),\beta X_{1}+(1-\beta)Y_{t}\right\rangle=0.$ (64)

Therefore, $\beta X_{1}+(1-\beta)Y_{t}$ is a strict feasible solution of (29).

Case 2. $v=0.$
In this case, $\sigma_{v}=-1$ and $aa^{\top}-\sigma_{v}\tau\cdot{\rm
diag}\left(a\circ d\right)=aa^{\top}-\tau{\rm diag}\left(a\right).$ We can
construct a strictly feasible solution in the same way as in Case 1.