Quadratic knapsack problem

Application

As one might expect, QKP 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. Applications of QKP in the field of computer science is more common after the early days: compiler design problem, clique problem, very large scale integration (VLSI) design. Additionally, pricing problems appear to be an application of QKP as described by Johnson et al.

Computational complexity

In general, the decision version of the knapsack problem (Can a value of at least V be achieved under a restriction of a certain capacity W?) is NP-complete. Thus, a given solution can be verified in polynomial time while no algorithm can identify a solution efficiently.

The optimization knapsack problem is NP-hard and there is no known algorithm that can solve the problem in polynomial time.

As a particular variation of the knapsack problem, the 0-1 quadratic knapsack problem is also NP-hard.

While no available efficient algorithm exists in the literature, there is a pseudo-polynomial time based on dynamic programming and other heuristic algorithms that can always generate “good” solutions.

Solving

While the knapsack problem is one of the most commonly solved operation research (OR) problems, there are limited efficient algorithms that can solve 0-1 quadratic knapsack problems. Available algorithms include but are not limited to brute force, linearization, and convex reformulation. Just like other NP-hard problems, it is usually enough to find a workable solution even if it is not necessarily optimal. Heuristic algorithms based on greedy algorithm, dynamic programming can give a relatively “good” solution to the 0-1 QKP efficiently.

Brute force

The brute-force algorithm to solve this problem is to identify all possible subsets of the items without exceeding the capacity and select the one with the optimal value. The pseudo-code is provided as follows:

<span id="mwXg"></span><span class="c1" id="mwXw">// Input:</span> <span class="c1" id="mwYA">// Profits (stored in array p)</span> <span class="c1" id="mwYQ">// Quadratic profits (stored in matrix P)</span> <span class="c1" id="mwYg">// Weights (stored in array w)</span> <span class="c1" id="mwYw">// Number of items (n)</span> <span class="c1" id="mwZA">// Knapsack capacity (W)</span> <span class="kt" id="mwZQ">int</span><span class="w" id="mwZg"> </span><span class="n" id="mwZw">max</span><span class="w" id="mwaA"> </span><span class="o" id="mwaQ">=</span><span class="w" id="mwag"> </span><span class="mi" id="mwaw">0</span> <span class="k" id="mwbA">for</span><span class="w" id="mwbQ"> </span><span class="n" id="mwbg">all</span><span class="w" id="mwbw"> </span><span class="n" id="mwcA">subset</span><span class="w" id="mwcQ"> </span><span class="n" id="mwcg">S</span><span class="w" id="mwcw"> </span><span class="k" id="mwdA">do</span> <span class="w" id="mwdQ"> </span><span class="kt" id="mwdg">int</span><span class="w" id="mwdw"> </span><span class="n" id="mweA">value</span><span class="p" id="mweQ">,</span><span class="w" id="mweg"> </span><span class="n" id="mwew">weight</span><span class="w" id="mwfA"> </span><span class="o" id="mwfQ">=</span><span class="w" id="mwfg"> </span><span class="mi" id="mwfw">0</span> <span class="w" id="mwgA"> </span><span class="k" id="mwgQ">for</span><span class="w" id="mwgg"> </span><span class="n" id="mwgw">i</span><span class="w" id="mwhA"> </span><span class="n" id="mwhQ">from</span><span class="w" id="mwhg"> </span><span class="mi" id="mwhw">0</span><span class="w" id="mwiA"> </span><span class="n" id="mwiQ">to</span><span class="w" id="mwig"> </span><span class="n" id="mwiw">S</span><span class="p" id="mwjA">.</span><span class="n" id="mwjQ">size</span><span class="mi" id="mwjg">-1</span><span class="w" id="mwjw"> </span><span class="k" id="mwkA">do</span><span class="o" id="mwkQ">:</span> <span class="w" id="mwkg"> </span><span class="n" id="mwkw">value</span><span class="w" id="mwlA"> </span><span class="o" id="mwlQ">=</span><span class="w" id="mwlg"> </span><span class="n" id="mwlw">value</span><span class="w" id="mwmA"> </span><span class="o" id="mwmQ">+</span><span class="w" id="mwmg"> </span><span class="n" id="mwmw">p</span><span class="p" id="mwnA">[</span><span class="n" id="mwnQ">i</span><span class="p" id="mwng">]</span> <span class="w" id="mwnw"> </span><span class="n" id="mwoA">weight</span><span class="w" id="mwoQ"> </span><span class="o" id="mwog">=</span><span class="w" id="mwow"> </span><span class="n" id="mwpA">weight</span><span class="w" id="mwpQ"> </span><span class="o" id="mwpg">+</span><span class="w" id="mwpw"> </span><span class="n" id="mwqA">w</span><span class="p" id="mwqQ">[</span><span class="n" id="mwqg">i</span><span class="p" id="mwqw">]</span> <span class="w" id="mwrA"> </span><span class="k" id="mwrQ">for</span><span class="w" id="mwrg"> </span><span class="n" id="mwrw">j</span><span class="w" id="mwsA"> </span><span class="n" id="mwsQ">from</span><span class="w" id="mwsg"> </span><span class="n" id="mwsw">i</span><span class="o" id="mwtA">+</span><span class="mi" id="mwtQ">1</span><span class="w" id="mwtg"> </span><span class="n" id="mwtw">to</span><span class="w" id="mwuA"> </span><span class="n" id="mwuQ">S</span><span class="p" id="mwug">.</span><span class="n" id="mwuw">size</span><span class="mi" id="mwvA">-1</span><span class="w" id="mwvQ"> </span><span class="k" id="mwvg">do</span><span class="o" id="mwvw">:</span><span class="w" id="mwwA"> </span> <span class="w" id="mwwQ"> </span><span class="n" id="mwwg">value</span><span class="w" id="mwww"> </span><span class="o" id="mwxA">=</span><span class="w" id="mwxQ"> </span><span class="n" id="mwxg">value</span><span class="w" id="mwxw"> </span><span class="o" id="mwyA">+</span><span class="w" id="mwyQ"> </span><span class="n" id="mwyg">P</span><span class="p" id="mwyw">[</span><span class="n" id="mwzA">i</span><span class="p" id="mwzQ">][</span><span class="n" id="mwzg">j</span><span class="p" id="mwzw">]</span> <span class="w" id="mw0A"> </span><span class="k" id="mw0Q">if</span><span class="w" id="mw0g"> </span><span class="n" id="mw0w">weight</span><span class="w" id="mw1A"> </span><span class="o" id="mw1Q"><=</span><span class="w" id="mw1g"> </span><span class="n" id="mw1w">W</span><span class="w" id="mw2A"> </span><span class="n" id="mw2Q">then</span><span class="o" id="mw2g">:</span> <span class="w" id="mw2w"> </span><span class="k" id="mw3A">if</span><span class="w" id="mw3Q"> </span><span class="n" id="mw3g">value</span><span class="w" id="mw3w"> </span><span class="o" id="mw4A">></span><span class="w" id="mw4Q"> </span><span class="n" id="mw4g">max</span><span class="w" id="mw4w"> </span><span class="n" id="mw5A">then</span><span class="o" id="mw5Q">:</span><span class="w" id="mw5g"> </span> <span class="w" id="mw5w"> </span><span class="n" id="mw6A">max</span><span class="w" id="mw6Q"> </span><span class="o" id="mw6g">=</span><span class="w" id="mw6w"> </span><span class="n" id="mw7A">value</span>

Given n items, there will be at most ${\displaystyle 2^n}$ subsets and for each legal candidate set, the running time of computing the values earned is ${\displaystyle O(n^2)}$. Thus, the efficiency class of brute-force algorithm is ${\displaystyle (2^n{n}^2)=\lambda (2^n)}$, being exponential.

Linearization

Problems of such form are difficult to solve directly using standard solvers and thus people try to reformulate it as a linear program using auxiliary variables and constraints so that the problem can be readily solved using commercial packages. Two well-known linearization approaches for the 0-1 QKP are the standard linearization and Glover's linearization.

Standard linearization

The first one is the standard linearization strategy, as shown below:

LP1: maximize

${\displaystyle \sum _{i=1}^n{p}_i{x}_i+\sum _{i=1}^n\left(\sum _{j=1,i<j}^n(P_{ij}+P_{ji})z_{ij}\right).}$

subject to

${\displaystyle z_{ij}\le x_i}$ for all ${\displaystyle (i,j),i<j}$

${\displaystyle z_{ij}\le x_j}$ for all ${\displaystyle (i,j),i<j}$

${\displaystyle x_i+x_j-1\le z_{ij}}$ for all ${\displaystyle (i,j),i<j}$

${\displaystyle z_{ij}\ge 0}$ for all ${\displaystyle (i,j),i<j}$

${\displaystyle x\in X,x}$ binary

In the formulation LP1, we have replaced the x_ix_j term with a continuous variable z_ij. This reformulates the QKP into a knapsack problem, which we can then solve optimally using standard solvers.

Glover's linearization

The second reformulation, which is more concise, is called Glover's linearization. The Glover formulation is shown below, where L_i and U_i are lower and upper bounds on ${\displaystyle \sum _{j=1,i\ne j}^n{P}_{ij}{x}_j}$, respectively:

LP2: maximize

${\displaystyle \sum _{i=1}^n{p}_i{x}_i+\sum _{i=1}^n{z}_i}$

subject to

${\displaystyle L_i{x}_i\le z_i\le U_i{x}_i}$ for ${\displaystyle i=1,\dots ,n}$

${\displaystyle \sum _{j=1,i\ne j}^n{P}_{ij}{x}_j-U_i(1-x_i)\le z_i\le \sum _{j=1,i\ne j}^n{P}_{ij}{x}_j-L_i(1-x_i)}$ for ${\displaystyle i=1,\dots ,n}$

${\displaystyle x\in X,x}$ binary

In the formulation LP2, we have replaced the expression ${\displaystyle \sum _{j=1,i\ne j}^n{P}_{ij}{x}_i{x}_j}$ with a continuous variable z_i. Similarly, we can use standard solvers to solve the linearized problem. Note that Glover's linearization only includes ${\displaystyle n}$ auxiliary variables with ${\displaystyle 2n}$ constraints while standard linearization requires ${\displaystyle (\genfrac{}{}{0ex}{}{n}{2})}$ auxiliary variables and ${\displaystyle 3(\genfrac{}{}{0ex}{}{n}{2})}$ constraints to achieve linearity.

Convex quadratic reformulation

Note that nonlinear programs are hard to solve due to the possibility of being stuck at a local maximum. However, when the program is convex, any local maximum is the global maximum. A convex program is to maximize a concave function or minimize a convex function on a convex set. A set S is convex if ${\displaystyle \mathrm{\forall }u,v\in S}$, ${\displaystyle \lambda u+(1-\lambda )v\in S}$ where ${\displaystyle \lambda \in [0,1]}$. That is to say, any point between two points in the set must also be an element of the set. A function f is concave if ${\displaystyle f[\lambda u+(1-\lambda )v]\le \lambda f(u)+(1-\lambda )f(v)}$. A function f is convex if ${\displaystyle f[\lambda u+(1-\lambda )v]\ge \lambda f(u)+(1-\lambda )f(v)}$. Informally, a function is concave if the line segment connecting two points on the graph lies above or on the graph, while a function is convex if below or on the graph. Thus, by rewriting the objective function into an equivalent convex function, we can reformulate the program to be convex, which can be solved using optimization packages.

The objective function can be written as ${\displaystyle c^{T}x+x^{T}Cx}$ using linear algebra notation. We need to make P a positive semi-definite matrix in order to reformulate a convex function. In this case, we modify the objective function to be ${\displaystyle p^{T}x+x^{T}Px+\sum _{i=1}^n\left(\sum _{j=1,j\ne i}^n|P_{ij}|\right)(x_i^2-x_i)}$ by applying results from linear algebra, where P is a diagonally dominant matrix and thus a positive semi-definite. This reformulation can be solved using a standard commercial mixed-integer quadratic package.

Greedy heuristic algorithm

George Dantzig proposed a greedy approximation algorithm to unbounded knapsack problem which can also be used to solve the 0-1 QKP. The algorithm consists of two phrases: identify an initial solution and improve it.

First compute for each item, the total objective contribution realizable by selecting it, ${\displaystyle p_i+\sum _{i\ne j}^n{P}_{ij}}$, and sort the items in decreasing order of the potential value per unit of weight, ${\displaystyle (p_i+\sum _{i\ne j}^n{P}_{ij})/w_i}$. Then select the items with the maximal value-weight ratio into the knapsack until there is no space for more, which forms the initial solution. Starting with the initial solution, the improvement is conducted by pairwise exchange. For each item in the solution set, identify the items not in the set where swapping results in an improving objective. Select the pair with maximal improvement and swap. There are also possibilities that removing one from the set or adding one to the set will produce the greatest contribution. Repeat until there is no improving swapping. The complexity class of this algorithm is ${\displaystyle O(2^n)}$ since for the worst case every possible combination of items will be identified.

Quadknap

Quadknap is an exact branch-and-bound algorithm proposed by Caprara et al., where upper bounds are computed by considering a Lagrangian relaxation which approximate a difficult problem by a simpler problem and penalizes violations of constraints using Lagrange multiplier to impost a cost on violations. Quadknap releases the integer requirement when computing the upper bounds. Suboptimal Lagrangian multipliers are derived from sub-gradient optimization and provide a convenient reformulation of the problem. This algorithm is quite efficient since Lagrangian multipliers are stable, and suitable data structures are adopted to compute a tight upper bound in linear expected time in the number of variables. This algorithm was reported to generate exact solutions of instances with up to 400 binary variables, i.e., significantly larger than those solvable by other approaches. The code was written in C and is available online.

Dynamic programming heuristic

While dynamic programming can generate optimal solutions to knapsack problems, dynamic programming approaches for QKP can only yield a relatively good quality solution, which can serve as a lower bound to the optimal objectives. While it runs in pseudo-polynomial time, it has a large memory requirement.

Dynamic programming algorithm

For simplicity, assume all weights are non-negative. The objective is to maximize total value subject to the constraint: that the total weight is less than or equal to W. Then for each ${\displaystyle w\le W}$, define ${\displaystyle f(m,w)}$ to be the value of the most profitable packing of the first m items found with a total weight of w. That is, let

${\displaystyle f(m,w)=max\left\{\sum _{i=1}^m{p}_i{x}_i+\sum _{i=1}^m\sum _{j=1,i\ne j}^m{P}_{ij}{x}_i{x}_j:\sum _{i=1}^m{w}_i=w,1\le i\le m\right\}.}$

Then, ${\displaystyle f(m,w)}$ is the solution to the problem. Note that by dynamic programming, the solution to a problem arises from the solution to its smaller sub-problems. In this particular case, start with the first item and try to find a better packing by considering adding items with an expected weight of 𝑤. If the weight of the item to be added exceeds 𝑤, then ${\displaystyle f(m,w)}$ is the same with ${\displaystyle f(m-1,w)}$. Given that the item has a smaller weight compared with the desired weight, ${\displaystyle f(m,w)}$ is either the same as ${\displaystyle f(m-1,w)}$ if adding makes no contribution, or the same as the solution for a knapsack with smaller capacity, specifically one with the capacity reduced by the weight of that chosen item, plus the value of one correct item, i.e. ${\displaystyle f(m-1,w-w_m)+p_m+\sum _{i=1}^{m-1}{P}_{im}{x}_i}$. To conclude, we have that

${\displaystyle f(m,w)=\{\begin{array}{ll}maxf(m-1,w),f(m-1,w-w_m)+p_m+\sum _{i=1}^{m-1}{P}_{im}{x}_i& \text{if }{w}_m\le w\\ f(m-1,w)& \text{otherwise}\end{array}}$

Note on efficiency class: Clearly the running time of this algorithm is ${\displaystyle O(W{n}^2)}$, based on the nested loop and the computation of the profit of new packing. This does not contradict the fact the QKP is NP-hard since W is not polynomial in the length of the input (and furthermore it provides no guarantee to reach the optimal solution).

Revised dynamic programming algorithm

Note that the previous algorithm requires ${\displaystyle O(W{n}^2)}$ space for storing the current packing of items for all m,w, which may not be able to handle large-size problems. In fact, this can be easily improved by dropping the index m from ${\displaystyle f(m,w)}$ since all the computations depend only on the results from the preceding stage.

Redefine ${\displaystyle f(w)}$ to be the current value of the most profitable packing found by the heuristic. That is,

${\displaystyle f(w)=max\left\{\sum _{i=1}^m{p}_i{x}_i+\sum _{i=1}^m\sum _{j=1,i\ne j}^m{P}_{ij}{x}_i{x}_j:\sum _{i=1}^m{w}_i=w,m\le n\right\}.}$

Accordingly, by dynamic programming we have that

${\displaystyle f(m)=\{\begin{array}{ll}maxf(w),f(w-w_m)+p_m+\sum _{i=1}^{m-1}{P}_{im}{x}_i& \text{if }{w}_m\le w,\\ f(w)& \text{otherwise.}\end{array}}$

Note this revised algorithm still runs in ${\displaystyle O(W{n}^2)}$ while only taking up ${\displaystyle O(Wn)}$ memory compared to the previous ${\displaystyle O(W{n}^2)}$.

Related research topics

Researchers have studied 0-1 quadratic knapsack problems for decades. One focus is to find effective algorithms or effective heuristics, especially those with an outstanding performance solving real world problems. The relationship between the decision version and the optimization version of the 0-1 QKP should not be ignored when working with either one. On one hand, if the decision problem can be solved in polynomial time, then one can find the optimal solution by applying this algorithm iteratively. On the other hand, if there exists an algorithm that can solve the optimization problem efficiently, then it can be utilized in solving the decision problem by comparing the input with the optimal value.

Another theme in literature is to identify what are the "hard" problems. Researchers who study the 0-1 QKP often perform computational studies to show the superiority of their strategies. Such studies can also be conducted to assess the performance of different solution methods. For the 0-1 QKP, those computational studies often rely on randomly generated data, introduced by Gallo et al. Essentially every computational study of the 0-1 QKP utilizes data that is randomly generated as follows. The weights are integers taken from a uniform distribution over the interval [1, 50], and the capacity constraints is an integer taken from a uniform distribution between 50 and the sum of item weights. The objective coefficients, i.e. the values are randomly chosen from [1,100]. It has been observed that generating instances of this form yields problems with highly variable and unpredictable difficulty. Therefore, the computational studies presented in the literature may be unsound. Thus some researches aim to develop a methodology to generate instances of the 0-1 QKP with a predictable and consistent level of difficulty.