On the Rectangular Knapsack Problem - Approximation of a Speciﬁc Quadratic Knapsack Problem

In this article, we introduce the rectangular knapsack problem as a special case of the quadratic knapsack problem consisting in the maximization of the product of two separate knapsack proﬁts subject to a cardinality constraint. We propose a polynomial time algorithm for this problem that provides a constant approximation ratio of 4.5. Our experimental results on a large number of artiﬁcially generated problem instances show that the average ratio is far from theoretical guarantee. In addition, we suggest reﬁned versions of this approximation algorithm with the same time complexity and approximation ratio that lead to even better experimental results.


Introduction
In contrast to classical (linear) knapsack problems, the profit of a collection of items in a quadratic knapsack problem is not only determined by the sum of individual profits, but also by profits generated by pairwise combinations of items. This can be used to model the fact that two items may complement each other such that their profit is increased if both of them are selected.
The quadratic objective still allows to model that the profit of two items is independent of each other by setting the combined profit to 0. In this case, including both items does not increase the profit over the sum of the individual profits. Furthermore, a negative combined profit value can model the fact that both items together are less profitable than the sum of individual profits. This might be the case, if both items are substitutes for each other and including both items is as profitable as including one.
The formulation of quadratic knapsack problems (QKP) is very general and, therefore, its range of applications is quite wide. For example, Johnson et al. [1993] present a problem in the context of compiler construction that can be formulated as a quadratic knapsack problem. Moreover, QKP have been discussed in the context of the location of airports, freight handling terminals, railway stations, and satellite stations [Rhys, 1970, Witzgall, 1975.
We present a variant of QKP which we call the rectangular knapsack problem (RKP). The profit matrix is built by the product of two vectors and the constraint is a cardinality constraint.
The main motivation for problem RKP arises when solving the cardinality constrained bi-objective knapsack problem (2oKP) where a, b ∈ N n , with a, b = 0 n = (0, 0, . . . , 0) ∈ N n , and k ∈ N, k < n. Instead of computing the set of efficient solutions for this bi-objective optimization problem, we want to find one (or several) representative nondominated point(s). Originally proposed by Zitzler and Thiele [1998] in the context of evolutionary algorithms, the hypervolume indicator is often used as a versatile quality measure of representation of the efficient set in multiobjetive optimization [c.f. Kuhn et al., 2016]. The problem of finding one solution of 2oKP that maximizes the hypervolume, considering (0, 0) as reference point, is equivalent to RKP. The structure of RKP allows to formulate a polynomial time 4.5-approximation algorithm. An algorithm is called a polynomial time ρ-approximation algorithm, if it computes a feasible solution in run time being polynomial in the encoding length of the input such that The remainder of this article is organized as follows: In Section 2, we give an introduction to quadratic knapsack problems. We introduce the rectangular knapsack problem in Section 3 and present upper and lower bounds. These bounds motivate an approximation algorithm that is formulated in Section 4, for which a constant approximation ratio ρ is proven. Furthermore, we also introduce improved implementations of this approximation method. In Section 5 we present a computational study of these algorithms and compare the realized approximation ratios to the theoretical bound of 4.5. Section 6 concludes this article.
2 Quadratic knapsack problems Gallo et al. [1980] first introduced the binary quadratic knapsack problem (QKP). It is a variant of the classical knapsack problem and can be concisely stated as follows: Given n items, the profit for including item i is given by the coefficient p ii . Additionally, a profit p ij + p ji is generated if both items i and j are selected. The values p ij (which are often assumed to be non-negative integers) can be compactly written in a profit matrix The profits p ij and p ji are either both realized, i. e., if items i and j are selected, or both not realized, i. e., if item i or item j is not selected. Hence, p ij and p ji can be assumed to be equally valued, which results in a symmetric matrix P .
As for the classical knapsack problem, each item i has a positive integral weight w i and the goal is to select a subset of items that maximizes the overall profit while the sum of weights does not exceed the given capacity W . As usual, the binary decision variable x i indicates if item i is selected, x i = 1, or not, x i = 0. Thus, QKP can be defined as follows: An illustrative interpretation of QKP can be given based on graphs. We define a complete undirected graph G = (V, E) with vertex set V = {1, . . . , n}, i. e., each vertex corresponds to one item of QKP. Each vertex has assigned a profit value p ii and a weight value w i and each edge (i, j) has assigned a profit value p ij + p ji . The task is to select a clique S ⊂ V with maximal profit which does not exceed the capacity W . The overall profit of S consists of the profit of the vertices in S and edges (i, j) connecting two vertices i, j ∈ S. It is well known that the quadratic knapsack problem is N P-complete in the strong 4 B. Schulze et al.
sense, which can be shown by a polynomial reduction from the Clique-problem [Garey andJohnson, 1979, Pisinger, 2007]. The quadratic knapsack problem has been widely studied in the literature, see Pisinger [2007] for a comprehensive survey. Exact solution algorithms are mainly based on branch-and-bound (B&B) schemes. Besides the model of QKP, Gallo et al. [1980] also presented the first B&B algorithm for this optimization problem. Their approach makes use of upper bounds that are computed by a relaxed version of QKP. The objective function is replaced by an upper plane, that is a linear function g such that g(x) ≥ x P x for any feasible solution x of QKP. Solving this new optimization problem, which is a classical knapsack problem due to the linear objective function, yields an upper bound on QKP. Billionnet and Calmels [1996] linearize the quadratic problem to compute upper bounds for a B&B algorithm. Caprara et al. [1999] use upper planes for computing upper bounds. In addition, they present a reformulation to an equivalent instance of QKP using Lagrangian relaxation. The reformulated instance provides a tight upper bound at the root node of the B&B algorithm. Billionnet et al. [1999] use upper bounds based on Lagrangian decomposition and Billionnet and Soutif [2004] introduce algorithms for fixing variables based on these bounds. The aim of these algorithms is to reduce the size of the problem before applying the B&B. Helmberg et al. [2000] apply several semidefinite relaxation techniques to QKP and include cutting planes to improve their basic approaches.  introduce an aggressive reduction algorithm for large instances of QKP, where aggressive means that a large effort is spent to fix as many of the decision variables as possible at their optimal value such that the final optimization is done rather easy. The authors apply the bounds of Caprara et al. [1999] and Billionnet et al. [1999] for the reduction of QKP. Rodrigues et al. [2012] present a linearization scheme for QKP that provides tight upper bounds for a B&B algorithm.
However, few results are known about the approximation of QKP. Since the problem is strongly N P-hard, a fully polynomial time approximation scheme (FPTAS) cannot be expected unless P = N P. Furthermore, it is unknown whether there exists an approximation with a constant approximation ratio for QKP. Taylor [2016] present an approximation algorithm based on an approach for the densest k-subgraph problem. They show that for ε > 0, QKP can be approximated with an approximation ratio in O(n 2 /5+ε ) and a run time of O(n 9 /ε ). Rader and Woeginger [2002] prove that for a variant of QKP, where positive as well as negative profit coefficients p ij are considered, there does not exist any polynomial time approximation algorithm with finite worst case guarantee unless P = N P.
Other approximation results concentrate on special cases of QKP where the underlying graph G = (V, E), with E = {(i, j) : i, j ∈ V, i = j, p ij = 0}, has a specific structure. Pferschy and Schauer [2016] present an FPTAS for QKP on graphs of bounded tree width, which includes series-parallel graphs [see Bodlaender and Koster, 2008]. Furthermore, the authors introduce a polynomial time approximation scheme (PTAS) for graphs that do not contain any fixed graph H as a minor, which includes planar graphs. As negative results, Pferschy and Schauer [2016] show that QKP on 3-book embeddable graphs is strongly N P-hard and Rader and Woeginger [2002] prove that QKP on vertex series-parallel graphs is strongly N P-hard. Kellerer and Strusevich [2010] introduce a very special variant of QKP: the symmetric quadratic knapsack problem. In addition to assigning a profit to pairs of items that both have been selected, this variant also assigns a profit to pairs of items that both have not been selected. Furthermore, the profits p ij are built as a multiplicative of two coefficients of which one also defines the weight of the constraint. Kellerer and Strusevich [2010] introduce an FPTAS to solve the problem, which is further improved by Xu [2012].

Rectangular knapsack problems
The (cardinality constrained) rectangular knapsack problem (RKP) is a variant of QKP which can be written as follows: where a, b ∈ N n , with a, b = 0 n , and k ∈ N, k < n. Note, that P = a b , i. e., rank(P ) = 1, with p ij = a i b j and p ji = a j b i , i. e., in general, P is not symmetric. We assume that k ≥ 2. Otherwise, i. e., if k = 1, the problem reduces to finding the largest coefficient a i b i , for i ∈ {1, . . . , n}.
The rectangular objective function is formulated in analogy to the so-called Koopmans-Beckmann form of the quadratic assignment problem, see Burkard et al. [1998], which is also a particular case of the more general Lawler formulation. In both cases, the two respective four dimensional arrays of profit/cost coefficients are given as a product of two lower dimensional parameter matrices or vectors, respectively.

Illustrative interpretation
The denotation rectangular knapsack problem is motivated by the special structure of P given by the coefficients a i b j . Each coefficient can be interpreted as the area of a rectangle. Accordingly, for fixed itemî ∈ {1, . . . , n}, all rectangles corresponding to coefficients aîb j , j = 1, . . . , n, have the same width, and all rectangles corresponding to coefficients a j bî, j = 1, . . . , n, have the same height. Note that the objective function can be rewritten as B. Schulze et al.
which can be interpreted as choosing a subset S ⊂ {1, . . . , n} of items such that the area of the rectangle with width i∈S a i and height i∈S b i is maximized.
Example 1 We consider the following instance of RKP: max (4, 5, 2, 12, 7) x · (6, 3, 8, 5, 10) x The corresponding rectangles are plotted in Figure 1. Each rectangle has the same position in the overall rectangle as the corresponding coefficient p ij = a i b j in the profit matrix P . The optimal solution x = (0, 1, 0, 0, 1) generates an objective function value that corresponds to the highlighted area in the figure.

Bounds
The structure of the profit matrix P implies an easy computation of bounds for RKP. In the following, we assume that all instances are defined or reordered such that a 1 ≥ . . . ≥ a n and, in case of ties, i. e., if Let S n denote the symmetric group of order n and π ∈ S n denote a permutation of {1, . . . , n}. More specifically, consider π such that b π(1) ≥ . . . ≥ b π(n) and in case of ties, i. e., if b π(j) = b π(j+1) for j ∈ {1, . . . , n − 1}, such that a π(j) ≥ a π(j+1) .
Using the sorted coefficients a i , b π(j) of the objective function, one can compute an upper bound for RKP in a straightforward way.
Lemma 1 For every feasible solution x ∈ {0, 1} n of RKP, the following inequality holds: (1) Note that, in general, this upper bound does not correspond to a solution of RKP since the value of a variable x i may be differently defined w. r. t. the respective sorting of the coefficients. As soon as Equation (1) holds, the upper bound U corresponds to a feasible solution of RKP and this solution is optimal.
Proof We consider the objective function of RKP: The cardinality constraint restricts the number of selected items to k. Due to the ordering of coefficients a i , it is for every feasible solution x of RKP. Analogously, due to the definition of the permutation π, we know that for every feasible solution x of RKP. Thus,

Preprint -Preprint -Preprint -Preprint -Preprint -Preprint
. . , k}, the upper bound is based on the selection of the k items 1, . . . , k: The solution x is feasible and realizes U. Hence, x is optimal and U is a tight upper bound.
A lower bound on RKP can also be obtained by using the sorting of the coefficients. Letx andx ∈ {0, 1} n be defined as follows: For notational convenience, let κ . . = k 2 , κ . . = k 2 , and κ . . = k 2 . If k is even, the equality κ = κ = k 2 holds, i. e.,x andx are identical. Remark 1 The definition ofx guarantees that at least the product is realized in the objective function. Due to the ordering of the coefficients a i and b π(j) , this is the maximal possible value that a product of κ coefficients a i and κ coefficients b j can achieve. The same holds analogously forx. This property is important to prove an approximation quality in the following, see the proof of Theorem 1.
We define L . . = f (x) and L . . = f (x). The following example shows a connection between the bound computation and the visualization of RKP as a selection of a subset of rectangular areas.
In this context, we show that the following inequality holds: Referring to the description of Example 2, the right-hand side of this inequality corresponds to the area of the κ · κ largest rectangles in the left upper part of the overall rectangle (see also Remark 1). We partition the area corresponding Fig. 3 Area that defines the lower bound L = L ( ) for Example 2. The assignment of labels I to IV is relevant for the proof of Equality (5).
to the lower bound into four distinct areas to show that the inequality holds: In the context of Example 2, the four terms resulting from this partition correspond, in this order, to the four areas (I to IV) in Figure 3. Analogously, using the definition ofx, it holds that: 12 B. Schulze et al.

Approximation algorithms
The results of Section 3.2 naturally motivate an approximation algorithm, see Algorithm 1. It computes the solutionsx andx and outputs the better alternative as an approximate solution.
Proof Algorithm 1 returns a feasible solution in polynomial time O(n log n).

Case 1: k even
Since the coefficients a i , b π(j) are in non-increasing order, it holds that

Preprint -Preprint -Preprint -Preprint -Preprint -Preprint
On the Rectangular Knapsack Problem 13 Case 2: k odd In analogy to case 1 we again use the fact that the coefficients a i , b π(j) are in non-increasing order. We can assume without loss of generality that: This inequality is equivalent to: Thus, the following inequality holds. Note that we use Equations 5 and 6 to bound several terms. In summary, this yields the approximation factor:

Preprint -Preprint -Preprint -Preprint -Preprint -Preprint
As presented in the proof of Theorem 1, we can guarantee better results for even values of k. Also, if k is odd the quality of the approximation increases for increasing values of k.

Remark 3
-If k is even, the result of Theorem 1 improves to a 4-approximation algorithm. -For fixed odd values of k, Algorithm 1 is a polynomial time ρ-approximation algorithm for RKP with (cf. Equation (8) However, in the worst case the approximation ratio is tight as is shown in the following example.
Thus, for increasing values of M the approximation ratio tends towards for even values of k and lim M →∞ for odd values of k ≤ 3, respectively. Note that, for fixed values of k, ρ odd exactly matches the approximation ratios given in Remark 3.

Improvements of the approximation algorithm
In practice, Algorithm 1 can be improved in two different ways. A first observation is that, due to the definition of the lower bound solutionx (c.f. (2)), we do not use the full capacity of RKP, if the sets {1, . . . , κ} and {π (1), . . . , π(κ)} are not disjoint, i. e., if n i=1x i < k. Hence, it is possible to increase the lower bound value by including further items. Algorithm 2 demonstrates a possible procedure to compute an improved lower bound L impr that takes this into account.
An additional parameter k , which we call adaptive capacity, is introduced to increase the sets {1, . . . , κ} and {π(1), . . . , π(κ)}, and, therefore, increase the number of selected items, without violating the constraint. In the beginning, k is set to k. After computing the lower bound solutionx as defined in (2), the algorithm tests whether k items are selected or not. In the latter case, the adaptive capacity k is increased by the difference k − n i=1x i . A re-computation ofx, using k as capacity, allows to include more items in accordance with the ordering of the respective coefficients a i or b π(i) which compensates for the fact that the original sets are not disjoint. Subsequently, it is tested again if the constraint is satisfied with equality. If not, the adaptive capacity k is further increased. Otherwise, the algorithm continues by computingx using the current value of the parameter k as capacity and testing which of the lower bound values is larger.
Lemma 3 If the solutionx allows to increase the adaptive capacity k to k + (k − n i=1x i ) (in Step 10 of Algorithm 2), then this increase is also feasible for the computation ofx. Proof For ease of notation, we assume that we are examining the iteration where the adaptive capacity k is increased for the first time from the capacity k to k = k + (k − n i=1x i ). The following discussion can be applied in an 16 B. Schulze et al.
If k is even, we know thatx =x and the statement is trivially true. Otherwise, i. e., if k is odd, we can take advantage of the fact that the solutionx orx uses less than k items if: It holds that n i=1x i = |Ĩ| and that n i=1x i = |Î|. Furthermore, we know that (1), . . . , π(κ)}, i, e., ifĨ = J |J | + 1 , else .
Lemma 4 Let n be the number of items and let k be the capacity of RKP. Algorithm 2 terminates, has a worst case time complexity of O(n log n), and a worst case approximation ratio of 4.5.
Proof The while-loop for computing the solutionx with n i=1x i = k is critical for the termination of Algorithm 2. In the first iteration, at least κ variables are set to 1. The parameter k is increased by at least 1 in each consecutive iteration and, thus, in at least every second iteration an additional entry ofx B. Schulze et al.
is set to 1. Hence, after at most 2 · κ + 1 = k iterations k variables have been selected forx and the loop terminates.
We take advantage of the ordering of the coefficients to set only new variables to 1 if the adaptive capacity is increased. Thus, the execution of the while loop requires O(k). The complexity of Algorithm 2 is determined by the sorting algorithm (cf. Algorithm 1), i. e., Algorithm 2 has a worst case time complexity of O(n log n).
The approximation ratio is at most 4.5, since Algorithm 2 computes the same optimal objective function value for the RKP instance of Example 3 as Algorithm 1.
Example 4 We apply the improved approximation algorithm, Algorithm 2, on the instance of RKP of Example 2. The solutionx is defined by the set with |Ĩ| = n i=1x i = 4 < 5. Thus, the adaptive capacity can be set to k . . = 5 + (5 − 4) = 6. The re-computation ofx leads tõ with |Ĩ| = n i=1x i = 5. Hence, the cardinality constraint is tight and, since k is even, the solution x = (1, 1, 1, 1, 0, 1, 0, 0) generates the improved lower bound L = f (x) = 920. The optimal solution of the instance x * is identical to the improved lower bound solution x * = x.
As proven above, setting the adaptive capacity to k . . = 6 is also feasible forx. The solutionx defined by k = 5 corresponds to the set with |Î| = n i=1x i = 4 < 5. Hence, one additional item can be included, resulting again in the same lower bound solution x = (1, 1, 1, 1, 0, 1, 0, 0) (with k = 6). The areas of rectangles corresponding to the lower bounds L and L based on the first computations ofx andx (Algorithm 1), respectively, and the improved lower bound L (Algorithm 2) are shown in Figure 4.
A second improvement of Algorithm 1 is motivated differently: Without an analysis of the input data, the distribution of the entries of the coefficient vectors a and b is unknown, and, thus, there might be better selections then deciding equally according to both of the orderings. One possible approach is to compute various alternative solutions, still based on the sorting of the coefficients, and select the best solution. The alternatives can be defined by setting the variables only corresponding to the sorting of a, to the sorting of b, and by all alternatives in between, i. e., x 1 = . . . = x j = 1, x π(1) = . . . = x π(k−j) = 1 and x i = 0 for all remaining indices, for 0 ≤ j ≤ k.
This approach is formalized in Algorithm 3. Following the idea of Algorithm 2, it includes a test if k items have been selected in the current solution. If not, further items are included according to the ordering of the coefficients b π(i) . Other strategies are possible: include further items according to the ordering of the coefficients a i , alternate between both orderings or include arbitrarily chosen items. The quality of those strategies strongly depends on the given problem instance.
In practice, it is quite intuitive to assume that Algorithm 3 leads to better approximation results than that of the basic version of Algorithm 1. However, the theoretical approximation ratio is the same.
Theorem 2 Algorithm 3 is a polynomial time 4.5-approximation algorithm for the rectangular knapsack problem.
Proof Algorithm 3 returns a feasible solution in polynomial time, since each alternative solutionx and the corresponding objective function value can be computed in linear time and there are linearly many alternatives (c. f. Theorem 1).
The approximation ratio is at least 4.5, since the solutionsx andx of Algorithm 1 are included in the set of alternatives. The approximation ratio is at most 4.5, since Algorithm 3 computes the same optimal objective function value for the RKP instance of Example 3 as Algorithm 1.
Three different classes of instances were generated to test the algorithms: Uncorrelated instances The coefficients a i , b i are generated according to a uniform distribution within the range [0, 100].
Positively correlated instances The coefficients a i are generated according to a uniform distribution within the range [0, 100] and b i = a i +n(i) where n(i) is a value generated according to a uniform distribution within the range [−5, 5].
Negatively correlated instances The coefficients a i are generated according to a uniform distribution within the range [0, 100] and b i = max{100 − a i + n(i), 0} where n(i) is a value generated according to a uniform distribution within the range [−5, 5].
For each type of instances, four different constraint slacknesses c k , with k = c k · n , were chosen: c k = 0.1, c k = 0.25, c k = 0.5, and c k = 0.75. The instance sizes were n = 100, 200, 300, 400, except for the negatively correlated instances, with n = 25, 50, 75. For the latter instance class, the QKP solver was not able to solve instances with n ≥ 75 and k ≥ 14 within one hour of CPUtime. For each combination of instance class, size and constraint slackness, 10 instances were generated. Noteworthy, all approximation algorithms required at most 0.01 seconds for all instances tested.  Tables 1 to 3 present the average results obtained for the three classes of instances where columns z * /L• refer to the average approximation ratios obtained by the four algorithm variants and column U/L * gives an upper bound on the approximation ratio. In general, the results indicate that the approximation quality of all algorithm variants is much better than the guaranteed approximation ratio of 4.5 and that the improved versions yield even better  Table 3 Results for negatively correlated instances.
results than the basic variant, except on negatively correlated instances, for which all versions presented a similar performance. Moreover, the instance size does not play a strong role on the approximation ratio. However, the performance of the four variants seem to be affected by the instance type. In the following, we discuss the results in more detail for each instance type.

Uncorrelated instances
The experimental results in Table 1 suggest that the improved variant performs better as the constraint slackness increases and that a larger capacity value k improves the approximation (see Remark 3). Differently, the basic variant does not seem to be affected by c k and presents the worst approximation ratio in all cases. The shifted and combined variants present the best approximation ratio for small c k . Both variants also improve the approximation for larger c k but not as much as for the improved variant, which gives the best approximation ratio.
Positively correlated instances Table 2 shows that the basic variant has the worst approximation ratio, although still far from the theoretical bound. For this variant, many items seem to be selected twice due to both orderings, which can be expected for positive correlated instances since the orderings of the coefficients to both objective functions should be rather similar. Noteworthy, all the three improved variants are close to an approximation ratio of 1.0. We observed that as soon as the equality in the constraint is ensured, all improved variants solve the problem to optimality.
Negatively correlated instances For this instance type, all variants present a similar approximation ratio; see Table 3. In fact, the basic variant generates very good approximation results with approximation ratios close to 1.0 for the tested instances. This is especially good, since the exact algorithm could not solve larger instances in less then one hour of CPU-time whereas the approximation algorithm can compute a high quality approximation in less then 0.01 seconds.

Conclusion
We presented a geometric interpretation of the rectangular knapsack problem. Upper and lower bounds for the problem can be computed directly by sorting the coefficients of the objective function. Based on these bound computations, we introduced a polynomial time approximation algorithm for RKP that provides an approximation ratio of 4.5. In practice, however, the algorithm can be further improved by selecting additional items if the cardinality constraint is not met with equality. Furthermore, the selection strategy for items can be modified We tested all algorithm variants on knapsack instances with three different correlation structures, up to 400 items, and four different constraint slacknesses. The approximations were computed in 0.01 seconds or less per instance. We observed that in practice the approximation ratios of all algorithms are much better than the theoretical ratio of 4.5. Thus, our approximation algorithms are an efficient tool to compute approximations of good quality for RKP.
In the future it would be interesting to integrate the bound computations in a branch-and-bound procedure to formulate an exact algorithm for RKP. Furthermore, the results seem to be transferable to higher dimensions, where we think of problems of the form x i ≤ k x i ∈ {0, 1}, i = 1, . . . , n.
The bound computations and algorithm formulations should be convertible without problems, whereas the proof of an approximation ratio may become more complicated due to more possible cases that may occur. We also suggested a field of application for RKP. Finding a representative solution of the bi-objective cardinality constrained knapsack problem that maximizes the hypervolume with the origin as reference point is modeled by the rectangular knapsack problem. It is, therefore, very interesting for future research.