\section{Problems and Motivation}
Bilevel optimization is a special kind of optimization where an (upper-level) optimization problem contains another (lower-level) optimization
problem as a constraint.
The first bilevel optimization problem  was  realized in the field of game theory by a German economist Heinrich Freiherr von Stackelberg in 1934 \cite{stack}.
Since then, these problems  have been found  in many real-world applications,
  including  artificial intelligence,   cryptography,  decision science, defensive force structure design,   game theory,     topology  design  of structures and  communication networks,   etc (see \cite{bc,cml,kjy,wz,yhl}).
It was discovered recently by Gao \cite{gao-to17}  that the well-studied  topology  optimization  in computational mechanics should be  formulated as
a   bilevel knapsack problem (BKP), i.e. the upper-level optimization is a knapsack problem, while the lower-level optimization is governed by the well-known minimum total potential energy principle.
The knapsack problem is one of the most fundamental problems in combinatorial optimization. It has been
 studied extensively  for more than a century  in multidisciplinary  fields (see \cite{k-p-p}).
The quadratic knapsack problem (QKP)
  is an extension of knapsack problem that allows for quadratic terms in the
target  function. QKP was  first introduced in 19th century (see  \cite{wit}),
 which has a wide range of applications including telecommunication, transportation network, computer science and economics. In fact, Witzgall first discussed QKP when selecting sites for satellite stations in order to maximize the global traffic with respect to a budget constraint. Similar model applies to problems like considering the location of airports, railway stations, or freight handling terminals \cite{rhy}.
Applications of QKP in the field of computer science is more common after the early days: compiler design problem \cite{h-r-w}, clique problem \cite{p-l-p},    very large scale integration (VLSI) design and manufacturing  \cite{f-m-s}.
In bilevel optimization of  multi-scale complex systems,
both linear and quadratic   knapsack problems can appear in either upper or lower optimization.
Indeed, over the last 20 years,  a variety of bilevel knapsack problems has been proposed (see \cite{bhm,cclw}).

 Generally speaking, bilevel knapsack    problems are  extremely difficult from
the computational point of view and cannot be expressed in terms of classical integer
programs (which can only handle a single level of optimization). Due to the integer constraint and bilevel coupling, bilevel knapsack  problems are  non-differentiable and   non-convex. It has been proven in \cite{vsj}
 that merely evaluating a solution for optimality is also a NP-hard task.
In fact,  it was shown in \cite{hanJS}  that even the linear bilevel problems are  strongly NP-hard.  The proof for the non-existence of a polynomial
time algorithm for linear bilevel problems can be found in \cite{deng}.
Classical  theories and methods  can't be used directly for solving  this type of problems.
Therefore, the general bilevel knapsack problem
 could be one of the most challenging problems in computer science and global optimization.

During the past 40 years  many approximation methods have been proposed for numerically  solving bilevel
optimization  problems. Impressive collections of these methods and applications can be found in the books
 \cite{bard,dem1,z-l-g} as well as the review article \cite{dem2}.

Evolutionary method \cite{smfb,wjl} is a   popular approach, where the lower level optimization
problem is solved corresponding to each and every upper
level member. One of the first evolutionary algorithms for solving bilevel
optimization problems was proposed in the early 1990s.
Mathieu et al. \cite{mpa}  used a nested approach with genetic
algorithm at the upper level, and linear programming at the
lower level. One of the early works on discrete bilevel optimization was
by Vicente et al. (1996) \cite{vsj}, which focused on discrete linear
bilevel programs, and analyzed the properties and existence
of the optimal solution for different kinds of discretization
arising from the upper and lower level variables. The authors
have shown in the paper that certain compactness conditions
guarantee the existence of optimal solution in continuous
linear bilevel programs, discrete-continuous linear
bilevel programs and discrete-discrete linear bilevel programs.
The conditions are equivalent to stating that the inducible
region is non-empty. However, the existence conditions in
the case of continuous-discrete linear bilevel programs are
not straightforward.
In the classical literature, branch-and-bound and branch-and-cut are some of the commonly used deterministic methods to handle discreteness in variables.
But  these well-studied  methods can be computational expensive, and can be used only for solving very small-sized problems. Therefore, global
optimization problems with 200 variables are referred to as ``medium scale", problems with 1000 variables as
``large scale", and the so-called ``extra-large scale" is only around 4000 variables. However, any simple problem
in topology  design can easily have  millions of variables \cite{gao-to18}.
Although  there are some fully polynomial time approximation schemes for discrete optimization,
they are not always applicable in practice due to memory requirements \cite{whs}.
Also these algorithms are based on  heuristic techniques,
 the obtained solutions come with no guarantee of global optimality and may get stuck in local minima.
Therefore,  it  is truly important to develop a powerful deterministic method for solving large-scale  general bilevel  knapsack problems.

Canonical duality theory  (CDT) is a  precise  methodological theory, which can be used not only for modeling complex systems within a unified framework, but also for solving a large class of challenging problems in nonconvex analysis and global optimization \cite{gao-book00}.
This theory is particularly powerful for solving integer constrained problems.
 In 2007  Gao discovered  that  by simply using a standard  canonical measure $\Lambda(\mathbf{z}) = \mathbf{z}^2 - \mathbf{z}$, the 0-1 constrained problems can be equivalently converted to a unified
continuous concave maximization problem in a convex  feasible space, which can be solved deterministically via well-developed convex optimization techniques \cite{gao-jimo07}.  This method has been generalized for solving general discrete optimization problems \cite{gao-cace,gao-watsonetal,ruan-gao-minl,wang-etal}.
 Applications of the canonical duality theory to multidisciplinary study
 was given recently \cite{g-l-r-17} and an analytic solution to linear knapsack problem has been obtained \cite{gao-to17,gao-to18}.

This paper deals with a general  bilevel knapsack problem (BKP), in which,  the upper-level optimization is a quadratic knapsack problem coupled with the continuous follower variables, while the lower-level optimization is a mixed integer nonlinear minimization problem involves both leader and follower variables
(since an alternative iteration method is used in this paper, this order of upper and lower levels  can reversed).
The main contributions of this paper include:
\begin{verse}
1. Canonical dual solutions to quadratic   knapsack problems;\\
2. Improved analytic solutions to  knapsack problems and  criteria for NP-hardness;\\
3. A volume reduction method combined with an alternative iteration for solving general BKP;\\
4. A powerful  CDT algorithm for solving general BKP  with an application to topology design.
\end{verse}
The   rest of the paper is organized as follows.
 The next section provides the mathematical formulation of the  bi-level knapsack problem and its special linear case.
 A decoupled alternative iteration is suggested.
 The   theoretical results of this paper are presented in  Section 3, including a canonical penalty-duality method,
 a general  analytic  solution form, as well as the existence and uniqueness of this solution to the quadratic knapsack problem.
In Section \ref{sec4} , the NP-hardness of the knapsack problems is addressed.
Improved analytic solutions are provided for both quadratic and linear knapsack problems.
 Section 5 proposes a volume reduction method and a CPD algorithm for solving  BKP.
 In Section 6, applications are  illustrated by  both 2D and 3D benchmark problems
 in   topology design. The performance and novelty of the CPD algorithm are verified by computational
results. Finally, Section 7 presents conclusions and open problems  for future research.

\section{Bilevel Optimization Modeling and  Knapsack Problems}
The bilevel optimization  problem proposed to solve is given below:
\begin{eqnarray}
 &({\cal P}):    &   \min_{  \mathbf{x} \in {\cal X}_a ,   {\bf z}} %{\mbox{\boldmath$\rho$} \in {\cal Z}_a} \left\{P(\mathbf{x},  \mathbf{z} ) =   \frac{1}{2}  {\bf z}} %{\mbox{\boldmath$\rho$}^T \mathbf{Q}( \mathbf{x}) \mathbf{z} - \mathbf{z}^T {\bf c}(\mathbf{x})  \right \} \;\;\;\; \label{eq-upp}\\
&& s.t.    \;\; \mathbf{x} \in \arg \min_{  \mbox{\boldmath$\chi$} \in {\cal X}_a}  \left\{ \Pi( \mbox{\boldmath$\chi$},  \mathbf{z})  = G(\mathbf{D} \mbox{\boldmath$\chi$}, \mathbf{z}) - \mbox{\boldmath$\chi$}^T \mathbf{f}  \right \} \;\;\;\;\;\; \;\;\;\; \label{eq-lowp}
\end{eqnarray}
  where $P: {\cal X}_a\times {\cal Z}_a \subset \mathbb{R}^m\times \mathbb{R}^n \rightarrow \mathbb{R} $ is  the upper-level target
 function,
  $\Pi :{\cal X}_a \times {\cal Z}_a   \rightarrow \mathbb{R} $ is  the lower-level target
 function, $\mathbf{Q}(\mathbf{x}) $ and ${\bf c}(\mathbf{x})$ are two pre-given matrix-valued and vector-valued  functions of $\mathbf{x} \in {\cal X}_a$, respectively; $\mathbf{f} \in \mathbb{R}^m$ is a given vector;
  $\mathbf{D}:{\cal X}_a \rightarrow {\cal G}_a \subset  \mathbb{R}^p$ is a linear operator, and $G:{\cal G}_a\times {\cal Z}_a \rightarrow \mathbb{R}$ is a given  function, which is usually a nonconvex objective function
   of $\mathbf{g} = \mathbf{D} \mathbf{x}$, but linearly depends on $\mathbf{z} \in {\cal Z}_a$.
  The upper-level (or leader) variable $\mathbf{z}$ is a discrete vector, whose domain is a subset of Boolean space (i.e. a knapsack):
 \begin{equation}
{\cal Z}_a = \{ \mathbf{z} = \{ z_i\} \in  \{ 0, 1 \}^n| \;\; \mathbf{v}^T \mathbf{z} \le V_c \} ,
\end{equation}
where $V_c > 0$ is a given  size of the knapsack, $\mathbf{v} = \{ v_i\}  \in \mathbb{R}^n_+$ is a given positive  vector.
The lower-level (or the follower) variable $\mathbf{x}$ is a continuous vector, whose domain ${\cal X}_a$ is a convex set of $\mathbb{R}^m$, in which, only linear equality constraints are included.

  Clearly,  the problem $({\cal P})$  is a  bilevel knapsack problem.
   Due to the integer constraint, the nonlinearity/nonconvexity of $\Pi(\mathbf{x}, \mathbf{z})$, and the strong coupling between the upper and  lower level  problems,
    the proposed
    problem $({\cal P})$ could be one of the most  challenging problems in global optimization and computer science.
    It is fundamentally difficult to solve this problem directly.
    Therefore,  a   decoupled alternative iteration (DAI) method  will be adopted in this paper, i.e., the  coupled bilevel optimization
    is   split into two decoupled subproblems by alternative iteration:
      \begin{verse}
  (i) For a given initial   ${\bf z}} %{\mbox{\boldmath$\rho$}_{k-1} \in {\cal Z}_a$, to solve the lower-level problem (\ref{eq-lowp})  first for
  \begin{equation}
  \mathbf{x}_k =  \arg \min \{ \Pi( \mathbf{x}, {\bf z}} %{\mbox{\boldmath$\rho$}_{k-1})   \;\; | \;\; \mathbf{x}\in {\cal X}_a \} . \label{eq-lowk}
  \end{equation}

  (ii) For the fixed  $\mathbf{x}_k \in {\cal X}_a$, to solve the upper-level knapsack problem (\ref{eq-upp}) for
  \begin{equation}
\mathbf{z}_k = \arg   \min \{  P(   \mathbf{x}_k, \mathbf{z})    | \; \;\; \mathbf{z}\in {\cal Z}_a \} \label{eq-uppk}
  \end{equation}
  \end{verse}

The upper-level minimization \eqref{eq-uppk}  is the well-known quadratic knapsack problem (QKP).
The  binary variable $\mathbf{z} = \{ z_i\}$  represents whether item $i$ is included in the knapsack ${\cal Z}_a$.
The  given vector  $\mathbf{c}(\mathbf{x}) = \{c_i(\mathbf{x}) \} \in \mathbb{R}^n_+$ must be
 positive  for any given $\mathbf{x} \in {\cal X}_a$ such that $c_i$
  is the profit earned by selecting item $i$.
 The given matrix $\mathbf{Q}(\mathbf{x}) \in \mathbb{R}^{n\times n} $ must be symmetrical,  its  diagonal elements are usually
$Q_{ii} = 0, \;\; i=1, \dots, n$  and $- Q_{ij} $ is the profit achieved if both item $i$ and $j$ are added.

  An important special case is  the following linear knapsack problem (LKP):
 \begin{equation}
 ({\cal P}_{\ell}): \;\; \min \{ P_{\ell} (\mathbf{z}) = - {\bf c}^T \mathbf{z} | \;\; \mathbf{v}^T \mathbf{z} \le V_c,  \;\; \mathbf{z} \in \{ 0, 1\}^n\}.
 \end{equation}
  A 1998 study of the Stony Brook University Algorithm Repository showed that, out of 75 algorithmic problems, the knapsack problem was the 19th most popular and the 3rd most needed after suffix trees and the bin packing problem.
 It is well-known that the    knapsack problem is NP-hard and there is no known algorithm that can solve the problem in polynomial time.
Actually, even the  linear knapsack problem is listed as one of Karp's 21 NP-complete problems~\cite{karp}.
  In this paper, the  canonical  duality method will be addressed for solving this challenging problem.

For a given design variable $\mathbf{z}_{k-1}$, the lower-level minimization  \eqref{eq-lowk}
 is the  general nonlinear/nonconvex  optimization problem proposed by the author for multi-scale systems \cite{gao-aip,gao-amma18}.

\begin{remark}[Objectivity,  Modeling, and Canonical Duality]
Objectivity is a central concept in our daily life, related to
reality and truth. In science, the objectivity is often attributed to the property of scientific
measurement, as the accuracy of a measurement can be tested independent from the individual scientist who first
reports it. In continuum physics, a real-valued function $G:{\cal G}_a \rightarrow \mathbb{R}$ is called objective if and only if it is an invariant under orthonormal transformation (see Definition 6.1.2 \cite{gao-book00}), i.e.
\begin{equation}
G({\bf R } \mathbf{g}) = G(\mathbf{g}) \;\; \forall \mathbf{g} \in {\cal G}_a , \;\; \forall   {\bf R } \in SO(p) ,
\end{equation}
where $SO(p) = \{ {\bf R } \in \mathbb{R}^{p\times p}  | \;\;  {\bf R }^T = {\bf R }^T, \; \det {\bf R } =1 \}$ is a special orthogonal  group in $\mathbb{R}^p$.
Physically speaking, an objective function is governed by the intrinsic physical law of the system,
which does not depend on observers. Therefore, the objectivity is essential for any real-world
mathematical models. In continuum physics, it  is also called the principle of material
frame indifference \cite{oden}.  Geometrically, an objective function does not depend on rigid rotation of the
system considered, but only on certain measure of its variable. In the Euclidean space ${\cal G}_a \subset \mathbb{R}^p$, the simplest objective function is the $\ell_2$-norm
$\| \mathbf{g} \|$ since $\| {\bf R } \mathbf{g}\|^2 = \mathbf{g}^T {\bf R }^T {\bf R } \mathbf{g} = \| \mathbf{g} \|^2 \;\; \forall {\bf R } \in SO(p)$. Therefore, an objective function must be nonlinear.

Correspondingly, the linear term $F(\mathbf{x}) = \mathbf{x}^T \mathbf{f}  $ in the lower-level problem \eqref{eq-lowp}  is called a {\em subjective function} \cite{gao-aip,gao-to18}, where the given input $\mathbf{f}$ and the constraints in ${\cal X}_a$   depend only   on each given problem.
Thus, if the operator ``$-$" is considered as the predicate, the
 difference  $\Pi(\mathbf{x}) = G(D \mathbf{x}) - F(\mathbf{x})$ between  object and subject
 forms a complete and grammatically correct mathematical formulation.
  In continuum physics, if $G(\mathbf{g})$ is the free (or internal) energy, $F(\mathbf{x})$ is the input (or external)
 energy, then $\Pi(\mathbf{x})$ is the total potential energy and the minimum total potential energy principle leads to
a general variational/optimization problem in mathematical physics \cite{gao-strang89}.
By the fact that ${\cal G}_a \subset \mathbb{R}^p$ and ${\cal X}_a \subset \mathbb{R}^m$ can be  in different dimensional spaces with different measures,
  the lower-level minimization
$\min \Pi(\mathbf{x}) $ presented in this paper covers general constrained nonconvex/nonsmooth  optimization
problems in multi-scale complex systems \cite{gao-book00,gao-amma18,g-l-r-17}.

According to  \cite{ciarlet},  $G(\mathbf{g})$ is an objective function if and only if
there exists an objective measure $\mbox{\boldmath$\epsilon$} = \mathbf{g}^T \mathbf{g}$ and a real-valued function $\Phi(\mbox{\boldmath$\epsilon$})$ such that $G(\mathbf{g}) =\Phi(\mathbf{g}^T \mathbf{g})$.
In continuum physics  and differential geometry, if $\mathbf{g} = D \mathbf{x}$ is the deformation gradient, then the objective measure $\mbox{\boldmath$\epsilon$} = \mathbf{g}^T \mathbf{g}$ is the well-known right Cauchy-Green tensor. By the fact that the  free energy $\Phi(\mbox{\boldmath$\epsilon$})$ is usually convex (say the St Venant-Kirchhoff material \cite{gao-book00}),
 the duality relation $\mbox{\boldmath$\epsilon$}^* = \nabla \Phi(\mbox{\boldmath$\epsilon$})$ is invertible. This one-to-one constitutive relation is called the canonical duality and $\Phi(\mbox{\boldmath$\epsilon$})$ is called the canonical function. These basic truths in continuum physics laid a foundation for the
canonical duality theory, in which, the key idea  of the
{\em canonical transformation} is to choose a nonlinear operator (not necessary to be objective)
$\mbox{\boldmath$\epsilon$}= \Lambda(\mathbf{x})$ and a canonical function $\Phi(\mbox{\boldmath$\epsilon$})$ such that
the  given nonconvex function $G(D\mathbf{x})$ can be written in the canonical form $G(D\mathbf{x}) = \Phi(\Lambda(\mathbf{x}))$.
By the fact that the  objectivity plays a fundamental role in mathematical modeling,   the canonical
duality theory can be powerfully used for solving
many  challenging problems in multidisciplinary fields \cite{g-l-r-17}.

 However, this  important concept of the objectivity  has been extensively  misused in optimization literature such that the general problem in nonlinear optimization (or programming)
has been  formulated as
\begin{equation}
\min f(x) \;\; s.t. \;\;  g(x) = 0, \;\; h(x) \le 0 ,
\end{equation}
where $f(x)$ is called the ``objective function"\footnote{This terminology is used mainly in English literature. The function $f(x)$
 is correctly called the target function in all Chinese and Japanese literature, or the goal
function in some Russian and German literature by serious researchers.}, which is allowed to be any arbitrarily given function, even a linear
function. Clearly, this mathematical problem is too abstract. Although it enables one to model a very wide range of
problems, it comes at a price: many global optimization problems are considered to be NP-hard. Without detailed
information on these arbitrarily given functions, it is impossible to have a powerful theory for solving the artificially
given constrained problem. Also, due to this conceptual mistake,   the canonical
duality theory has been mistakenly challenged
in several publications\footnote{See  M.D. Voisei, C. Zalinescu, Some remarks concerning Gao-Strang’s complementary gap function, {\em Appl. Anal.} 90 (2011) 1111-1121.}.

The conceptual
mistake  in these false challenges revealed a big gap   between physics and optimization.
Interested readers are recommended to read the recent papers \cite{gao-opl16,g-r-l-17} for details.  \hfill $\diamondsuit$
\end{remark}

 The canonical duality theory for solving the lower-level  nonconvex continuous optimization problem (\ref{eq-lowk})
has been studied extensively during the past 15 years    (see \cite{gao-jogo00, gao-cace,g-l-r-17}).
 As long as the nonconvex function $G(\mathbf{D} \mathbf{x}, \mathbf{z})$ in $\Pi(\mathbf{x}, \mathbf{z})$  can be written in the
canonical form $\Phi(\Lambda(\mathbf{x}), \mathbf{z})$, this problem can be  solved easily by the canonical duality theory
  to obtain both global and local optimal solutions (see many real-world applications in \cite{g-l-r-17}).
This paper will show how to use the canonical duality theory for solving the challenging upper-level knapsack problem \eqref{eq-uppk}.

 \section{Canonical  Duality  Solution to Upper-Level Knapsack Problem}
\label{Pure}
 For a fixed $\mathbf{x} \in {\cal X}_a$,  the upper-level  optimization (\ref{eq-uppk})  can be written in  the typical  quadratic knapsack  problem:
  \begin{equation}
  ({\cal P}_q): \;\; \min_{  {\bf z}} %{\mbox{\boldmath$\rho$}\in \{ 0,1\}^n } \left\{ P_q({\bf z}} %{\mbox{\boldmath$\rho$}) = \frac{1}{2}    {\bf z}} %{\mbox{\boldmath$\rho$}^T \mathbf{Q} \mathbf{z}  - {\bf c}^T {\bf z}} %{\mbox{\boldmath$\rho$}  \;\; |\;\; \mathbf{v}^T {\bf z}} %{\mbox{\boldmath$\rho$} \le V_c  \right\}.
  \end{equation}

    According to the canonical duality theory for mathematical modeling \cite{gao-opl16,gao-to17}, the inequality
    ${\bf z}} %{\mbox{\boldmath$\rho$} \cdot \mathbf{v} \le V_c$ is a geometrical constraint, while the  integer constraint ${\bf z}} %{\mbox{\boldmath$\rho$} \in \{ 0,1\}^n$ in   $({\cal P}_q)$
 is a constitutive   condition \cite{gao-to18,gao-amma18}, which can be equivalently replaced by  the canonical constraint
 $\mathbf{z}^2 - \mathbf{z} = {\bf 0} \in \mathbb{R}^n \;\;\forall \mathbf{z} \in \mathbb{R}^n$, where $\mathbf{z}^2 = \mathbf{z}\circ \mathbf{z}  = \{ z_i^2\}  \in \mathbb{R}^n$.
 Therefore, by introducing a penalty parameter $\beta > 0$ and let  ${\cal Z}_b = \{    \mathbf{z} \in \mathbb{R}^n  \; | \;\; \mathbf{v}^T \mathbf{z} \le V_c  \}$,  the problem $({\cal P}_q)$ can be relaxed by  the canonical  penalty function:
\begin{equation}
({\cal P}_\beta): \;\;
\min_{ \mathbf{z} \in {\cal Z}_b} \left\{ P_\beta(\mathbf{z}) = \frac{1}{2} \mathbf{z}^T \mathbf{Q} \mathbf{z} +  \frac{1}{2} \beta \|\mathbf{z}^2 - \mathbf{z} \|^2  -  {\bf c}^T  {\bf z}} %{\mbox{\boldmath$\rho$}
\right\} \label{eq-pda} .
\end{equation}
Clearly, we have
\begin{equation}
\{ P_q(\mathbf{z})  | \;\; \mathbf{z} \in {\cal Z}_a\} = \lim_{\beta \rightarrow \infty} \{   P_\beta({\bf z}} %{\mbox{\boldmath$\rho$} )   | \;  \mathbf{z} \in {\cal Z}_b  \}.
\end{equation}
Although the integer constraint is relaxed by the   penalty function,
 \begin{equation}
 \Psi_\beta({\bf z}} %{\mbox{\boldmath$\rho$}) =  \frac{1}{2} \beta \|\mathbf{z}^2 - \mathbf{z} \|^2   ,
 \end{equation}
  the    problem (\ref{eq-pda}) is a nonconvex minimization in $\mathbb{R}^n$, which
 is still considered to be NP-hard   by traditional methods. This is the reason why the traditional external penalty method can be used mainly for linear constrained problems \cite{l-g-opl}.
  To solve this problem by  the canonical duality theory,  we need first to introduce the following canonical transformation:
  \[
  \Psi_\beta({\bf z}} %{\mbox{\boldmath$\rho$}) = \Phi_\beta(\Lambda({\bf z}} %{\mbox{\boldmath$\rho$})),
\]
\[
   \Phi_\beta(\bveps) =   \frac{1}{2}  \beta \| \bveps\|^2, \;\;
 \bveps = \Lambda({\bf z}} %{\mbox{\boldmath$\rho$}) = {\bf z}} %{\mbox{\boldmath$\rho$}^2 - \mathbf{z}  \in    \mathbb{R}^n.
  \]
Clearly,  $\Phi_\beta: \mathbb{R}^n  \rightarrow \mathbb{R}$ is a convex quadratic function. Its canonical dual can be simply given by the Legendre transformation
\begin{equation}
\Phi^*_\beta(\mbox{\boldmath$\sigma$}) = \max_{\bveps \in {\cal E} } \{ \bveps^T \mbox{\boldmath$\sigma$} - \Phi_\beta(\bveps)   \} =
\frac{1}{2} \beta^{-1}   \|\mbox{\boldmath$\sigma$}\|^2 .
\end{equation}
 Thus, replacing $\Psi_\beta({\bf z}} %{\mbox{\boldmath$\rho$}) $ in $P_\beta({\bf z}} %{\mbox{\boldmath$\rho$} )$  by the Fenchel-Young equality $\Phi_\beta(\Lambda({\bf z}} %{\mbox{\boldmath$\rho$}) = \Lambda({\bf z}} %{\mbox{\boldmath$\rho$})^T \mbox{\boldmath$\sigma$} - \Phi^*_\beta(\mbox{\boldmath$\sigma$})$,
 and introducing a Lagrange multiplier $\tau \ge 0$,
the canonical penalty function $P_\beta(\mathbf{z})$ can be reformed as  the   Gao-Strang total complementary function
  \begin{equation}
   \Xi_\beta({\bf z}} %{\mbox{\boldmath$\rho$}, \mbox{\boldmath$\sigma$},\tau) =  \frac{1}{2}   {\bf z}} %{\mbox{\boldmath$\rho$}^T \mathbf{G}(\mbox{\boldmath$\sigma$} ) \mathbf{z}  -  \frac{1}{2}  \beta^{-1} \| \mbox{\boldmath$\sigma$}\|^2
  -  {\bf z}} %{\mbox{\boldmath$\rho$}^T \mbox{\boldmath$\psi$} (\mbox{\boldmath$\sigma$}, \tau  )  -  \tau V_c ,
\end{equation}
  where
\begin{equation}
\mathbf{G} (\mbox{\boldmath$\sigma$}) = \mathbf{Q} + 2 \mbox{Diag}(\mbox{\boldmath$\sigma$}), \;\; \mbox{\boldmath$\psi$} (\mbox{\boldmath$\sigma$}, \tau  ) = {\bf c}  - \tau \mathbf{v}   + \mbox{\boldmath$\sigma$}.
\end{equation}
 By introducing a  canonical dual feasible  space
\begin{equation}
 {\cal S}_a^+ = \{ \mbox{\boldmath$\zeta$} = ( \mbox{\boldmath$\sigma$}, \tau) \in \mathbb{R}^{n+1} | \;\; \mathbf{G}(\mbox{\boldmath$\sigma$})  \succ 0 , \;\; \tau >  0 \},
\end{equation}
 the function  $  \Xi_\beta({\bf z}} %{\mbox{\boldmath$\rho$}, \mbox{\boldmath$\zeta$}) $  is convex in ${\bf z}} %{\mbox{\boldmath$\rho$} \in \mathbb{R}^n$ for   any given $\mbox{\boldmath$\zeta$} \in {\cal S}^+_a$.
The canonical penalty-duality function $ P^d_\beta : {\cal S}^+_a \rightarrow \mathbb{R}$  can be defined by
\begin{eqnarray*}
P_\beta^d(\mbox{\boldmath$\zeta$}) &=& \min_{{\bf z}} %{\mbox{\boldmath$\rho$} \in \mathbb{R}^n} \{  \Xi_\beta ({\bf z}} %{\mbox{\boldmath$\rho$}, \mbox{\boldmath$\zeta$}) |\;\; \mbox{\boldmath$\zeta$} \in {\cal S}^+_a \} \\
&=&  - \frac{1}{2}  \mbox{\boldmath$\psi$}(\mbox{\boldmath$\zeta$})^T \mathbf{G}(\mbox{\boldmath$\sigma$}) ^{-1} \mbox{\boldmath$\psi$}(\mbox{\boldmath$\zeta$}) -
\frac{1}{2} \beta^{-1}  \| \mbox{\boldmath$\sigma$}  \|^2 -\tau V_c .
\end{eqnarray*}
Thus, the  canonical penalty-duality  problem  can be proposed in the following
 \begin{equation}
 ({\cal P}^d_\beta): \;\;\;\; \max \left\{ P^d_\beta (\mbox{\boldmath$\zeta$}) | \;\;\mbox{\boldmath$\zeta$} \in {\cal S}^+_a \right\}.
 \end{equation}
  \begin{thm}[Complementary-Dual Principle]\label{thm1}
 For any given $\beta  > 0$, if $  (\mathbf{z}_\beta,  \mbox{\boldmath$\zeta$}_\beta)    \in  \mathbb{R}^n \times   {\cal S}^+_a$ is a KKT point of $\Xi_\beta(\mathbf{z}, \mbox{\boldmath$\zeta$})$,     then
 $\mathbf{z}_\beta$  is a global minimum solution to  the canonical  penalty   problem $({\cal P}_\beta)$,
 $\mbox{\boldmath$\zeta$}_\beta$ is a  solution to $({\cal P}^d_\beta)$, and
  \begin{equation}\label{eq-cdp}
P_\beta(\mathbf{z}_\beta) = \min_{\mathbf{z} \in {\cal Z}_b} P_\beta(\mathbf{z}) =    \Xi_\beta (\bar{\brho}, \bar{\bzeta} )=
\max_{\mbox{\boldmath$\zeta$} \in {\cal S}^+_a} P^d_\beta(\mbox{\boldmath$\zeta$}) =    P_\beta^d(\mbox{\boldmath$\zeta$}_\beta).
\end{equation}
\end{thm}
 {\bf Proof}. By the fact that $\Xi_\beta :\mathbb{R}^n \times {\cal S}^+_a \rightarrow \mathbb{R}$ is a saddle function,  if $  (\mathbf{z}_\beta,  \mbox{\boldmath$\zeta$}_\beta)    \in  \mathbb{R}^n \times   {\cal S}^+_a$ is a KKT point of $\Xi_\beta(\mathbf{z}, \mbox{\boldmath$\zeta$})$,   it must be a saddle point of $\Xi_\beta(\mathbf{z}, \mbox{\boldmath$\zeta$})$.
 Then by the definition of $P_\beta^d(\mbox{\boldmath$\zeta$})$, we have
 \begin{eqnarray}
\Xi_\beta(\mathbf{z}_\beta, \mbox{\boldmath$\zeta$}_\beta) &=&  \min_{\mathbf{z} \in \mathbb{R}^n} \max_{\mbox{\boldmath$\zeta$} \in {\cal S}^+_a} \Xi(\mathbf{z}, \mbox{\boldmath$\zeta$}) =   \max_{\mbox{\boldmath$\zeta$} \in {\cal S}^+_a}  \min_{\mathbf{z} \in \mathbb{R}^n}\Xi(\mathbf{z}, \mbox{\boldmath$\zeta$}) \nonumber  \\
&  =&  \max_{\mbox{\boldmath$\zeta$} \in {\cal S}^+_a}  P^d_\beta(\mbox{\boldmath$\zeta$}) = P^d_\beta(\mbox{\boldmath$\zeta$}_\beta).
 \end{eqnarray}
   By the KKT conditions
   \[
\mathbf{v}^T \mathbf{z}_\beta  \le V_c, \;\; \tau_\beta \ge 0 , \;\;
 \tau_\beta(  \mathbf{v}^T \mathbf{z}_\beta -  V_c) = 0 ,
 \]
and the condition $\mbox{\boldmath$\zeta$}_\beta = (\mbox{\boldmath$\sigma$}_\beta, \tau_\beta) \in {\cal S}^+_a$,
we have  $\tau_\beta > 0$ and  $\mathbf{v}^T \mathbf{z}_\beta -  V_c = 0$.
By the
   convexity of the penalty function $\Psi_\beta (\mbox{\boldmath$\epsilon$})$, we have
$ \Psi^{**}_\beta (\mbox{\boldmath$\epsilon$}) = \Psi_\beta(\mbox{\boldmath$\epsilon$}) $.
 Thus, for any given $\beta > 0$
\begin{equation}
P_\beta(\mathbf{z}_\beta) = \min_{\mathbf{z} \in {\cal Z}_b} P_\beta ({\bf z}} %{\mbox{\boldmath$\rho$} ) = \min_{\mathbf{z} \in{\cal Z}_b} \max_{\mbox{\boldmath$\zeta$} \in {\cal S}^+_a}  \Xi_\beta ({\bf z}} %{\mbox{\boldmath$\rho$}, \mbox{\boldmath$\zeta$})   = \Xi_\beta (\mathbf{z}_\beta, \mbox{\boldmath$\zeta$}_\beta) .
\end{equation}
Thus, $\mathbf{z}_\beta$ is a global optimal solution to $({\cal P}_\beta)$.  \hfill $\Box$\\

Theorem 1 shows that the nonconvex minimization problem $({\cal P}_\beta)$ is canonically dual to a concave maximization problem  $({\cal P}^d_\beta)$ in
a convex space ${\cal S}^+_a$, which can be  solved easily by well-developed convex minimization methods.

 \begin{remark}[Canonical Penalty-Duality    and $\beta$-Perturbation]
By the facts that
\begin{equation}\label{eq-pane}
  \Phi(\bveps) = \lim_{\beta \rightarrow \infty}  \Phi_\beta(\bveps)    = \left\{ \begin{array}{ll}
0 & \mbox{ if } \bveps  = {\bf 0} \in \mathbb{R}^n\\
+ \infty & \mbox{ otherwise} \end{array} \right.
 \end{equation}
and
\begin{equation}\label{eq-pane*}
\Phi^*(\mbox{\boldmath$\sigma$}) = \sup_{\bveps \in \mathbb{R}^n}  \{ \bveps^T \mbox{\boldmath$\sigma$} - \Phi(\bveps) \} =
0  \;\; \forall \mbox{\boldmath$\sigma$} \in \mathbb{R}^n,
 \end{equation}
we have
\begin{eqnarray*}
P^d_q(\mbox{\boldmath$\zeta$})  & = & \lim_{\beta \rightarrow + \infty} P_\beta^d (\mbox{\boldmath$\zeta$}) \\
&=&  - \frac{1}{2}  \mbox{\boldmath$\psi$}(\mbox{\boldmath$\zeta$})^T \mathbf{G}(\mbox{\boldmath$\sigma$})^{-1} \mbox{\boldmath$\psi$}(\mbox{\boldmath$\zeta$})
  -\tau V_c  - \Phi^*(\mbox{\boldmath$\sigma$}),
  \end{eqnarray*}
  which is exactly the canonical dual function to the primal function $P_q(\mathbf{z})$.
  Therefore, $({\cal P}^d_\beta)$ is actually the so-called $\beta$-perturbation of the canonical dual problem $({\cal P}^d_q)$:
  \begin{equation}
  ({\cal P}^d_q): \;\; \max \{ P^d_q (\mbox{\boldmath$\zeta$}) | \;\; \mbox{\boldmath$\zeta$} \in {\cal S}^+_a \} .
  \end{equation}
  The penalty-duality method  was first proposed  by  Gao for solving convex variational problems \cite{gao-thesis,gao-cs88}.
 The $\beta$-perturbation method for nonconvex  integer constrained problems was  first proposed
 in \cite{gao-to17,gao-ruan-jogo10}.
 It was proved by   Theorem 7  in \cite{gao-ruan-jogo10} that there exists a $\beta_c>0$ such that for any given $\beta \ge \beta_c$,
both  the $\beta$-perturbed canonical dual problem
$({\cal P}^d_\beta)$
 and the problem $({\cal P}^d_q)$ have  the same solution set.
 This shows the relation between the canonical penalty-duality method proposed in this paper and the $\beta$-perturbation method proposed in \cite{gao-ruan-jogo10}. \hfill$\diamondsuit$
 \end{remark}

  \begin{thm}[Analytic  Solution to Knapsack Problem]\label{thm2}
 For any given  $\beta  > 0$, if  $   \mbox{\boldmath$\zeta$}_\beta =(\mbox{\boldmath$\sigma$}_\beta, \tau_\beta)  \in   {\cal S}^+_a$ is a   solution to $({\cal P}^d_\beta)$,
   then  the  vector defined by
 \begin{equation} \label{eq-solu}
 \mathbf{z}_\beta   =   \mathbf{G}(\mbox{\boldmath$\sigma$}_\beta)^{-1} \mbox{\boldmath$\psi$} (\mbox{\boldmath$\zeta$}_\beta)
 \end{equation}
is a global minimum solution to the canonical penalty   problem $({\cal P}_\beta)$.

Moreover, there exists   $\beta_c  \gg 0$ such that $\beta\ge \beta_c$ and $\mathbf{z}_\beta \in {\cal Z}_a$, then $\mathbf{z}_\beta $ is a global optimal solution to the knapsack problem $({\cal P}_q)$
and
 \begin{equation}\label{eq-cdp}
P_q(\mathbf{z}_\beta) = \min_{\mathbf{z} \in {\cal Z}_b} P_\beta(\mathbf{z}) =    \Xi_\beta (\bar{\brho}, \bar{\bzeta} )=
\max_{\mbox{\boldmath$\zeta$} \in {\cal S}^+_a} P^d_\beta(\mbox{\boldmath$\zeta$})  =   P_q^d(\mbox{\boldmath$\zeta$}_\beta).
\end{equation}
\end{thm}
 {\bf Proof}. It is easy to verify  that if  $(\mathbf{z}_\beta, \mbox{\boldmath$\zeta$}_\beta)$ is  a KKT point of $\Xi_\beta(\mathbf{z}, \mbox{\boldmath$\zeta$})$, the
 criticality condition $\nabla_\mathbf{z} \Xi_\beta(\mathbf{z}_\beta, \mbox{\boldmath$\zeta$}_\beta) = 0 $ leads to  $ \mathbf{z}_\beta   =   \mathbf{G}(\mbox{\boldmath$\sigma$}_\beta)^{-1} \mbox{\boldmath$\psi$} (\mbox{\boldmath$\zeta$}_\beta) $.
 By Theorem \ref{thm1} we know that if  $   \mbox{\boldmath$\zeta$}_\beta =(\mbox{\boldmath$\sigma$}_\beta, \tau_\beta)  \in   {\cal S}^+_a$ is a   solution to $({\cal P}^d_\beta)$, the vector defined by (\ref{eq-solu})
 must be a solution to $({\cal P}_\beta)$.
Since the penalty function $\Phi_\beta(\mathbf{z}) \ge 0 \;\; \forall \mathbf{z} \in \mathbb{R}^n$, there must exists a sufficiently big $\beta_c \gg  0$
such that $\Phi_\beta(\mathbf{z}_\beta)  = 0 \; \forall \beta \ge \beta_c$, i.e.   $\mathbf{z}_\beta \in {\cal Z}_a$.
Thus, $P_q(\mbox{\boldmath$\zeta$}_\beta) = P_\beta(\mathbf{z}_\beta) = \min \{ P_q(\mathbf{z}) | \; \mathbf{z} \in {\cal Z}_a\}$. \hfill $\Box$\\

 Since ${\cal S}^+_a$ is a convex open set, its boundary is
 \begin{equation}
 \partial {\cal S}^+_a = \{ \mbox{\boldmath$\zeta$} = (\mbox{\boldmath$\sigma$}, \tau)  \in \mathbb{R}^{n+1} | \;\;  \det  \mathbf{G}(\mbox{\boldmath$\sigma$}) = 0 , \;\; \tau = 0 \}.
 \end{equation}
 The following theorem is important for understanding the NP-hardness of the knapsack problem.
  \begin{thm}[Existence and Uniqueness for $({\cal P}_q)$]\label{thm3}
 For any given ${\bf c}  \in \mathbb{R}^n$, $ \mathbf{v} \in \mathbb{R}^n_+ $, and $\mathbf{Q} \in \mathbb{R}^{n\times n}$,   if ${\cal S}^+_a$ has a  non-empty interior $\mbox{\boldmath$\zeta$}_o \in {\cal S}^+_a$ such that
 $P_q^d(\mbox{\boldmath$\zeta$}_o) >  - \infty $ and
 \begin{equation}
 \lim_{\mbox{\boldmath$\zeta$} \rightarrow \partial {\cal S}^+_a } P^d_q(\mbox{\boldmath$\zeta$}) = - \infty \;\; \;\; \forall \mbox{\boldmath$\zeta$} \in {\cal S}^+_a ,\label{eq-corc}
 \end{equation}
 then  the  problem $({\cal P}_q^d)$ has a unique solution $\mbox{\boldmath$\zeta$}_c \in {\cal S}^+_a$ and $\mathbf{z}_c= \mathbf{G}(\mbox{\boldmath$\sigma$}_c) ^{-1}\mbox{\boldmath$\psi$}(\mbox{\boldmath$\zeta$}_c) $ is a unique global optimal solution to  $({\cal P}_q)$.
 \end{thm}
 {\bf Proof}. This theorem is a direct application of the canonical duality  theory (see Theorem 4 in  \cite{gao-cace}\footnote{The condition $P^d(\mbox{\boldmath$\sigma$}_o) < -\infty$ in  \cite{gao-cace}  is a typo, it should be  $P^d(\mbox{\boldmath$\sigma$}_o) > -\infty$.}).
  By the facts that $P^d_q(\mbox{\boldmath$\zeta$})$ is concave and its feasible set ${\cal S}^+_a$ is convex and nonempty,   if the condition (\ref{eq-corc})
 holds, then   the canonical penalty-dual function $P^d_q:{\cal S}^+_a \rightarrow \mathbb{R} $ is strictly concave, bounded up and
$-P^d_q(\mbox{\boldmath$\zeta$})$ is coercive on ${\cal S}^+_a$. Therefore, $({\cal P}^d_q)$ has a unique solution $\mbox{\boldmath$\zeta$}_c$. By Theorem \ref{thm2}, we know that the associated  $\mathbf{z}_c= \mathbf{G}(\mbox{\boldmath$\sigma$}_c) ^{-1}\mbox{\boldmath$\psi$}(\mbox{\boldmath$\zeta$}_c) $ is a unique global optimal solution to  $({\cal P}_q)$. \hfill $\Box$ \\

 As an important application, let us re-call  the   linear knapsack problem:
 \begin{equation}
 ({\cal P}_{\ell}): \;\; \min_{\mathbf{z} \in {\cal Z}_a}  \{ P_{\ell} (\mathbf{z}) = - {\bf c}^T \mathbf{z}  \} .
 \end{equation}
The associated canonical penalty problem is
\begin{equation}
({\cal P}_{\ell\beta}): \;\;  \min_{  \mathbf{z} \in {\cal Z}_b}  \left\{ P_{\ell\beta} (\mathbf{z}) = \frac{1}{2} \beta \|\mathbf{z}^2 - \mathbf{z} \|^2   - {\bf c}^T \mathbf{z} \;\;
  \right\} .
\end{equation}
 Since the matrix $\mathbf{Q} $ vanished in this case, the  canonical penalty-dual function  has the following simple form:
 \begin{equation}
  P^d_{\ell\beta} (\mbox{\boldmath$\zeta$}) = -  \frac{1}{2}  \sum_{i=1}^n   \left( \frac{1}{2}  \sigma^{-1}_i (\sigma_i + c_i - \tau v_i)^2 + \beta^{-1} \sigma_i^2 \right)   - \tau V_c  .
 \end{equation}
Let
$ \mbox{\boldmath$\sigma$}^{-1} = {  \mbox{Diag}}(\mbox{\boldmath$\sigma$})^{-1}  =\{ \sigma^{-1}_i\}$ and
\[
{\cal S}^+_a = \{ \mbox{\boldmath$\zeta$} = (\mbox{\boldmath$\sigma$}, \tau) \in \mathbb{R}^{n+ 1} | \;\; \mbox{\boldmath$\sigma$} > {\bf 0} , \;\; \tau > 0 \}.
\]
The   canonical penalty-dual problem for the linear knapsack problem is
\begin{equation}
({\cal P}_{\ell\beta}^d): \;\; \max \{   P^d_{\ell\beta} (\mbox{\boldmath$\zeta$}) | \;\; \mbox{\boldmath$\zeta$} \in {\cal S}^+_a\}.
\end{equation}

 \begin{Corollary}[Analytic  Solution to Knapsack Problem \cite{gao-to18}]\label{thm4}
 For any given ${\bf c}, \mathbf{v} \in \mathbb{R}^n_+$ and $V_c , \beta > 0$,
 if    $   \mbox{\boldmath$\zeta$}_\beta =(\mbox{\boldmath$\sigma$}_\beta, \tau_\beta)  \in   {\cal S}^+_a$ is a   solution to $({\cal P}^d_{\ell\beta})$
 and $\tau_\beta \mathbf{v} - {\bf c} \neq {\bf 0}$,
   then  the penalty   problem $({\cal P}_{\ell\beta})$ has a unique global minimum solution which is given analytically
   by
 \begin{equation} \label{eq-solul}
 \mathbf{z}_\beta   =  \frac{1}{2}  \mbox{\boldmath$\sigma$}_\beta^{-1} \circ   (\mbox{\boldmath$\sigma$}_\beta + {\bf c} - \tau_\beta \mathbf{v}) .
 \end{equation}
Moreover, there exists a $\beta_c \gg  0$ such that $\beta\ge \beta_c$ and  $\mathbf{z}_\beta \in {\cal Z}_a$, then $\mathbf{z}_\beta $ is a global optimal solution to the linear knapsack problem $({\cal P}_\ell)$.
\end{Corollary}

This special  result for linear knapsack problem was obtained recently in \cite{gao-to18}.  For $\beta = 0$, the canonical dual problem $ ({\cal P}_{\ell}^d)$ of $({\cal P}_\ell)$ is
 \begin{equation}
\max_{\mbox{\boldmath$\zeta$} \in {\cal S}^+_a } \left\{ P^d_{\ell } (\mbox{\boldmath$\zeta$}) = -  \frac{1}{2}  \sum_{i=1}^n   \left( \frac{1}{2}  \sigma^{-1}_i (\sigma_i + c_i - \tau v_i)^2   \right)   - \tau V_c \right\} .
 \end{equation}
It is easy to prove that if $\tau_\beta \mathbf{v} - {\bf c} \neq {\bf 0}$, the condition (\ref{eq-corc})  is reduced by
 \begin{equation}
 \lim_{\mbox{\boldmath$\sigma$} \rightarrow {\bf 0}^+ } P^d_{\ell} (\mbox{\boldmath$\sigma$}, \tau_\beta) = -\infty  .
 \end{equation}
 Thus, by Theorem \ref{thm3} we know that the vector (\ref{eq-solul}) is a unique solution to  the problem $({\cal P}_{\ell\beta})$ and it is a solution to the  linear knapsack problem
 $({\cal P}_\ell)$ if $\beta \ge \beta_c \gg 0$ such that $\mathbf{z}_\beta \in {\cal Z}_a$.

Actually, for any given $\beta > 0$, the criticality condition $\nabla P^d_{\ell\beta} = 0$ leads to the following algebraic equations:
\begin{equation}
  \beta^{-1} \sigma_i^3 + \sigma_i^2 = (\tau v_i - c_i)^2, \;\; i = 1, \dots, n, \label{eq-cdas}
\end{equation}
\begin{equation}
\sum_{i=1}^n \frac{1}{2} \frac{v_i}{\sigma_i} ( \sigma_i- v_i \tau + c_i ) - V_c = 0 .\label{eq-cdv}
\end{equation}
 It was  proved in \cite{gao-book00,gao-jimo07} that for any given $\beta > 0$,  $\tau \ge 0$ and ${\bf c} \in \mathbb{R}^n$ such that
 $\theta_i = \tau v_i - c_i \neq 0, \forall  i = 1, \dots, n$, the canonical dual algebraic equation (\ref{eq-cdas}) has a unique
positive real solution
\begin{equation}
\sigma_i  =  \frac{1}{6} \beta   [- 1 +  \phi_i(\tau ) + \phi_i^c(\tau )] > 0 , \;\; i = 1, \dots, n
\label{eq-solus}
\end{equation}
where
\[
\phi_i(\varsigma )  = \eta^{-1/3} \left [2 \theta_i^2 - \eta + 2 \theta_i  \sqrt{ \theta_i^2 -\eta } \right]^{1/3} ,
 \;\; \eta = \frac{4 \beta^2}{27},
\]
and $\phi_i^c $ is the complex conjugate of $\phi_i $, i.e. $\phi_i  \phi_i^c  = 1$.

On the other hand, for a given $\mbox{\boldmath$\sigma$} \in \mathbb{R}^n_+$, the Lagrange multiplier $\tau$ can be uniquely obtained by
\begin{equation}
        \tau = \frac{ \sum_{i=1}^n  v_i (1 + c_i  /\sigma_i) - 2 V_c }{\sum_{i=1}^n v_i^2/  \sigma_i}. \label{eq-solust}
\end{equation}
 It is easy to prove that for any given
$\beta > 0$ and $\tau_{k-1} > 0$,
the solution $\mbox{\boldmath$\sigma$}_k $ produced by \eqref{eq-solus} satisfies
$\mbox{\boldmath$\sigma$}_k < {\bf c}$.
Thus, for a given initial $\tau_0 > 0$, an alternative iteration can be used for solving \eqref{eq-solus} and \eqref{eq-solust} and the
   sequence  $(\mbox{\boldmath$\sigma$}_k, \tau_k) \in {\cal S}^+_a$ approaches to the global optimal solution of $({\cal P}^d_\ell)$ in polynomial time.

Theorem \ref{thm4} shows that although the canonical dual problem is a concave maximization in continuous space,  it produces the analytical solution  (\ref{eq-solu})
to  the well-known integer Knapsack problem $({\cal P}_q)$. This truth was first discovered  by Gao in 2007
for  general quadratic integer programming problems (see Theorem 3, \cite{gao-jimo07}).

\section{Improved Solutions and NP-Hardness  } \label{sec4}
 By the fact that $\mbox{\boldmath$\alpha$} \circ \mathbf{z}^2 = \mbox{\boldmath$\alpha$}\circ  \mathbf{z} \;\; \forall \mathbf{z} \in \{0,1\}^n, \;  \forall \mbox{\boldmath$\alpha$} \in \mathbb{R}^n$,
for any given symmetrical  $\mathbf{Q} \in \mathbb{R}^{n\times n}$ we can choose an $\mbox{\boldmath$\alpha$}$ such that   $\mathbf{Q}_\alpha= \mathbf{Q} +2  \mbox{Diag} (\mbox{\boldmath$\alpha$}) \succeq 0$. Thus, by $\mathbf{c}_\alpha = \mathbf{c} +  \mbox{\boldmath$\alpha$} $,
the problem $({\cal P}_q)$ can be equivalently written in the  so-called $\alpha$-perturbation form  \cite{gao-cace}:
\[
({\cal P}_\alpha): \;\;    \min_{ {\bf z}} %{\mbox{\boldmath$\rho$}\in \{ 0,1\}^n}
\left\{ P_{\alpha} ({\bf z}} %{\mbox{\boldmath$\rho$}) = \frac{1}{2}    {\bf z}} %{\mbox{\boldmath$\rho$}^T \mathbf{Q}_\alpha  \mathbf{z}  - {\bf c}_\alpha^T {\bf z}} %{\mbox{\boldmath$\rho$}  \;\; |\;\; \mathbf{v}^T {\bf z}} %{\mbox{\boldmath$\rho$} \le V_c   \right\} .
  \]
Since  $\mbox{rank} \mathbf{Q}_\alpha = r \le n$,   there must exist  (see \cite{strang}) a
 ${\bf L }  \in \mathbb{R}^{r\times n} $ and
${\bf H} \in \mathbb{R}^{r\times r}$ with $\mbox{rank} {\bf L } = \mbox{rank} {\bf H} = r$ and ${\bf H} \succ  0$
 such that $\mathbf{Q}_\alpha = 4 {\bf L }^T {\bf H}{\bf L }$.
Similar to  the $\alpha$-perturbed canonical dual problem $({\cal P}^g_{ip})$ given in \cite{gao-cace},
 the  canonical dual function in $({\cal P}^d_q)$ can be
reformulated as \cite{gao-amma18}:
\[
  P^g_{\alpha}   (\mbox{\boldmath$\sigma$},\tau) = - \frac{1}{2}  \mbox{Abs}[ \mbox{\boldmath$\phi$}(\mbox{\boldmath$\sigma$}, \tau)]
    - \frac{1}{2} \mbox{\boldmath$\sigma$}^T {\bf H}^{-1} \mbox{\boldmath$\sigma$} - \tau V_b  + d  ,
\]
where $V_b = V_c - \frac{1}{2} \sum_{i=1}^n v_i , \;\; d = \frac{1}{8}  \sum_{i=1}^n  (2 \alpha_i+ \sum_{j=1}^n Q_{ij})  - \frac{1}{2} \sum_{i=1}^n (c_i  + \alpha_i) $ are two constants, and
\begin{equation}
   \mbox{\boldmath$\phi$} ( \mbox{\boldmath$\zeta$} ) =     \mathbf{c}  -   \tau \mathbf{v} -  2  {\bf L }^T \mbox{\boldmath$\sigma$}  - \frac{1}{2}  \mathbf{Q} {\bf e}  , \;\; {\bf e} = \{ 1\}^n .
\end{equation}
The notation $\mbox{Abs}[ \mbox{\boldmath$\phi$}(\mbox{\boldmath$\sigma$},\tau)]$ denotes
$ \mbox{Abs}[ \mbox{\boldmath$\phi$}(\mbox{\boldmath$\sigma$},\tau) ]= \sum_{i=1}^n | \phi_i(\mbox{\boldmath$\sigma$},\tau) | $.
Let
\[
{\cal S}^+_c= \{ \mbox{\boldmath$\zeta$} =(\mbox{\boldmath$\sigma$}, \tau) \in \mathbb{R}^{r+1} | \;\;  \;\; \tau >  0 \}.
\]
Then the improved canonical dual problem to $({\cal P}_\alpha)$ can be proposed as
\begin{equation}
({\cal P}^g_{\alpha}): \;\;\; \max \left\{ P^g_{\alpha}   (\mbox{\boldmath$\sigma$},\tau) | \;\; \mbox{\boldmath$\zeta$} \in {\cal S}^+_c \right\}.
\end{equation}
In many real-world applications, we have $r\ll n$, thus the problem $({\cal P}^g_{\alpha})$ is much easier than $({\cal P}^d_q)$. Similar to Theorem 7 in \cite{gao-cace}, we have the following improved result.
 \begin{thm}[Improved Analytic  Solution to $({\cal P}_q)$]\label{thm-eqkp}
 For any given $V_c > 0$,   $\mathbf{v}, \mathbf{c}  \in \mathbb{R}^n_+ $, $\mathbf{Q} \in \mathbb{R}^{n\times n}$ and an $\mbox{\boldmath$\alpha$} \in \mathbb{R}^n $ such that
  $\mathbf{Q}_\alpha = \mathbf{Q} + 2 \mbox{Diag} (\mbox{\boldmath$\alpha$})  = 4 {\bf L }^T {\bf H}{\bf L }$,  and
${\bf H} \succ  0$, if
 $  \mbox{\boldmath$\zeta$}_c =\{ \mbox{\boldmath$\sigma$}_c,   \tau_c\} $ is a  solution to $({\cal P}^g_{\alpha})$ and
\begin{equation}
\phi_i( \mbox{\boldmath$\zeta$}_c)   \neq { 0} \;\; \forall i \in \{ 1, \dots, n \} ,\label{eq-exis}
\end{equation}
then the quadratic knapsack problem $({\cal P}_{q})$ has a unique  global optimal solution
\begin{equation}\label{eq-solu2}
\mathbf{z}_c = \frac{1}{2} \left\{ \frac{  \phi_i (  \mbox{\boldmath$\zeta$}_c) }{| \phi_i (  \mbox{\boldmath$\zeta$}_c) |} + 1 \right\}^n
\end{equation}
and
 \begin{equation}
 P_{q} (\mathbf{z}_c) = \min_{{\bf z}} %{\mbox{\boldmath$\rho$} \in {\cal Z}_a} P_{q} ({\bf z}} %{\mbox{\boldmath$\rho$}) = \max_{\mbox{\boldmath$\zeta$} \in {\cal S}^+_c}  P^g_{\alpha} (\mbox{\boldmath$\zeta$})
=   P^g_{\alpha}(\mbox{\boldmath$\zeta$}_c).
\end{equation}
  Otherwise,  if $\phi_i( \mbox{\boldmath$\zeta$}_c) = 0 $ for at least one $i \in \{1, \dots, n\}$, then
$({\cal P}_{q})$  has at least two  solutions.
\end{thm}

It is easy to prove that $P^g_\alpha(\mbox{\boldmath$\zeta$})$ is strictly concave on the convex feasible set
${\cal S}_c^+$. If $ \mbox{\boldmath$\zeta$}_c \in {\cal S}^+_c$ is a solution to $({\cal P}^g_\alpha)$, this solution must be unique.
 Since $\dim  {\cal S}^+_c = r+1 \le n+1$, the problem $({\cal P}^g_\alpha)$ is easier than $({\cal P}^d_\beta)$.
 In the case that  $({\cal P}^g_\alpha)$ does not have a solution to satisfy  \eqref{eq-exis},
we can oppositely  chose $\mbox{\boldmath$\alpha$} \in \mathbb{R}^n$  such that ${\bf H} \prec 0$. But in this case, $ P^g_\alpha(\mbox{\boldmath$\zeta$})$ is a d.c. function (difference of convex functions) and the corresponding problem should be \cite{gao-cace,gao-watsonetal}
\begin{equation}
({\cal P}^g_\alpha): \;\; \min \mbox{sta}  \left\{ P^g_\alpha(\mbox{\boldmath$\zeta$}) | \;\; \mbox{\boldmath$\zeta$} \in {\cal S}^+_c \right\},
\end{equation}
where $ \min \mbox{sta} \{ f(x)\} $ means to find the minimum stationary point of $f(x)$.
This is a nonsmooth d.c. programming problem.
A so-called  VTDIRECT parallel algorithm has been used successfully by Gao {\em et al}
 \cite{gao-watsonetal} for finding global optimum solutions to quadratic integer programming problems, but not in polynomial time.

For the  linear Knapsack problem, the improved canonical dual has a very simple form:
\begin{equation}
({\cal P}^g_{\ell }): \;\;\; \max_{ \tau \ge 0} \left\{  P^g_{\ell }(\tau) = -\frac{1}{2} \sum_{i=1}^n ( | c_i - \tau v_i |  -  \tau v_i)  - \tau V_c  \right\}. \label{eq-lknap}
\end{equation}
\begin{Corollary}[Improved Analytic  Solution  to $({\cal P}_\ell)$]
For any given  $V_c > 0$, $\mathbf{v}, \mathbf{c}  \in \mathbb{R}^n_+ $,    if $\tau_c >  0$ is a solution to $({\cal P}^g_{\ell })$
and $\theta_i= \tau_c   v_i - c_i  \neq 0 \;\forall i\in \{ 1,\dots, n\}$, then the linear knapsack problem $({\cal P}_\ell)$
has a  unique  global optimal solution
\begin{equation} \label{eq-rtau}
\mathbf{z}_c= \frac{1}{2} \left\{  \frac{c_i - \tau_c  v_i }{|c_i - \tau_c   v_i |} + 1 \right\}^n
\end{equation}
 and
$ P_{\ell} (\mathbf{z}_c) = P^g_{\ell } (\tau_c)$.
\end{Corollary}

Actually, from the proof of Corollary \ref{thm4} we know that if there exists a $\tau_\beta > 0$ such that
 $ \theta_i = \tau_\beta v_i - c_i\neq 0 \;\; \forall i \in \{ 1, \dots, n\}$, the canonical dual algebraic equation \eqref{eq-cdas} has a unique positive solution $\mbox{\boldmath$\sigma$}_\beta$
\cite{gao-book00,gao-na00,gao-jimo07}.
For a sufficiently big $\beta \gg 0$ we have $\mbox{\boldmath$\sigma$}_\beta = \mbox{\boldmath$\sigma$}_c$ and
 the solution $\tau_\beta$ by \eqref{eq-solust} is exactly the solution $\tau_c$ of the problem
 $({\cal P}^d_{\ell\alpha})$ (see \cite{gao-to18}). Therefore, $({\cal P}_\ell)$ has a unique solution defined by either \eqref{eq-solu} or \eqref{eq-rtau}.

\begin{remark} [Criteria for Non-Uniqueness and NP-Hardness]
Theoretical results presented so far show that if  the condition \eqref{eq-corc}  holds, the canonical dual problem $({\cal P}_q^d)$ has a unique solution
$\mbox{\boldmath$\zeta$}_c = (\mbox{\boldmath$\sigma$}_c, \tau_c) \in {\cal S}^+_a$, which can be obtained deterministically in polynomial time since  $({\cal P}_q^d)$ is equivalent to a convex minimization problem.
 In this case, both the quadratic and the  linear  knapsack problems are  not NP-hard and  their solutions $\mathbf{z}_c$ can be analytically given by Theorems \ref{thm2} and Corollary \ref{thm4}, respectively. Also, we must have $\mathbf{z}_c^T \mathbf{v} = V_c$ since $\tau_c > 0$.
 For the linear knapsack problem $({\cal P}_\ell)$, the  condition \eqref{eq-corc} is simply $\tau_c \mathbf{v} \neq \mathbf{c}$.

On the other hand, if the canonical dual problem $({\cal P}_q^d)$ has no solution in ${\cal S}^+_a$, the primal problem $({\cal P}_q)$ could be   NP-hard,
which is a conjecture first proposed by Gao in 2007 \cite{gao-jimo07}, i.e.
\begin{verse}
{\bf Conjecture of NP-Hardness}. A global optimization problem is NP-hard only if its canonical dual has no solution in ${\cal S}^+_a$.
\end{verse}
 It is also an open problem left in \cite{gao-cace,gao-ruan-jogo10}.
The reason for   NP-hard problems and possible   solutions  were discussed recently in \cite{gao-aip,gao-to18}.
For the linear knapsack problem $({\cal P}_\ell)$, as long as  $\tau_c v_i = c_i$ for any one $i \in \{1, \dots, n\}$,
  the canonical dual algebraic equation \eqref{eq-cdas} has at least  two repeated solutions $\mbox{\boldmath$\sigma$}_c$ located on the boundary of ${\cal S}^+_a$.
In this case, the primal  problem $({\cal P}_\ell)$ has multiple solutions, i.e. it is not well-posed \cite{gao-aip,gao-to18}.
A linear perturbation method for   solving this case was proposed recently in \cite{chen-gao-amma,gao-to18,wang-etal}.
\end{remark}

\section{Volume Reduction Method and Canonical  Duality Algorithm}

Theoretically speaking, for any given $V_c < V_0 = \sum_{i=1}^n v_i$,  the   canonical penalty-duality method can produce global optimal solution to the bilevel minimization
problem $({\cal P})$. However, if  $ \mu_c = V_c/V_0 \ll 1, $   any   iteration  method  could lead to unreasonable numerical solutions. In order to resolve this problem,
a volume reduction method  (VRM) method can be proposed:
\begin{verse}
(i) Introduce a volume reduction control parameter
$\mu \in (\mu_c,1)$ to produce a volume reduction sequence  $\{ V_{\gamma } = \mu V_{\gamma-1}\} $ ($\gamma =  1, \dots, \gamma_c$)
such that $V_{\gamma_c} = V_c$, ${\cal Z}_\gamma = \{ \mathbf{z} \in \{0,1\}^n| \;\;  \mathbf{v}^T {\bf z}} %{\mbox{\boldmath$\rho$} \le V_\gamma\}$; \\
(ii) For a  given $V_\gamma \in [V_c, V_0]$ and  initial values $ (\mathbf{x}^{\gamma-1}, {\bf z}} %{\mbox{\boldmath$\rho$}^{\gamma-1}) \in {\cal X}_a \times {\cal Z}_a$,
to find $(\mathbf{x}^\gamma, {\bf z}} %{\mbox{\boldmath$\rho$}^\gamma)$ such that
 \end{verse}
 \begin{eqnarray*}
({\cal P})^\gamma: \;\;    (\mathbf{x}^\gamma, {\bf z}} %{\mbox{\boldmath$\rho$}^\gamma) &   = &    \arg \min \left\{ P (\mathbf{x},{\bf z}} %{\mbox{\boldmath$\rho$}) |\; \mathbf{x}\in{\cal X}_a, \;\; \mathbf{z} \in {\cal Z}_\gamma \right\}, \\
 s.t. \;\; \mathbf{x}   & = & \arg \min \{\Pi(\mbox{\boldmath$\chi$}, \mathbf{z}) | \;\; \mbox{\boldmath$\chi$} \in {\cal X}_a \}.
\end{eqnarray*}

Generally  speaking,  for any given sequence $\{V_\gamma \}$ we should have
 \begin{equation}
 ({\cal P}) = \lim_{\gamma\rightarrow \gamma_c} ({\cal P})^\gamma.
 \end{equation}
Numerically, different volume sequence $\{V_\gamma\}$ may produce totally different numerical results as long as the decoupled  alternative iteration is used.
This is an  intrinsic difficulty for all   coupled bi-level  optimal design problems.
Based on the decoupled alternative iteration (DAI) and the volume reduction method (VRM),
 the canonical  duality  algorithm (CDT)  for solving the general bi-level optimization  problem $({\cal P})$ can be proposed.

\begin{algorithm}
\caption{Canonical  Dual Algorithm for Bi-Level Knapsack Problem (CDT)}
\begin{algorithmic} [1]
\State  \textbf{Input parameters:}  $\mu$,     error allowances  $\omega  > 0$,
  ${\bf z}} %{\mbox{\boldmath$\rho$}^0 = \{1\} \in \mathbb{R}^n$,   and     $V_0 = \sum_{i=1}^n v_i  $. Let $\gamma = 1 $.
\State  \textbf{Solve}  the lower-level optimization problem  \eqref{eq-lowk}:
  \begin{equation}
  \mathbf{x}^\gamma =  \arg \min \{ \Pi( \mathbf{x}, \mathbf{z}^{\gamma-1})   \;\; | \;\; \mathbf{x}\in {\cal X}_a \} .
  \end{equation}
\State \textbf{Compute} $\mathbf{Q}  = \mathbf{Q}(\mathbf{x}^{\gamma })$,  ${\bf c}  = {\bf c}(\mathbf{x}^\gamma)$, and $V_{\gamma } = \max\{ V_c,  \mu V_{\gamma-1}\} $.

  \State \textbf{Solve}  the canonical  dual  problem $({\cal P}^d_q)^\gamma$ to obtain
   \begin{equation}
\mbox{\boldmath$\zeta$}^\gamma =    (\mbox{\boldmath$\sigma$}^\gamma, \tau^\gamma) = \arg   \min \{  P^d_q (   \mbox{\boldmath$\zeta$})    | \; \;\; \mbox{\boldmath$\zeta$} \in {\cal S}^+_a \}
  \end{equation}

     \State   \textbf{Compute} the upper-level  solution   $\mathbf{z}^{\gamma}   $ by
        \[
        \mathbf{z}^{\gamma }   = \mathbf{G}(\mbox{\boldmath$\sigma$}^\gamma)^{-1} \mbox{\boldmath$\psi$}(\mbox{\boldmath$\zeta$}^\gamma)
        \]
      \State    {\bf  If }   $| P_q (  {\bf z}} %{\mbox{\boldmath$\rho$}^{\gamma })  -  P_q({\bf z}} %{\mbox{\boldmath$\rho$}^{\gamma-1})   | \le \omega $ and $V_\gamma  \le V_c$ ,
      then stop; Otherwise, continue.
                \State \textbf{Let}       $\gamma = \gamma+1 $. Go to Step 2.
\end{algorithmic} \label{alg-CDT}
\end{algorithm}
The canonical dual problem $({\cal P}^d_q)^\gamma$ in this algorithm could be either $({\cal P}^g_\alpha)^\gamma$
or $({\cal P}^d_\beta)^\gamma$  for quadratic knapsack problem $({\cal P}_q)$.
For linear knapsack problem $({\cal P}_\ell)^\gamma$, its canonical dual could be either
 $({\cal P}^g_{\ell })^\gamma$ or $({\cal P}^d_{\ell\beta})^\gamma$.

\section{Application  to  Optimal Topology Design}\label{sec-example}
Topology design is the arrangement of the various elements (links, nodes, volumes, etc.) of a complex network or a discretized continuous system
  in multidisciplinary fields of
   communication, electronics,    optics,  structural, bio- and nano-mechanics (cf. \cite{bc,cml,gao-to17,kjy}).
 It was discovered recently by Gao   \cite{gao-to17} that the  optimal topology design  for  general  elastic structures  should be
  formulated as a bilevel knapsack problem:
   \begin{eqnarray}
   ({\cal P}_{to}): &    & \min\{ P(\mathbf{x},\mathbf{z}) = - {\bf c}(\mathbf{x})^T \mathbf{z} | \;\; \mathbf{x} \in {\cal X}_a, \;\;  \mathbf{z} \in {\cal Z}_a   \} \;\;\;\;\;\; \;\;\;\\
& &  s.t.  \;\; \mathbf{x} \in \arg \min \left\{ \Pi( \mbox{\boldmath$\chi$},  \mathbf{z})   \;\;   | \;\;  \mbox{\boldmath$\chi$} \in {\cal X}_a  \right\} .\label{eq-lowto}
\end{eqnarray}
The leader design variable $\mathbf{z} = \{ 0 , 1\}^n $ in this problem is the so-called density distribution such that if  $z_e = 1$, then  the $e$-th element of the elastic structure is  solid, otherwise, this element is void.
The follower  variable $\mathbf{x} \in {\cal X}_a \subset \mathbb{R}^m$ represents the  displacement vector, whose domain ${\cal X}_a$ is a convex set, in which the   boundary condition is given (i.e. certain components of $\mathbf{x}$ are fixed). The vector ${\bf c}(\mathbf{x}) =\{ c_e (\mathbf{x}_e) \}   $ and its component  represents the strain energy   in the $e$-th element.
For linear elastic structural, the total potential energy $\Pi(\mathbf{x},\mathbf{z})$ is a quadratic function of $\mathbf{x}$:
\begin{equation}
\Pi(\mathbf{x},\mathbf{z}) =  \frac{1}{2} \mathbf{x}^T \mathbf{K}(\mathbf{z} ) \mathbf{x} - \mathbf{x}^T \mathbf{f} = \sum_{e=1}^n  \frac{z_e}{2} \mathbf{x}^T_e \mathbf{K}_e \mathbf{x}_e   - \mathbf{x}^T \mathbf{f},
\end{equation}
 where $\mathbf{K}_e$  and $\mathbf{x}_e$ are  respectively  the stiffness matrix  and the nodal displacement vector of the $e$-th element;
 $\mathbf{v} =\{ v_e \} \in \mathbb{R}^n_+ $ and $v_e \ge 0$ represents the volume of the $e$-th element;
$V_c$ is a desired volume.
Since $\mathbf{K}_e $ is positive definite, we have ${\bf c}(\mathbf{x}) = \frac{1}{2}  \{   \mathbf{x}^T_e \mathbf{K}_e \mathbf{x}_e  \}^n \in \mathbb{R}^n_+ $ and the global stiffness matrix  $\mathbf{K} (\mathbf{z}) = \{ z_e \mathbf{K}_e \} \in \mathbb{R}^{m\times m} $ is also positive definite for any given $\mathbf{z} \in {\cal Z}_a$.
The given vector $\mathbf{f} \in \mathbb{R}^m$ is  the  external force.

Clearly, for a given displacement vector $\mathbf{x} \in {\cal X}_a$, the upper-level problem is a typical linear knapsack problem, which can be solved analytically
by the canonical duality theory.
While for a fixed design variable $\mathbf{z} \in {\cal Z}_a$, the lower-level minimization is the  well-known minimum total potential energy principle for linear elastic structures.
Since the total potential energy $\Pi(\mathbf{x},\mathbf{z})$ is a quadratic function of $\mathbf{x}$, the lower-level solution can be given analytically by  $\mathbf{x} = \mathbf{K}(\mathbf{z})^{-1}  \mathbf{f}$.
Thus, this  challenging topology design problem can be solved by the combination of the DAI-VRM and the canonical  duality method.
Since the linear knapsack problem can be solved analytically by Corollary 1 or 2, the computational complexity for solving
the lower-level solution  $\mathbf{x} = \mathbf{K}(\mathbf{z})^{-1}  \mathbf{f}$ is about $O(m^3)$, the CDT is a polynomial time algorithm.

The proposed CDT algorithm  for topology design
 has been implemented in Matlab.
 The test examples are   the 2D and 3D   benchmark   cantilever beams    in optimal topology design of elastic structure
   (see  Fig.~\ref{fig-cant}).
\begin{figure}[tb]
  \centering
\subfigure[2D beam]{\includegraphics[width=.45\textwidth, height=50mm]{cant.jpg}    }
\subfigure[3D beam]{\includegraphics[width=.45\textwidth, height=50mm] {can3d.jpg}  }
     \caption{Design domains   for   long cantilever   beams with external load}  \label{fig-cant}
 \end{figure}
Performance of the  CDT   method  is first  tasted   for different mesh resolutions.
Results in Fig. \ref{fig1} show    that for any given mesh resolutions,
 the CDT   method produces precise integer solutions.
Clearly,   the finer the resolution, the smaller  the stored energy $C=\mathbf{z}^T \mathbf{c}(\mathbf{x})$  with better material distributions.
This means that the fine structure has more capacity for external load.
For the mesh resolution  $180\times 60$ elements, we have $n=180\times 60 = 10,800$ discrete variables and $m=2\times (180+1) \times (60+1) = 21,960$ continuous variables, but the total computing time by a   HP labtop computer (with Processor Intel
Core I7-4810, CPU @ 2.80 GHz and memory 2.80 GB) is only $t = 7.69$ seconds.

  \begin{figure} [h]
  \centering
  \subfigure [mesh=$40 \times 10$,  $C=416.577,$ Time $=1.2165$sec] {
\includegraphics[width=.45\textwidth, height=30mm]{fig1a.jpg}}
  \subfigure[mesh=$100 \times 30$,  $C=232.6$, Time$=2.96$sec] {
\includegraphics[width=.45\textwidth, height=30mm]{fig1b.jpg}}
\subfigure[mesh=$180 \times 60$,   $C=171.77,$ Time $=7.69$sec ] {\includegraphics[width=.45\textwidth, height=30mm]{fig1c.jpg}}
\caption{Optimal topology by different mesh resolutions with  $\mu_c = V_c/V_0 =  0.5$ and  $\mu=0.975$. Discrete variables $n=$mesh, continuous variables $m=2 (nex+1) (ney+1)$}
   \label{fig1}
\end{figure}

To compare with the commercial  codes  TOP88 (i.e. SIMP) and BESO \cite{Xie-Steven2}, we use the same mesh resolution of $180\times 60$.
The volume reduction rate $\mu=0.97$ for both BESO and CDT.
Computational results are reported in
Figure \ref{fig2}, which show clearly that the CDT method produces
geometrically simple and mechanically sound structure.
Since the SIMP method is  a continuous relaxation approach, which cannot produce integer solutions,
the structure obtained by the SIMP code is ugly with a large area of grayscaling checkerboard  patterns.
By the fact that the BESO code is an evolutionary method for solving the linear knapsack problem by simply using comparison algorithm, although it can produce integer solution similar to that by CDT, it is not a polynomial time algorithm \cite{gao-to18},
therefore, it takes much longer time  (more than 30 times) than the CDT  method.
 \begin{figure} [h]
  \centering
   \subfigure[computational result by SIMP] {\includegraphics[width=.45\textwidth, height=40mm] {TOP180-60.jpg} }
 \subfigure[computational result by BESO] {\includegraphics[width=.45\textwidth, height=40mm] {BESO180-60.jpg} }
 \subfigure[Computational result by CDT]   {   \includegraphics[width=0.45\textwidth, height=40mm]{Fig7c1.jpg} }
  \caption{2D beam structures by SIMP (a), BESO (b) and CDT (c) with  $n=180\times 60$, $m=2(180+1) \times (60+1)$.
The CDT is   twice faster than the SIMP, and about 47 times faster than BESO}
  \label{fig2}
\end{figure}

In order to see the structures' evolution  during the iterations,
we tested the 2D beam with $V_c=0.5, \;\;  nex = 150, \;\; ney=50$.
Fig. \ref{fig2d3} shows the outputs of the structures by all methods  at the fifth, 15th and the final iterations.
Since the initial volume is $V_0 = 1$, the volume reduction rate is $\mu=0.97$ for both CDT and BESO,
therefore, the volumes at the 5-th and the 15-th iteration should be $V_5=0.859$ and $V_{15} = 0.633$, respectively.
 Results  in Fig. \ref{fig2d3} (a) and (b) show clearly that the structures produced by CDT algorithm are
final knapsack solutions if the desired knapsacks $V_c$  are $V_5 $ and $V_{15}$, respectively.
However, the structures produced by BESO are   broken with disconnected branches Fig. \ref{fig2d3}(b),
the structures by SIMP are uncertain with large areas of grayscales.
\begin{figure*}
\begin{tabularx}{0.5\textwidth}{ccc}
\hspace{-1.3em
{
\begin{tabular}{lc}
\small  \subfigure[ CDT: $V_5=0.8587$, Time=4.97] {\includegraphics[scale=.35]{cdt1.eps}}\vspace{.4cm}\\
\small  \subfigure[ BESO: $V_5=0.8585$, Time=23.66] { \includegraphics[scale=.35]{beso1.eps}}\vspace{.3cm} \\
\small   \subfigure[SIMP:  Itration=5,  Time=3.07] {\includegraphics[scale=.35]{simp1.eps}}
\end{tabular}}
&
\hspace{-2.5em
{
\begin{tabular}{c}
\small \subfigure[ $V_{15}=0.63325$, Time=14.1421] {\includegraphics[scale=.35]{cdt2.eps}}\\
\small \subfigure[ $V_{15}=0.63343$, Time=73.812] {\includegraphics[scale=.35]{beso2.eps}}\\
\small \subfigure[ Iteration=15, Time=11.1871] {\includegraphics[scale=.35]{simp2.eps}}
\end{tabular}}&
\hspace{-2.5em
{
\begin{tabular}{c}
\small \subfigure[ $V_c=0.5$, Time=25.4657] {\includegraphics[scale=.35]{cdt3.eps}}\\
\small \subfigure[ $V_c=0.5$, Time=116.332] {\includegraphics[scale=.35]{beso3.eps}}\\
\small \subfigure[ Final Iteration, Time=104.0837] {\includegraphics[scale=.35]{simp3.eps}}
\end{tabular}}
\end{tabularx}
\caption{Topology evolutions for 2D  beam by CDT (top), BESO (middle), and SIMP (bottom)
with  $V_c=0.5$ and mesh resolution $150 \times 50$. The CDT is about four times faster than both BESO and SIMP.
The SIMP can produce only grayscale  topology since it is a continuous approach, i.e. $\mathbf{z} \in (0, 1)^n$}\label{fig2d3}
\end{figure*}

 In order to look the   energy density  distribution $\mathbf{c} = \{ c_e(\mathbf{x}) \}$ in the optimal structures, we use the mesh resolution  $nex\times ney = 80 \times 30$.
 Fig. \ref{fig23} shows clearly that the CDT can produce mechanically sound structure with homogeneous distribution of strain energy density. Since the material distribution produced by SIMP is not mechanically reasonable, the overall strain energy stored in each element is about five times  higher than that by CDT (see the color bars in Fig. \ref{fig23}).  If we consider the dark-red level (= 1)  is the elastic limit, then the structure produced by SIMP is far beyond this limit.
This shows that the optimal structure by CDT has much potential to support even more  external load.
Although the BESO can produce better result than the SIMP, it can be used only for small-sized problems since it is not a polynomial time algorithm.
  \begin{figure} [h]
  \centering
  \subfigure[SIMP: $C=178.688$] {
\includegraphics[width=.5\textwidth, height=40mm]{Fig3a.jpg}}
  \subfigure[BESO: $C=176.494$] {
\includegraphics[width=.5\textwidth, height=40mm]{Fig3b.jpg}}
\subfigure[CDT: $C=170.869$] {\includegraphics[width=.5\textwidth, height=40mm]{Fig3c.jpg}}
\caption{Structures and energy density distributions  by SIMP (a), BESO (b) and CDT (c) with $V_c=0.4$.}
  \label{fig23}
\end{figure}

For the 3D beam,  we only compare the CDT with the SIMP since the BESO is computationally expensive.
First, we use the mesh resolution  $40\times 15 \times 4$ and volume fraction  $V_c = 0.2$.
Fig. \ref{fig3D1} shows that the optimal topology produced by CDT is not only geometrically elegant, but also mechanically sound with homogeneous distribution of energy density and very few red elements. From the color bar we can see that the
 red-level is only at scale of 26.34. However, the topology produced by SIMP is  overstaffed with many red elements at scale of 69.4.   This shows that the structure produced by the CDT can support external load three times  than the structure by SIMP.

\begin{figure}
\subfigure[SIMP solution: C=4465.3]{\includegraphics[scale=0.45]{top3du3.jpg}}\quad
\subfigure[SIMP solution: Front view]{\includegraphics[scale=0.45]{top3du3f.jpg}}
 \subfigure[CDT solution: C=2792]{\includegraphics[scale=0.45]{cdt3dukp.jpg}}\quad
\subfigure[Front view of CDT solution]{\includegraphics[scale=0.4]{cdt3dukpf.jpg}}
\caption{Structures  for 3-D   beam with mesh=$40\times 15 \times 4$   and $V_c =0.2$. The SIMP beam is overstaffed with  many red elements at scale=69.4 (highly stressed).
The CDT beam is geometrically elegant with very few red elements at scale=26.34, so it can support more external load. }\label{fig3D1}
\end{figure}

Now we increase the mesh resolution to
   $60\times 20 \times 10$ and reduce the volume fraction to $V_c = 0.1$.
  Fig. \ref{fig4} shows results   obtained by  the commercial code TOP3D (SIMP)  and CDT3D     (with $\mu=0.93$).
Since  the  0-1 density distribution produced  by the  CDT is much more mechanically reasonable than the one by SIMP,
most of elements in Fig. \ref{fig4}(a)  are in blue color with only a few red elements at the scale= 88.
The beam by   SIMP has more red  elements  at the scale =117.
The structure shown in Fig. \ref{fig4}(c) is the out put of TOP3D at the  $It=50$-th iteration, which is almost the same as the final result
($It=199$) shown in Fig. \ref{fig4}(e). This result shows that the SIMP has a very slow rate of  convergence.
Also the CDT is more than 10 times faster than the TOP3D.
The uncolored Fig. \ref{fig5} shows  that most   elements  are grayscale (i.e. $0< z_i < 1$)  since the SIMP is a continuous approximation.

\begin{figure}[tb]
\centering
\subfigure[CDT solution: C=13192.7, Time=1.56min]{
\includegraphics[width=.45\textwidth] {fig4c.jpg} }\quad
\subfigure[Front view of CDT solution ]{
\includegraphics[width=0.4\textwidth, height=35mm] {fig4d.jpg} }
\qquad
\subfigure[SIMP solution at the $50$-th iteration: C=21960, Time=12.2min]{
\includegraphics[width=.45\textwidth] {fig4a.jpg} }\quad
\subfigure[Front view of SIMP solution ]{
\includegraphics[width=0.4\textwidth, height=45mm] {fig4b.jpg} }
\caption{3D beam structures  by CDT and SIMP with $V_c=0.1$
and mesh resolution  $ 60 \times 20 \times 10$. Since the strain energy stored in the CDT beam ($C=13192.7$)  is only
about a  half of that stored in the SIMP beam ($C=21960$), the CDT beam
 has much more potential to support  more external load.}
\label{fig4}
\end{figure}

\begin{figure}[tb]
\centering
\subfigure[SIMP solution ($It=199$): C=21426, Time=19.8min]{
\includegraphics[width=.40\textwidth] {fig4a1.jpg} }\quad
\subfigure[Front view of SIMP solution ]{
\includegraphics[width=0.4\textwidth, height=35mm] {fig4b1.jpg} }
\caption{3D beam structures  by  SIMP with $V_c=0.1$
and mesh resolution  $ 60 \times 20 \times 10$.  }
\label{fig5}
\end{figure}

Detailed study on the canonical duality theory for solving topology design  problems and a 66-line  Matlab code  of the CDT algorithm   will be  given   in
\cite{gao-r-l-to}.

\section{Concluding remarks  and open problems}\label{sec-con}
Based on a  combination of an alternative iteration, a volume reduction and    the novel canonical duality theory,
a new powerful   method  is proposed   for solving  bilevel optimization  governed by knapsack problems.
The main theoretical contributions are presented in Sections  3 and 4, i.e.  by using  the canonical  duality method,  the well-known knapsack problems can be  solved
 analytically in terms of their canonical dual solutions.
The existence and uniqueness of this general analytical solution are proved for both quadratic and linear knapsack problems.
It shows that for any given ${\bf c},\;\mathbf{v}\in \mathbb{R}^n$, if there exists a canonical dual solution  $\tau_c >0$ such that $\tau_c {\bf c} \neq \mathbf{v}$, the  linear knapsack problem
can be solved  uniquely in polynomial time.
Therefore, a polynomial  algorithm CDT  is proposed for solving the bilevel knapsack problems.
Its novelty is demonstrated by solving a benchmark  problem  in topology design.
Numerical results indicated   that the volume reduction method is essential for solving bilevel optimization  governed by knapsack problems.
Both theoretical and numerical results verified  that the   CDT  is a powerful method for studying
general bilevel  optimization problems.

Due to the intrinsic  coupling effect in  the  bilevel  optimization  problem $({\cal P})$,
  although  for each  given lower level variable $\mathbf{x}_k\in {\cal X}_a$,  the  upper level  knapsack problem  $\mathbf{z}_k = \arg \min \{ P(\mathbf{x}_k, \mathbf{z}) | \; \; \mathbf{z}   \in {\cal Z}_a  \} $
  can be solved analytically by the canonical duality theory, the solution sequence $\{\mathbf{x}_k,\mathbf{z}_k\}$  may not converge to  the global optimal solution of  $({\cal P})$
  since the decoupled alternative iteration method is adopted
 in Algorithm 1.
This is the main reason that  the numerical solutions produced by the CDT  depend  sensitively on  the  volume reduction  rate $\mu$.
 This is an open problem and deserves   theoretical study in the future.

\section*{Acknowledgement}
The author would like to sincerely acknowledge the important comments and suggestions from  the editor and anonymous
reviewers.
This paper is based on several  invited/keynote/plenary
lectures presented at international conferences including
the Symposia on Intelligent Technologies for Advancing and Safeguarding Australia, 15 August 2017, Deakin University,
and the  Second In. Conf.  on Modern Mathematical Methods and High Performance Computing in Science and Technology, 4-6 January, 2018, New Delhi, India.
 Organizers' hospitality and financial supports are sincerely acknowledged.
 This research is supported by US Air Force Office for Scientific Research (AFOSR)  under the grants   FA2386-16-1-4082 and FA9550-17-1-0151.