Approximation Schemes for Binary Quadratic Programming Problems with Low cp-Rank Decompositions
aa r X i v : . [ c s . D S ] N ov Approximation Schemes for Binary Quadratic ProgrammingProblems with Low cp-Rank Decompositions
Khaled Elbassioni ∗ Trung Thanh Nguyen † November 8, 2018
Abstract
Binary quadratic programming problems have attracted much attention in the last few decades due totheir potential applications. This type of problems are NP-hard in general, and still considered a challengein the design of efficient approximation algorithms for their solutions. The purpose of this paper is toinvestigate the approximability for a class of such problems where the constraint matrices are completelypositive and have low cp-rank . In the first part of the paper, we show that a completely positive rationalfactorization of such matrices can be computed in polynomial time, within any desired accuracy. We nextconsider binary quadratic programming problems of the following form: Given matrices Q , . . . , Q n ∈ R n × n + , and a system of m constrains x T Q i x ≤ C i ( x T Q i x ≥ C i ), i = 1 , . . . , m , we seek to find a vector x ∗ ∈ { , } n that maximizes (minimizes) a given function f . This class of problems generalizes manyfundamental problems in discrete optimization such as packing and covering integer programs/knapsackproblems, quadratic knapsack problems, submodular maximization, etc. We consider the case when m and the cp-ranks of the matrices Q i are bounded by a constant.Our approximation results for the maximization problem are as follows. For the case when the objec-tive function is nonnegative submodular, we give an (1 / − ǫ ) -approximation algorithm, for any ǫ > ;when the function f is linear, we present a PTAS. We next extend our PTAS result to a wider class ofnon-linear objective functions including quadratic functions, multiplicative functions, and sum-of-ratiofunctions. The minimization problem seems to be much harder due to the fact that the relaxation is not convex. For this case, we give a QPTAS for m = 1 . ∗ Department of Electrical Engineering and Computer Science, Masdar Institute of Science and Technology, Abu Dhabi, UAE † Department of Electrical Engineering and Computer Science, Masdar Institute of Science and Technology, Abu Dhabi, UAE Introduction
Binary quadratic programming is a classical class of combinatorial optimization problems containing a quadraticfunction either in its objective or in its constraints. The study of this problem has been of great theoreticalinterest in the last few decades due to its wide range of applications, for example, capital budgeting and fi-nancial analysis [Lau70, MY80], machine scheduling [AKA94], and traffic message management problems[GHS80]. One can mainly distinguish two types of special classes of binary quadratic programming that havebeen studied extensively in the literature. The first class involves maximizing a quadratic function without anystructural constraints (a.k.a unconstrained - quadratic problem ) which is is known to be APX-hard sinceit generalizes, for example, the maximum cut problem. The second class deals with the quadratic knapsackproblem (QKP) arising from the classical knapsack problem by maximizing a quadratic objective functionsubject to (linear non-negative) knapsack constraint(s) (see. e.g., [Pis07]). Very recently, Yu and Chau [YC13]introduced a new variant of the knapsack problem, called the complex-demand knapsack problem (CKP), thatarose from the allocation of power in AC (alternating current) electrical systems. In this setting, the constraintis exactly a sum of two squares of linear non-negative functions, and thus is a generalization of the classicalknapsack constraint.Binary programming with quadratic objective function and/or quadratic constraints is hard to approximatein general. For example, quadratic knapsack problem is NP-hard to approximate to within any finite worstcase factor [JW02]. Most algorithms for this kind of problems focus on handling various special cases of thegeneral problem [JW02, PS13, KS10, KS12, Xu12]. Those papers investigate the existence of approximationschemes by exploiting the graph-theoretic structure of QKPs or taking into account the special multiplicativestructure of the coefficients of the objective function. Another remarkable direction in studying the approx-imability of binary quadratic programming is to restrict the objective function to certain special types, inparticular, those having low-rank [KN07]. Intuitively, a function of low-rank can be expressed as a combina-tion of a (fixed) number of linear functions. Allemand et al. [AFLS01] study the unconstrained 0-1 quadraticmaximization problem and propose a polynomial-time algorithm when the symmetric matrix describing thecoefficients of the quadratic objective function is positive semidefinite and has fixed rank . Kern and Woegin-ger [KW07] and Goyal et al. [GKR11] give an FPTAS for minimizing the product of two non-negative linearfunctions when the convex hull of feasible integer solutions is given in terms of linear inequalities. Mittal andSchulz [MS13a, MS13b] and Goyal and Ravi [GR13] extend this result to the more general class of low-rankfunctions.In this paper, we consider a class of 0-1 quadratic programming, where a (nonnegative) linear or sub-modular objective function has to be minimized or maximized subject to a system of quadratic inequalitiesinvolving nonnegative semidefinite matrices with low completely positive ranks . Intuitively, the completelypositive rank (cp-rank) of a semidefinite matrix Q is the smallest positive integer r such that Q can be de-composed into a sum of r non-negative rank-one matrices. As a consequence, every quadratic constraint in-volving such matrices can be written as the sum of r squares of (non-negative) linear functions which we call cp-decomposition . Hence, our problem can be seen as a generalization of the submodular maximization with(linear) knapsack constraints which has been studied extensively in the literature [Svi04, KST09, LMNS09].Herein, we would like to emphasize that the cp-decomposition plays an important role not only in binaryquadratic programming [Bur09], but also in in many other areas such as block designs in combinatorial anal-ysis [Hal86], data mining [ZS05, HXZ +
11] and economics [VHVD08]. Several attempts have been devotedin the last decades to settling the complexity of determining if a given semidefinite positive matrix has fi-nite cp-rank. Surprisingly, however, this has been settled only very recently in [DG14], where it was shownthat this problem is NP-hard. In this paper we restrict out attention to the case of low cp-rank and study thepolynomial-time algorithm for finding a cp-decomposition of a matrix whose cp-rank is fixed. A function f : R n → R is said to be of rank k , if there exists a continuous function g : R k → R and linearly independent vectors a , . . . , a k ∈ R n such that f ( x ) = g ( a T x, . . . , a Tk x ) , (see [KN07]). Note that in this case the objective function can be written as a sum of the squares of not more than k linear functions, where k is the rank of the coefficient matrix. .1 Our Contribution. We summarize our contributions in this paper as follows: • We give an n O ( r ) deterministic algorithm that given a positive semidefinite matrix Q of size n decidesif a nonnegative factorization of inner dimension at most r exists. Furthermore, a rational approximationto the solution to within an additive accuracy of δ can be computed in time poly ( L , n O ( r ) , log δ ) , where L is the bit length of the input. • We present a PTAS for maximizing a linear function subject to a fixed number of quadratic packingconstraints involving positive semidefinite matrices with low cp-ranks. We provide an (1 − ǫ )( − ǫ ) -approximation algorithm for the problem of maximizing a submodular function, for any ǫ > . For thecase when the submodular function is monotone, the approximation factor is (1 − ǫ )(1 − e − ǫ ) , forany ǫ > . Our results generalize the ones of Sviridenko [Svi04] and of Kulik et al. [KST09]. • We extend our PTAS result to a wider class of non-linear objective functions including quadratic func-tions, multiplicative functions, and sum-of-ratio functions. • We give a quasi-PTAS for minimizing a linear function subject to one quadratic covering constraintinvolving a positive semidefinite matrix with low cp-rank.
The problem of maximizing a nonnegative submodular function under knapsack constraints has been studiedextensively in the literature. Sviridenko [Svi04] studied the problem for monotone submodular functions andgave a (1 − /e ) -approximation for the case of one knapsack constraint. Kulik et al. [KST09] consideredthe problem with any constant number of knapsack constraints, and gave a (1 − /e − ǫ ) -approximation, forany ǫ > . Lee et al. [LMNS09] investigated the problem of maximizing a general submodular functionunder a constant number of knapsack constraints, and presented algorithms that achieve approximation ratioof / − ǫ . A better approximation factor of / − ǫ was obtained in [KST09]. In this paper, we will extendthese results to the case with quadratic knapsack constraints having low cp-rank.Chau et al. [CEK14] considered the problem CKP with linear objective function and gave a PTAS , whichsettles the complexity of the problem given the strong NP-hardness in [YC13]. As CKP is a special case of theproblem of maximizing a linear function subject to a single quadratic inequality whose constraint matrix hascp-rank 2, our second result above can be thought of as a generalization of the PTAS in [CEK14] in severaldirections: we allow the cp-rank to be any constant, we consider any constant number of inequalities, andconsider more general objective functions. In fact, our result for submodular functions is based essentially onthe same geometric idea used in [CEK14]. A positive semi-definite matrix Q (cid:23) is said to have completely positive rank r , denoted cp - rank( Q ) = r ,if r is the smallest number of non-negative rank-one matrices into which the matrix Q can be decomposedadditively, i.e., Q = P rj =1 q j ( q j ) T , where q j ∈ R n + are non-negative vectors. If no such r exists then cp - rank( Q ) = + ∞ . A notion related to the cp-rank is the nonnegative rank of a nonnegative (not necessarilysymmetric) matrix: a positive integer r is called the nonnegative rank of a matrix P ∈ R n × m + if and only ifit is the smallest number such that there exist two matrices U ∈ R n × r + , V ∈ R r × m + such that P = U V . Itis worth noting that, in general, the nonnegative rank and the cp-rank can be different from each other. Infact, for a completely positive matrix Q , it is known that the nonnegative rank is always less than or equalto the cp-rank. We refer to the textbook [BSM03] and the references therein for more discussion about thecp-rank. The complexity of deciding the existence of such a cp-decomposition for a given matrix remainedopen for a long time. Very recently, Dickinson and Gijben [DG14] proved that this problem is NP-hard but2eft open the question whether or not it belongs to NP. The similar hardness result was also proved for thecase of nonengative rank by Vavasis [Vav09]. Therefore, it is natural to pay attention to special cases in whichthe factorization problem can be solved efficiently. Arora et al. [AGKM12] and Moitra [Moi13] proposedexact polynomial-time algorithms for computing a nonnegative factorization when the nonnegative rank ofthe input matrix is constant. Their basic idea is to transform the factorization problem into a problem offinding a solution of a system of polynomial equations which is known to be solved by an algorithm (e.g,[Ren92, BPR96]) whose running-time is polynomial in the number of equations but is exponential in thenumber of variables. Based on this idea, Arora et al. [AGKM12] gave a transformation to a system with nm equations and r r variables. As a result, checking if the nonnegative rank of P is at most r can be donein O (( nm ) r r ) time. This has then been improved by Moitra [Moi13] who gave a method to exponentiallyreduce the number of variables to r , thus yielding a faster algorithm with running-time ( nm ) O ( r ) . It wasalso shown by Arora et al. [AGKM12] that achieving an exact ( nm ) o ( r ) -algorithm is impossible under theExponential Time Hypothesis [IP01]. In the following, we will show that one can obtain an exact algorithmfor the cp-rank case by employing the same idea as for the nonnegative rank. Assume Q is rational and let L be the maximum bit length of any element in Q . The result is stated in Theorem 2.1 below. Theorem 2.1.
There is an n O ( r ) time deterministic algorithm that given a positive semi-definite matrix Q ofsize n produces a nonnegative factorization Q = U U T of inner dimension at most r if such a factorizationexists. Furthermore, a rational approximation to the solution to within an additive accuracy of δ can becomputed in time poly( L , n O ( r ) , log δ ) .Proof. We proceed essentially along the same lines as in the proof of Theorem 2.1 in [AGKM12] (for findinga simplicial factorization). Let Q be a nonnegative positive semi-definite matrix of size n . Assume that Q hasa decomposition Q = U U T , where U is a nonnegative matrix of size n × r . A basic fact from linear algebrais that Q, U and U T have the same rank: rank( Q ) = rank( U ) = rank( U T ) = s ≤ r .Fix an arbitrary basis V of the column vectors of Q , let B be an n × s matrix corresponding to this basisand let Q V be the matrix of size s × n corresponding to the (unique) representation of Q in the basis V , thatis, BQ V = Q . To prove the claim of the theorem, it will suffice to prove the following. Claim 2.2. Q has a non-negative factorization Q = U U T of inner-dimension at most r if and only if there isan r ′ × s matrix H , with r ′ ≤ r , satisfying the following two conditions:(i) HQ V is a nonnegative matrix,(ii) ( HQ V ) T ( HQ V ) = Q .Proof. ( ⇐ ) Suppose that the conditions (i) and (ii) are satisfied. Let U = ( HQ V ) T . This matrix is nonnega-tive and has size of n × r and thus, Q has cp-rank at most r ′ ≤ r . ( ⇒ ) Now suppose that there is a non-negative factorization Q = U U T , where U is an n × r ′ -matrix, with r ′ ≤ r . The singular value decomposition of U has the form U = LSR T , where L, R are orthogonal matrices( L − = L T and R − = R T ) of size n × n and r ′ × r ′ , respectively, and S is an n × r ′ diagonal matrix. Thefirst s (non-zero) diagonal entries of S are exactly the singular values of U (corresponding to the first s (non-zero) diagonal entries of Q ). We can rewrite Q in the form Q = ( LSR T )( LSR T ) T = ( LSR T )( RS T L T ) .Let N = RS + L T (an r ′ × n matrix), where S + is the pseudo-inverse of matrix S . More precisely, S + is an r ′ × n matrix such that: SS + = (cid:18) I s × s
00 0 (cid:19) n × n and S + S = (cid:18) I s × s
00 0 (cid:19) r ′ × r ′ , where I s × s is an identity matrix of size s .Define an r ′ × s matrix H = N B . We now prove that this matrix T satisfies the two conditions (i) and(ii). Indeed, we have: HQ V = N BQ V = ( RS + L T ) BQ V = ( RS + L T ) Q = ( RS + L T )( LSR T )( RS T L T )= RS + SS T L T = RS T L T = U T , cp - rank( Q ) ≤ r is equivalent to determining if a system of at most nr linearinequalities and n quadratic equations on at most r variables is feasible. This decision problem can be solvedin n O ( r ) time using quantifier elimination algorithms [BPR96]. Furthermore, a rational approximation to thesolution to within an additive accuracy of δ can be computed in time poly( L , n O ( r ) , log δ ) ; see [GV88,Ren92]. Corollary 2.3.
For a real matrix M let k M k ∞ := max i,j | M ij | . Given a rational positive semi-definite matrix Q of size n such that cp - rank( Q ) = r , and ǫ > , one can compute a rational nonnegative n × r -matrix e U such that e U e U T ≥ Q and k Q − e U e U T k ∞ ≤ ǫ , in time poly( L , n O ( r ) , log ǫ ) .Proof. Let s = rank( Q ) and B and Q V be the matrices corresponding to basis V of Q , as in the proof ofTheorem 2.1. Then from this proof, it follows that we can discover in n O ( r ) time that there is an r × s -realmatrix H such that conditions (i) and (ii) of claim 2.2 hold. Furthermore, in time poly( L , n O ( r ) , log δ ) ,we can find a rational matrix b H , such that k H − b H k ∞ ≤ δ , for any desired accuracy δ > . Let B ′ be anon-singular s × s -submatrix of B , obtained by selecting s linearly independent rows of B , and let Q ′ be thecorresponding s × n submatrix of Q . Then Q V = ( B ′ ) − Q ′ . Let us choose δ := min ( ǫ k ( B ′ ) − k ∞ k Q k / ∞ rs (4 + 3 s k Q k / ∞ ) , · max { s k ( B ′ ) − k ∞ , } ) , (1)and set δ ′ := sδ k ( B ′ ) − k ∞ and e H := b H + δ ′ E r × s B ′ , where E r × s is the r × s -matrix with all-ones. Define U := ( HQ V ) T , b U := ( b HQ V ) T , and e U := ( e H Q V ) T . Then, it follows that e U T = U T + [ δ ′ E r × s + ( e H − b H )( B ′ ) − ] Q ′ ≥ U T + ( δ ′ − sδ k ( B ′ ) − k ∞ ) E r × s Q ′ = U T , where the first inequality follows from the fact that Q ′ ≥ , and the second equality follows from our choiceof δ ′ . Let us further note that e U e U T − U U T = ( Q ′ ) T ∆ T U T + U ∆ Q ′ + ( Q ′ ) T ∆ T ∆ Q ′ , (2)where ∆ := δ ′ E r × s + ( e H − b H )( B ′ ) − . To bound k Q − e U e U T k ∞ , we bound the norm of each term in (2).We first note that k ∆ k ∞ ≤ δ ′ + δs k ( B ′ ) − k ∞ = 2 δ ′ , and thus k U ∆ Q ′ k ∞ ≤ δ ′ k U E r × s Q ′ k ∞ ≤ δ ′ rs k U k ∞ k Q k ∞ ≤ δ ′ rs k Q k / ∞ , (3)where the last inequality follows from the fact that k U k ∞ ≤ k Q k / ∞ . Now we write ∆ T ∆ = ( δ ′ ) E s × r E r × s + δ ′ ( E s × r X + X T E r × s ) + X T X, (4)where X := ( e H − b H )( B ′ ) − . From (4) and k X k ∞ ≤ δ ′ follows the inequality k ∆ T ∆ k ∞ ≤ r ( δ ′ ) , andthen the inequality k ( Q ′ ) T ∆ T ∆ Q ′ k ∞ ≤ r ( δ ′ ) k ( Q ′ ) T E s × s Q ′ k ∞ ≤ r ( sδ ′ ) k Q k ∞ . (5)Using (3) and (5) in (2), we get k e U e U T − Q k ∞ ≤ δ ′ rs k Q k / ∞ + 3 r ( sδ ′ ) k Q k ∞ ≤ ǫ, by our selection (1) of δ . Finally note that log δ ≤ poly( L , r log ǫ ) .4 Approximation Algorithms for Binary Quadratic Programming
Binary programming with quadratic constraints (BQC) is the problem of maximizing (or minimizing) a func-tion of a set of binary variables, subject to a finite set of quadratic constraints. A BQC problem takes theform:(P
ACKING -BQC) max f ( x ) (6)subject to x T Q i x ≤ C i , i ∈ [ m ] x ∈ { , } n , (C OVERING -BQC) min f ( x ) (7)subject to x T Q i x ≥ C i , i ∈ [ m ] x ∈ { , } n , where f is a nonnegative function, and Q i is (without loss of generality) an n × n nonnegative symmetricmatrix, for all i ∈ [ m ] := { , , . . . , m } . For simplicity, we fix some useful notation. For any integer n ≥ ,we denote by CP n the set of completely positive matrices (i.e., with finite cp-rank), and by CP ∗ n the subsetof CP n that includes matrices of fixed cp-rank. Furthermore, we identify a vector x ∈ { , } n with a subset S ⊆ [ n ] , i.e, write S = S ( x ) = { i ∈ [ n ] | x i = 1 } . Hence, for a function f defined on the power set [ n ] , f ( x ) ≡ f ( S ) . We shall consider non-negative functions f with a special structure: f is linear if f ( x ) := u T x for some vector u ∈ R n + ; f is submodular if f ( S ∪ T ) + f ( S ∩ T ) ≤ f ( S ) + f ( T ) for all subsets S, T ⊆ [ n ] ; f is quadratic if f ( x ) = x T Qx + u T x , for some Q ∈ CP ∗ n and u ∈ R n + ; we refer to the correspondingpacking (covering) problems as P ACK -L IN , P ACK -S UB , and P ACK -Q UAD (C OVER -L IN , C OVER -S UB , andC OVER -Q UAD -QC), respectively. As standard in the literature, we assume that the submodular function f isgiven by a value oracle that always returns the value f ( S ) for any given set S ⊆ [ n ] .For α > , a vector x ∈ { , } n (or the corresponding set S ⊆ [ n ] ) is said to be α -approximate solution for problem P ACKING -BQC (resp., C
OVERING -BQC), if x is a feasible solution satisfying f ( x ) ≥ α · O PT (resp., f ( x ) ≤ α · O PT ), where O PT is the value of an optimal solution; for simplicity, if α = 1 − ǫ (resp., ǫ ), for ǫ > , x will be called ǫ -optimal .It can be seen that the problem (BQC) is NP-hard as it includes among others the multi-dimensionalKnapsack packing and covering problems as special cases, when we set f to be a linear function and set eachmatrix Q i to be a rank-one non-negative matrix. A fully polynomial-time approximation scheme (FPTAS) for the knapsack problem was given in [IK75], and a PTAS for the m -dimensional knapsack problem wasgiven in [CHW76, FC84], and these results are best possible, assuming P = NP. When f is submodular, it wasshown that there exist constant factor approximation algorithms for the packing problem, subject to a constantnumber of knapsack constraints (e.g., [Svi04, KST09]). We extend these results as follows. Theorem 3.1. P ACK -L IN admits a PTAS when m and each cp - rank( Q i ) is fixed. Theorem 3.2.
For any ǫ > , there is a (1 − ǫ ) α -approximation algorithm for P ACK -S UB when m and each cp - rank( Q i ) is fixed, where α = (1 − e − ǫ ) , if the objective f is monotone, and α = − ǫ , otherwise. Theorem 3.3. C OVER -L IN admits a QPTAS when m = 1 , and cp - rank( Q ) is fixed. The next sections are devoted to present two methods for designing approximation algorithms for P
ACKING -BQC. The first method relies on the polynomial solvability of the convex relaxation, whereas the second onetakes into account the geometric properties of the region of feasible solutions. The first method works for thecase of linear objective function and allows also for an additive linear term in each quadratic constraint; itcan also be extended to handle more general classes of objective functions such as quadratic functions withfixed cp-rank, multiplicative functions, and sums of ratios of linear functions, as well as to packing constraints Note that a quadratic function f : 2 [ n ] → R + is supermodular , that is, f ( S ∪ T ) + f ( S ∩ T ) ≥ f ( S ) + f ( T ) for all subsets S, T ⊆ [ n ] . A PTAS is an algorithm that runs in time polynomial in the input size n , for every fixed ǫ , and outputs an ǫ -optimal solution; anFPTAS is a PTAS where the running time is polynomial in ǫ ; a QPTAS is similar to a PTAS but the running time is quasi-polynomial(i.e., of the form n polylog n ), for every fixed ǫ . ℓ p -norm for p ≥ . However, we are not able to show that this method works for submodularobjective functions. The second method works for the latter case, but it does not seem to have the flexibilityof the first method in handling more general objective functions and/or additional linear terms in the quadraticconstraints.For the minimization problem C OVER -L IN , we use the idea of geometric partitioning combined togetherwith a greedy covering approach to arrive at the result of Theorem 3.3; see Section 3.3. We prove Theorem 3.2in Section 3.2 and Theorem 3.1 for the linear objective case in Section 3.1; we extend the result to the quadraticobjective case in Section 4.As a technical remark, we note that the decomposition of each Q i as given by Corollary 2.3 involves anerror term ∆ i . To simplify the presentation, we will assume first that an exact decomposition Q i = U i U Ti isgiven, then show in Section 3.4 to see how to deal with the error term. In this section, we present a polynomial-time approximation scheme for P
ACKING -BQC with linear objectivefunction. This problem has the form (6) with f ( x ) = u T x , where Q i ∈ CP ∗ n and u ∈ R n + . For convenience,we will sometimes write u ( S ) := P k ∈ S u k , for a set S ⊆ [ n ] . The method we use here is based on solvingthe convex relaxation obtained by replacing the constraint x ∈ { , } n by the weaker one x ∈ [0 , n . Thecrucial point of our method is that an upper bound on the optimal value can be computed easily by solving(approximately) the relaxed problem. The obtained solution to the convex program defines a point in a certainpolytope, which is then rounded to an integral solution without losing much on the quality of the solution.The details of our method is described in Algorithm 1. Algorithm 1 P ACK -L IN -PTAS ( u, { Q i , C i } i ∈ [ m ] , ǫ ) Input:
Utility vector u ∈ R n + ; matrices Q i ∈ CP ∗ n , capacities C i ∈ R + , for i ∈ [ m ] ; and accuracy parameter ǫ Output: An ǫ -optimal solution S to P ACK -L IN v ← (0 , . . . , For all i ∈ [ m ] , decompose Q i as Q i = U i U Ti , where U i ∈ Q n × r i ⊲ r i = cp-rank ( Q i ) ¯ r ← P mi =1 ( r i + 1); λ ← ¯ rǫ for each subset U of [ n ] of cardinality at most λ do V ← { k ∈ [ n ] \ U | u k > min { u k ′ | k ′ ∈ U }} if T U Q i [ U ; U ] U ≤ C i for all i ∈ [ m ] then ⊲ (CP1 [ U ] ) is feasible Let x ∗ be an ǫ -optimal solution to the convex program (CP1 [ U ] ) t i ← U Ti [ ∗ ; N ] x ∗ N ; t ′ i ← T U Q i [ U ; N ] x ∗ N for all i ∈ [ m ] ⊲ t i is a vector of dimension r i Find a BFS y of the polytope P ( U ) such that u T y ≥ u T x ∗ N : x ← { ( x , . . . , x n ) | x k = ⌊ y k ⌋ for k ∈ [ n ] \ ( U ∪ V ) , x k = 1 , for k ∈ U , x k = 0 , for k ∈ V} ⊲ rounding down y k if u T x > u T v then v ← x return S ( v ) Let ǫ > be a fixed constant, and define λ = ¯ rǫ , where ¯ r = P mi =1 ( r i + 1) . The idea, extending the one in[FC84], is to try to guess λ items of highest utility in the optimal solution. This can be done by consideringall possibilities for choosing a set of cardinality at most λ . We denote by X ⊆ [ n ] the set of all such sets.Note that the size of X is bounded by O ( n λ ) and thus is polynomial in the size of input for every constant λ .For each U ∈ X , we define a set V , which contains all items that are not in the optimal solution, giventhat U is a subset of the optimal solution (these are the items k ∈ [ n ] \ U , with u k > u k ′ for some u k ′ ∈ U ).Let N := [ n ] \ ( U ∪ V ) . Using polynomial-time algorithms for convex programming (see, e.g., [NT08]), we6an find an ǫ -optimal solution to the convex relaxation (CP1 [ U ] ) defined as follows:(CP1 [ U ] ) max u T x subject to x T Q i x ≤ C i , i ∈ [ m ] ,x k = 1 for k ∈ U ; x k = 0 for k ∈ V ,x k ∈ [0 , for k ∈ N. We then define a polytope P ( U ) ⊆ [0 , N by replacing each constraint x T Q i x ≤ C i by the twoconstraints U Ti [ ∗ ; N ] y ≤ t i and T U Q i [ U ; N ] y ≤ t ′ i , for all i ∈ [ m ] , where t i := U Ti [ ∗ ; N ] x ∗ N ; t ′ i := T U Q i [ U ; N ] x ∗ N ; here, x U denotes the restriction of the vector x ∈ [0 , n to the set of components indexed by U , and Q i [ U ; N ] denotes the restriction of the matrix Q i to the columns and rows indexed by the sets U and N , respectively; similarly, U Ti [ ∗ ; N ] means the restriction of U Ti to the set of columns defined by N : P ( U ) := { y ∈ [0 , N | U Ti [ ∗ ; N ] y ≤ t i , T U Q i [ U ; N ] y ≤ t ′ i , for i ∈ [ m ] } . We can find a basic feasible solution (BFS) y in this polytope such that u T y ≥ u T x ∗ N , by solving at most n linear systems (see standard texts on Linear Programming, e.g., [Sch86]). Next, an integral solution x is obtained from y by dropping all fractional components to and setting x k ∈ { , } according to theassumption k ∈ U ∪ V . For each set U ∈ X such that (CP1 [ U ] ) is feasible, the algorithm outputs an integralsolution x . Let v be the solution of maximum value amongst all such solutions x . The algorithm produces aset S ( v ) of items that corresponds to v . Lemma 3.4.
Assume Q i ∈ CP ∗ n for all i ∈ [ m ] and m = O (1) . Then, for any fixed ǫ > , Algorithm 1 runsin polynomial time and produces an ǫ -optimal solution to the input instance.Proof. From the argument above, it can be easily seen that the running time of the Algorithm 1 is polynomialin size of the input, for any fixed constant ǫ . We now argue that the solution S returned by the algorithmis ǫ -optimal. Indeed, let S ∗ be the optimal solution to P ACK -L IN of utility u ( S ∗ ) = O PT . If λ ≥ | S ∗ | ,then S is an exact optimum solution to the P ACK -L IN instance. Suppose that λ < | S ∗ | . Let U be the set ofhighest-utility λ items in S ∗ . Then, the feasibility of S ∗ for (CP1 [ U ] ) guarantees that T U Q i [ U ; U ] U ≤ C i ,and hence an ǫ -optimal solution x ∗ for (CP1 [ U ] ) is found in step 7 of the algorithm. Assume w.l.o.g. that T U Q i [ U ; U ] U ≤ C i . From the feasibility of x ∗ for (CP1 [ U ] ) and the factorization of Q i , it follows that, forall i ∈ [ m ] , ( x ∗ N ) T U i [ N ; ∗ ] U Ti [ ∗ ; N ] x ∗ N + 2 · T U Q i [ U ; N ] x ∗ N ≤ C i − T U Q i [ U ; U ] U , implying that t i + 2 t ′ i ≤ C i − T U Q i [ U ; U ] U . (8)The definition of t i , t ′ i implies that x ∗ N ∈ P ( U ) . Thus, there is a BFS y of P ( U ) such that u T y ≥ u T x ∗ N .Let x ∈ { , } n be the rounded solution obtained from y in step 10. Since any BFS of P ( U ) has at most ¯ r fractional components and u k ≤ min k ′ ∈U u k ′ for all k ∈ N , it follows that u T x = u T x N + u T U ≥ u T y − ¯ r · u T U |U | + u T U ≥ u T x ∗ N + (1 − ¯ rλ ) u T U ≥ (1 − ǫ ) u T x ∗ ≥ (1 − ǫ ) O PT . It remains to argue that x is feasible for P ACK -L IN . Since x N ∈ P ( U ) , we have x T Q i x = x TN U i [ N ; ∗ ] U Ti [ ∗ ; N ] x N + 2 · T U Q i [ U ; N ] x N + T U Q i [ U ; U ] U ≤ t i + 2 t ′ i + T U Q i [ U ; U ] U ≤ C i , In fact, such algorithms can find a feasible solution x ∗ to the convex relaxation such that u T x ∗ ≥ opt − δ , in time polynomial inthe input size (including the bit complexity) and log δ , where opt is the value of the fractional optimal solution. We may assume that Tk Q i k ≤ C i , for all i ∈ [ m ] , where k is the k th unit vector; otherwise item k can be removed form the (P ACK -L IN ) problem.This implies that opt ≥ O PT ≥ ¯ u := max k u k . Now setting δ to ǫ · ¯ u assures that u T x ∗ ≥ (1 − ǫ ) opt . Remark . Note that the result presented in this section can be extended easily to the following problem: max u T x subject to x T Q i x + q Ti x ≤ C i , i ∈ [ m ] Ax ≤ b,x ∈ { , } n , where Q i ∈ CP ∗ n , q i ∈ R n + , for all i ∈ [ m ] ; A ∈ R d × n + , b ∈ R d + ; m, d are constants. In this section, we present a constant factor approximation algorithm for P
ACK -S UB . This problem has theform (6) where f ( x ) is a nonnegative submodular function, Q i ∈ CP ∗ n , and m is constant.If one attempts to use the technique of the previous section, (s)he runs into the difficulty of how to re-lax the objective function; the two well-known relaxations (see, e.g., [Dug09] for a recent survey) are: the Lov´asz extension which is known to be convex [GLS88], and the multi-linear extension [Von08] which is ingeneral neither convex nor concave. While it is not known how to solve the maximization problems for therelaxations corresponding to these two types of functions over a polytope in polynomial time, it is known[Von08, VCZ11] how to approximate within a constant factor the maximum of the multi-linear extensionover a packing polytope. It is not clear if this result can be extended to the case when the feasible region isdescribed by convex quadratic constraints of the form x T Q i x ≤ C i , where Q i ∈ CP ∗ n . While this remainsan interesting open question, we will rely here, instead, on a geometric-based approach that uses the results in[VCZ11, KST09] as a black-box. The key idea is to make use of the geometry of the problem to reduce it intoa multi-dimensional knapsack problem, which can be solved using enumeration and dynamic programmingfor the linear objective case, or LP-rounding for the submodular case.Let I = ( f, { Q i , C i } i ∈ [ m ] ) be an instance of the problem P ACK -S UB , and ǫ > be a fixed constant.We will construct a polynomial-time algorithm which produces a constant factor-approximate solution to theinstance I . Write Q i = U i U Ti , for some U i ∈ Q n × r i + , where r i = cp - rank( Q i ) . For k ∈ [ n ] , define the vector q ik ∈ Q r i + to be the k th column of U Ti . Since ¯ r , max i ∈ [ m ] r i is bounded by a constant, we may assumewithout loss of generality that ǫ < r . For r ∈ Z + and C ∈ R + , define B ( r, C ) , { ν ∈ R r + : k ν k ≤ C } to be the non-negative sector of the ball in R r of radius C centered at the origin. Then problem P ACK -S UB amounts to finding an S ⊆ [ n ] , maximizing f ( S ) , such that P k ∈ S q ik ∈ B ( r i , C i ) for all i ∈ [ m ] .Given a feasible set T ⊆ [ n ] (that is, k d T k ≤ C i for all i ∈ [ m ] ), we define R iT as the conic regionbounded as the following: R iT , n ν ∈ R r i + : k ν k ≤ C i , ν ≥ q iT o , (9)where q iT , P k ∈ T q ik . Write q ik , ( q ijk : j ∈ [ r i ]) ∈ R r i . Given R iT , we define an r i -dimensional grid in theregion R iT by interlacing equidistant ( r i − -dimensional parallel hyperplanes with inter-separation ǫ r i w ijT ,for each j ∈ [ r i ] , and with the j th principal axis as their normal, where w ijT , s C i − X j ′ = j ( q ij ′ T ) − q ijT (10)is the distance along the j th axis from q iT to the boundary of the ball . For j ∈ [ r i ] , let H ij be such hyperplanethat is furthest from the origin, perpendicular to the j th principal direction, and let H ij + be the half-space To maintain finite precision , we need to round down w ijT in (10) to the nearest integer multiple of D , where D is the commondenominator of all the rational number q ijk ; this should be done also for any computations involving square roots in the sequel. Thisintroduces an error that can be dealt with in a similar way as we deal with the error caused by the ∆ i ’s. To keep the presentationsimple, we will assume infinite precision. H ij that includes the origin. This grid defines at most r i ( r i ǫ ) r i − lines that intersect the sphericalboundary of region R iT at a set of points P iT ( ǫ ) . Let us call a cell (that, is a hypercube) of the grid boundary if it overlaps with the surface of the ball B ( r i , C i ) and denote by B iT the set of all such boundary cells. Theconvex hull of the set of points P iT ( ǫ ) defines a polytope Q iT with the following properties; the first two areeasily verifiable; the last one follows form the upper bound Theorem [McM70] and the algorithms for convexhull computation in fixed dimension (e.g. [Cha93]):(I) All points in P iT ( ǫ ) are vertices of Q iT .(II) Let H iT be the set of half-spaces including the origin and defined by the facets of Q iT . Then for each H + ∈ H iT , the facet defining H + is the convex hull of a set of vertices from P iT ( ǫ ) that are containedcompletely inside a boundary cell, and its normal is a nonnegative vector.(III) The number of facets of Q iT is at most ( r i ǫ ) r i / and they can be found in time O (( r i ǫ ) r i − ( r i log r i ǫ +( r i ǫ ) (( r i − +1) / )) .Let P iT ( ǫ ) be the polytope defined by the intersection \ H + ∈H iT H + ∩ \ j ∈ [ r i ] H ij + ∩ { ν ∈ R r i : ν ≥ } , and denote by m iT ( ǫ ) the number of facets of P iT ( ǫ ) . By construction, m iT ( ǫ ) ≤ ( r i ǫ ) r i / + 2 r i .Consider a collection of feasible subsets T , { T , . . . , T m } , T i ⊆ [ n ] , to P ACK -S UB . We define thefollowing approximate version of problem P ACK -S UB , which is obtained by approximating the ball sectors B ( r i , C i ) by inscribed polytopes.(PQP T ) max x ∈{ , } n f ( x ) subject to X k ∈ [ n ] q ik x k ∈ P iT i ( ǫ ) , for all i ∈ [ m ] . Given two vectors µ, ν ∈ R r i , we denote the (vector) projection of µ on ν by Pj ν ( µ ) , ν k ν k µ T ν . Giventhe polytope P iT ( ǫ ) , we define a set of m iT ( ǫ ) vectors { σ T,iℓ } in R r i , each of which is perpendicular to a facetof P iT ( ǫ ) , starting at the origin, and ending at the facet.For a collection of feasible subsets T , { T , . . . , T m } to P ACK -S UB , we define a submodular functionmaximization problem subject to ¯ m knapsack constraints based on { σ T i ,iℓ } , where ¯ m , P i ∈ [ m ] m iT i ( ǫ ) :( ¯ m DKS-S UB { σ T i ,iℓ } ) max x ∈{ , } n f ( x ) subject to X k ∈ [ n ] k Pj σ Ti,iℓ ( q ijk ) k x k ≤ k σ T i ,iℓ k , for all ℓ = 1 , . . . , m iT i ( ǫ ) . The following lemma follows straightforwardly from the convexity of the polytopes P iT i ( ǫ ) . Lemma 3.6.
Given a collection of feasible solutions T , { T , . . . , T m } to P ACK -S UB , problems PQP T and ¯ m DKS-S UB { σ T i ,iℓ } are equivalent. Our approximation algorithm for P
ACK -S UB is described in Algorithm 2 below, which enumerates everycollection of sets T = { T , . . . , T m } s.t. | T i | ≤ ǫ , then finds a near optimal solution for PQP T using anapproximation algorithm for ¯ m DKS-S UB (e.g., the algorithm in [KST09]).9 lgorithm 2 P ACK -S UB -A PPROX ( f, { Q i , C i } i ∈ [ m ] , ǫ, α ( ǫ )) Input:
A submodular function f : { , } n → R + ; matrices Q i ∈ CP ∗ n , capacities C i ∈ R + , for i ∈ [ m ] ;accuracy parameter ǫ ; and an α ( ǫ ) -approximation algorithm for ¯ m DKS-S UB with ¯ m = O (1) Output: (1 − ǫ ) m α ( ǫ ) -approximate solution ˆ S to P ACK -S UB For all i ∈ [ m ] , decompose Q i as Q i = U i U Ti , where U i ∈ Q n × r i + ⊲ r i = cp-rank ( Q i ) ˆ S ← ∅ for each collection of sets T = { T , . . . , T m } s.t. T i ⊆ [ n ] and | T i | ≤ ǫ do Set q iT i ← P k ∈ T i q ik , and define the corresponding vectors { σ T i ,iℓ } Obtain an α ( ǫ ) -approximate solution S to ¯ m DKS-S UB { σ T i ,iℓ } if u ( ˆ S ) < u ( S ) then ˆ S ← S return ˆ S Theorem 3.7.
For any fixed ǫ > , if the problem of maximizing a submodular function subject to multi-ple knapsack constraints can be approximated within some constant factor α ( ǫ ) , then Algorithm 2 runs inpolynomial in size of the input and finds a (1 − ǫ ) m α ( ǫ ) -approximate solution to P ACK -S UB .Proof. It is easy to see that Algorithm 2 runs in polynomial time in the number of items n . Let S ∗ be anoptimal solution. To establish the approximation ratio for Algorithm 2, we explicitly construct a set S ⊆ S ∗ in Lemma 3.8 below such that u ( S ) ≥ (1 − ǫ ) m u ( S ∗ ) , and S is feasible to PQP T for some collection T = { T , . . . , T m } s.t. T i ⊆ [ n ] and | T i | ≤ ǫ . By Lemma 3.6, invoking the α ( ǫ ) -approximation algorithmfor ¯ m DKS-S UB { σ T i ,iℓ } gives an α ( ǫ ) -approximation b S to PQP T . Then it follows that u ( b S ) ≥ α ( ǫ ) u ( S ) ≥ (1 − ǫ ) m α ( ǫ ) O PT . (11)The proof is thus completed by Lemmas 3.8-3.10 below. Lemma 3.8.
Consider a feasible solution S ∗ to P ACK -S UB . Algorithm 3 returns a subset S ⊆ S ∗ and acollection T = { T , . . . , T m } , such that (i) T i ⊆ S ∗ and | T i | ≤ ǫ ; (ii) S is a feasible solution to PQP T and u ( S ) ≥ (1 − ǫ ) m u ( S ∗ ) .Proof. Starting from S ∗ , where P k ∈ S ∗ q k ∈ B ( r , C ) , Algorithm 3 first finds sets T , S ⊆ S ∗ such that | T | ≤ ǫ , u ( S ) ≥ (1 − ǫ ) u ( S ∗ ) , and the set of vectors in S can be packed inside the polytope P T ( ǫ ) .The way for finding such sets is as follows. Let ℓ and T ℓ be the values of ℓ and T at the end of the repeat-until loop (line 8). The algorithm first constructs a nested sequence T ⊂ T ⊂ · · · ⊂ T ℓ , such that a vector q k is included in each iteration if it has a ”large” component q jk for some j ∈ [ r ] . The iteration proceeds untila sufficiently large number of vectors have been accumulated (namely, | T ℓ | ≥ ǫ ), or no vectors with largecomponents remain. At the end of the iteration, if the condition in line 9 holds, then set T = T ℓ ; otherwise,the algorithm finds a subset of S ∗ that belongs to P iT ℓ ( ǫ ) . To do so, the set S ∗ \ T ℓ of remaining vectors ispartitioned into at least ǫ − groups such that, along one of the principal axes, the total sum in each group is”large”. Then some vector or group of vectors is dropped, ensuring that the set S of remaining vectors canbe packed inside P T ( ǫ ) .We now have to consider two cases (line 12): (i) | T ℓ | ≥ ǫ , or (ii) S ℓ = ∅ . For case (i), the algorithmproceeds to line12 - combining the vectors in S \ T ℓ into one group V . For case (ii), the set S \ T ℓ can bepartitioned into at least ǫ − groups V , . . . , V h due to Lemma 3.9, where each group V p , p ∈ [ h ] , has alarge total component along one of principal axes (precisely, greater than ǫ r w jT ℓ for some j ∈ [ r ] ). Wedefine S by deleting some vector or group of vectors from S (lines 22 and 25). Hence, in both cases,10 ∗ = S k ∈ T ℓ { k } ∪ S p ∈ [ h ] V p induces a partition of S ∗ into | T ℓ | + h ≥ ǫ subsets. By Lemma 3.10 below, wehave f ( S ) ≥ (1 − | T ℓ | + h ) f ( S ∗ ) ≥ (1 − ǫ ) f ( S ∗ ) . Note that removing any one vector k ∈ T ℓ or group V p from S ∗ will insure that the resulting set S satisfies q S ∈ P T ℓ ( ǫ ) , since the lengths w jT ℓ are monotone decreasing in ℓ , for all j ∈ [ r ] , and the droppedvector or group of vectors has total length at least ǫ r i w jT along one of the principal directions j ∈ [ r ] ; onthe other hand, the boundary cell that contains the original point q S ∗ has length at most ǫ r i w jT along the j thdirection. It follows from property (II) stated above that reducing j th component of q S ∗ by ǫ r i w jT guaranteesthat the resulting vector lies in the polytope P T ℓ ( ǫ ) .Now starting from S which satisfies P k ∈ S q k ∈ B ( r , C ) , we find sets T , S ⊆ S such that | T | ≤ ǫ , f ( S ) ≥ (1 − ǫ ) f ( S ) ≥ (1 − ǫ ) f ( S ∗ ) and the set of demands in S can be packed inside the polytope P T ( ǫ ) , and so on. Finally, we find the set { T , . . . , T m , S ′ } that satisfies the lemma. The details are givenin Algorithm 3. Algorithm 3
QC-C
ONSTRUCT ( f, { q ik } k ∈ S ∗ , { C i } i ∈ [ m ] , ǫ ) Input:
A submodular function f : { , } n → R + ; vectors { q ik } k ∈ S ∗ , i ∈ [ m ] ; capacities { C i } i ∈ [ m ] ; accuracyparameter ǫ Output:
A pair ( T = { T , . . . , T m } , S ) of a collection of sets T i ⊆ S ∗ and a subset S ⊆ S ∗ S ← S ∗ for i = 1 , . . . , m do ℓ ← ; T ← ∅ repeat ⊲ Find a subset of large vectors T w.r.t. to the i th constraint ℓ ← ℓ + 1 S ℓ ← { k ∈ S \ T | ∃ j ∈ [ r i ] : q ijk > ǫ r i w ijT } , where w ijT is given by (10) T ← T ∪ S ℓ until | T | ≥ ǫ or S ℓ = ∅ or S \ T = ∅ if q iS ∈ P iT ℓ ( ǫ ) then T i ← T else ⊲ Find a subset of S that belongs to P iT ℓ ( ǫ ) if | T | ≥ ǫ then T ← the set of the first ǫ elements added to T h ← ; V ← S \ T else T i ← T Find a partition V , . . . , V h of S \ T s.t. ∃ j ∈ [ r i ] : P k ∈ V p q ijk ≥ ǫ r i w ijT ∀ p ∈ [ h ] K ← { k ∈ T | f ( S \{ k } ) ≥ (1 − ǫ ) f ( S ) } P ← { p ∈ [ h ] | f ( S \ V p ) ≥ (1 − ǫ ) f ( S ) } if K = ∅ then Pick arbitrary ˆ k ∈ K S ← S \{ ˆ k } else Pick arbitrary ˆ p ∈ P S ← S \ V ˆ p return ( { T , . . . , T m } , S ) emma 3.9. Consider sets T ⊆ S ⊆ [ n ] such that (C1) P k ∈ S q ik ∈ B ( r i , C i ) \ P iT ( ǫ ) ; (C2) q ijk ≤ ǫ r i w ijT ,for all j ∈ [ r i ] and k ∈ S \ T . Then there exist j ∈ [ r i ] , h ∈ [ ǫ − , r i ǫ ) , and a partition { V , . . . , V h } of S \ T such that P k ∈ V s q ijk ≥ ǫ r i w ijT for all s ∈ [ h ] .Proof. We define η , P k ∈ T q ik and κ , P k ∈ S \ T q ik . Since η + κ = P k ∈ S q ik ∈ B ( r i , C i ) \ P iT ( ǫ ) , we claimthat there is a j ∈ [ r i ] such that κ j > w ijT r i . Indeed, consider the ( r i − -dimensional simplex whose j th vertex,for j ∈ [ r i ] , is ( η , . . . , η j − , η j + w ijT , η j +1 , . . . , η r i ) . Then the point ξ ∈ R r i , defined by ξ j = η j + w ijT r i for j ∈ [ r i ] , lies on the simplex, implying that the whole box C := { ν ∈ R r i : η ≤ ν ≤ ξ } lies on the same side,including η , of the hyperplane H defined by the simplex. On the other hand, any facet of P iT ( ǫ ) is defined bythe convex hull of a set of points on the ball, which lie on the other side of H . Consequently, the point η + κ strictly lies on this other side, and hence outside the box C , implying the claim.Let us fix j ∈ [ r i ] as in the claim, and pack consecutive vectors q ik , k ∈ S \ T , into batches such that thesum in each batch has the j th component of length in the interval ( ǫ r i w ijT , ǫr i w ijT ] . More precisely, we fix anorder on S \ T := { , . . . , ℓ } , and find indices k < k < · · · < k h ′ < k h ′ +1 = ℓ + 1 such that k l +1 − X k = k l q ijk ≤ ǫr i w ijT , for l = 1 , . . . , h ′ (12)and k l +1 X k = k l q ijk > ǫr i w ijT for l = 1 , . . . , h ′ − (13)the existence of such indices is guaranteed by (C2). It follows from (13) that P k l +1 − k = k l q ijk > ǫ r i w ijT for l = 1 , . . . , h ′ − , since q ijk l +1 ≤ ǫ r i w ijT . It also follows that ǫ ≤ h ′ < r i ǫ + 1 , since summing (12) for ℓ = 1 , . . . , h ′ yields ǫr i h ′ w ijT ≥ h ′ X l =1 k l +1 − X k = k l q ijk = ℓ X k =1 q ijk = κ j ≥ w ijT r i ; (14)the last inequality follows from our assumption. Similarly, summing (13) for ℓ = 1 , . . . , h ′ − yields ( h ′ − ǫr i w ijT < h ′ − X l =1 k l +1 X k = k l q ijk ≤ r i X k =1 q ijk = 2 κ j ≤ w ijT , (15)where the last inequality follows from (C1). Setting V l , { k l , k l + 1 , . . . , k l +1 − } , for l = 1 , . . . , h ′ − , V h ′ − = { k h ′ − , . . . , r } , and h , h ′ − will satisfy the claim of the Lemma. Lemma 3.10.
Let f : 2 [ n ] → R + be a nonnegative submodular function and let S ⊆ [ n ] be a non-empty set.If { S , . . . , S k } is a partition of S into k disjoint subsets, then f ( S \ S i ) ≥ (1 − k ) f ( S ) for some i ∈ [ k ] .Proof. We prove the claim by contradiction. Suppose that f ( S \ S i ) < (1 − k ) f ( S ) holds for all i ∈ [ k ] . Weclaim by induction on i = 1 , , . . . , k that f ( S \ i [ j =1 S j ) < (cid:18) − ik (cid:19) f ( S ) , (16)12hich, when applied with i = k , would give the contradiction f ( ∅ ) < . The claim is true for i = 1 byassumption. Let us assume it is true up to i − . By the submodularity of f , we have: f (( S \ i − [ j =1 S j ) ∪ ( S \ S i )) + f (( S \ i − [ j =1 S j ) ∩ ( S \ S i )) ≤ f ( S \ i − [ j =1 S j ) + f ( S \ S i ) implying by the induction hypothesis and the assumption that f ( S \ S i ) < (1 − k ) f ( S ) that f ( S ) + f ( S \ i [ j =1 S j ) < (cid:18) − i − k (cid:19) f ( S ) + (1 − k ) f ( S ) , proving the claim.The following results are straightforward consequences of Theorem 3.7 and the results in [KST09]. Corollary 3.11.
There is a (1 − ǫ ) m ( − ǫ ) -approximation algorithm for the P ACK -S UB problem, for any ǫ > , when m and each r i , cp - rank( Q i ) are fixed. For the case when the submodular function is monotone,the approximation factor is (1 − ǫ ) m (1 − e − ǫ ) , for any ǫ > .Remark . Note that a slight modification of Algorithm 2 and Algorithm 3 yields a PTAS for the P
ACK -L IN problem. To show that, we do some following changes. Firstly, in the step 5 of Algorithm 2, the solution S willbe computed by a PTAS for the m -dimensional knapsack problem instead of using the α ( ǫ ) -approximationalgorithm for the submodular objective case. As a result, we have u ( ˆ S ) ≥ (1 − ǫ ) u ( S ′ ) ≥ (1 − ǫ ) m +1 O PT .Secondly, for Algorithm 3, we do the same steps 1-25. Then, the solution S i will be obtained from S i − bydropping the smallest utility-demand or group of demands with large component, ensuring that the set S i ofremaining demands can be packed inside P iT i ( ǫ ) . Hence, in the former case, u ( S i ) ≥ (1 − ǫ ) u ( S i − ) , and inthe latter case, u ( S i ) ≥ (1 − h ) S i − ≥ − ǫ − ǫ S i − ≥ (1 − ǫ ) S i − . − a Greedy-based Approach In this section we consider problem C
OVER - LIN with one quadratic constraint. We want to minimize a linearfunction f ( x ) = u T x , subject to x T Qx ≥ C and x ∈ { , } n , where Q ∈ CP ∗ n and u ∈ R n + . Note that theconvex programming-based method can not be applied here since the relaxed problem, which is obtained byconsidering x i ∈ [0 , for all i ∈ [ n ] , is non-convex , and thus we do not know if it can be solved efficiently.Furthermore, one can easily show that this programming relaxation has a bad integrality gap (of at least C ),and thus is not a good choice for approximation. Instead, we will follow a geometric approach.Let I = ( u, Q, C ) be an instance of the problem C OVER -L IN , and ǫ > be a fixed constant. We willconstruct a quasi-polynomial-time algorithm which produces an (1 + ǫ ) -approximate solution to the instance I . By guessing, we may assume that B ≤ O PT < (1 + ǫ ) B , where B := (1 + ǫ ) i min k u k for some i ∈ Z + (the number of possible guesses is O (log ǫ ( n · max k u k / min k u k )) ). Let T := { k | u k ≥ (1 + ǫ ) B } and V := { k | u k < ǫn · B } ; we set x k = 1 for all k ∈ T , and x k = 0 for all k ∈ V , and assume therefore that weneed to optimize over a set N := [ n ] \ ( T ∪ V ) , for which u k ∈ [ ǫn · B, (1 + ǫ ) B ) for k ∈ N . Note that suchrestriction increases the cost of the solution obtained by at most ǫ · O PT .As before, write Q = U U T , for some U ∈ Q n × r + , where r = cp - rank( Q ) , and for k ∈ [ n ] define thevector q k ∈ Q r + to be the k th column of U T . Define the conic region R T as in (9) (with the index i dropped).Then the problem now amounts to finding S ⊆ N s.t. q S := P k ∈ S q k is not in the interior of R T .13e begin by partitioning the set of vectors (indices) in N into h := r (cid:16) √ rǫ (cid:17) r − space-classes N , . . . , N h ,with the following property: for all s ∈ [ h ] , there exists ξ ( s ) ∈ R r + such that for all k ∈ N s , it holds q Tk ξ ( s ) k q k k k ξ ( s ) k ≥ − ǫ. (17)Condition (17) says that there is a fixed direction ξ ( s ) such that the angle that any vector q k , k ∈ N s , makeswith ξ ( s ) is sufficiently small.We will rely on the following geometric facts in our analysis of the algorithm. Fact 3.13.
Let a, b, ξ ∈ R r + be such that a T ξ k a k k ξ k ≥ − ǫ and b T ξ k b k k ξ k ≥ − ǫ . Then ( a + b ) T ξ k a + b k k ξ k ≥ − ǫ .Proof. Indeed, ( a + b ) T ξ k a + b k k ξ k = k a k k a + b k · a T ξ k a k k ξ k + k b k k a + b k · b T ξ k b k k ξ ( s ) k ≥ k a k + k b k k a + b k (1 − ǫ ) ≥ − ǫ, by the triangular inequality.Note that Fact 3.13 implies that, for any S ⊆ N s , q S := P k ∈ S q k also satisfies the condition (17). Fact 3.14.
Let a, b, ξ ∈ R r + be such that a T ξ k a k k ξ k ≥ − ǫ and b T ξ k b k k ξ k ≥ − ǫ . Then a T b k a k k b k ≥ − ǫ .Proof. Let ¯ a := a k a k , ¯ b := b k b k , and ¯ ξ := ξ k ξ k . If a, b, ξ all lie in the same subspace, then the claim followssince the angle between a and b is no more than the sum of the angles between a and ξ , and b and ξ , whichis at most is cos − ((1 − ǫ ) − ǫ (2 − ǫ )) ≤ cos − (1 − ǫ ) . Otherwise, let ¯ b = b b + ˜ b be the orthogonaldecomposition of ¯ b with respect to the 2-dimensional subspace formed by the two vectors ¯ a and ¯ ξ , where b b is the projection of ¯ b into this space, and ˜ b is the orthogonal component. Then b b T ¯ ξ k b b k = ¯ b T ¯ ξ k b b k ≥ − ǫ k b b k ≥ − ǫ (since k b b k ≤ k ¯ b k = 1 ), which also implies that k b b k ≥ − ǫ . Since a , b b , and ξ lie in the same subspace, itfollows by the above argument that b b T ¯ a k b b k ≥ (1 − ǫ ) , and hence, ¯ b T ¯ a = b b T ¯ a ≥ (1 − ǫ ) k b b k ≥ (1 − ǫ )(1 − ǫ ) ,implying the claim. Fact 3.15.
Let a, b ∈ R be such that a T b k a k k b k ≥ − ǫ and k b b k = λ k a k , where λ ≥ and b b := Pj a ( b ) isthe projection of b on a . Then for any vector η ∈ R r , it holds that k Pj η ( b ) k ≥ λ k Pj η ( a ) k − p ǫ (2 − ǫ )1 − ǫ ! k a k . Proof.
Since the statement is invariant under rotation, we may assume w.l.o.g. that η = j , the j th-dimensional unit vector in R r . Write b := b b + ˜ b , where ˜ b is the vector orthogonal to a in the subspacespanned by a and b . Then k Pj η ( b ) k = b j is the j th component of b , and k Pj η ( b b ) k = b b j . Since k b b k ≥ a T b b k a k = a T b k a k ≥ (1 − ǫ ) k b k , it follows that k ˜ b k = q k b k − k b b k ≤ p ǫ (2 − ǫ ) k b k ≤ p ǫ (2 − ǫ )1 − ǫ k b b k = λ p ǫ (2 − ǫ )1 − ǫ k a k . Since | b j − b b j | = | ˜ b j | ≤ k ˜ b k , and b b = λ k a k · a , the claim follows.14ondition (17) together with Facts 3.13 and 3.14 imply that for any two sets S, S ′ ⊆ N s , we have q TS q S ′ k q S k k q S ′ k ≥ − ǫ .The space partitioning can be done as follows. Let q T := P k ∈ T q k . We partition the region R T intodisjoint regions R T (1) , . . . , R T ( h ) , obtained as follows. Let µ = ¯ w T · , where ¯ w T := max j ∈ [ r ] w jT (recallthe definition of w jT from (10)), and define the r -dimensional box C T := { ν ∈ R r + | q T ≤ ν ≤ µ } . Notethat R T ⊆ C T . We grid the r facets of C T that do not contain the point q T by interlacing equidistant ( r − -dimensional parallel hyperplanes with inter-separation ǫ ¯ w T √ r , for each j ∈ [ r ] , and with the j th principalaxis as their normal. Note that the total number of grid cells obtained is h ; let us call them C , . . . , C h (these are ( r − -dimensional hypercubes). We then define the region R ( s ) as the r -dimensional simplex R ( s ) := conv( { q T } ∪ C s ) and ξ ( s ) = c ( s ) − q T as the designated vector, where c ( s ) is the vertex centerof the cell C s . Note that the angle condition (17) is satisfied. Indeed, consider any vector q k such that q T + q k ∈ R T ( s ) . Let the ray { q T + λq k : λ ≥ } hit the boundary cell C s in the point x . Consider thetriangle formed by the three points q T , c ( s ) and x . Then, by construction, the distances between q T and both c ( s ) and x are at most ¯ w , whereas the distance between c ( s ) and x is at most ǫ ¯ w √ r . It follows that the anglebetween the two vectors q k and ξ ( s ) is no more than sin − ǫ , implying (17).Finally, we define N s := { k | q T + q k ∈ R T ( s ) } . This gives the required space partitioning of the vectors.Next, we group the set of items in N into ℓ := 1 + log ǫ nǫ (some possibly empty) utility-classes N , . . . , N ℓ , where N l = { k | u k ∈ [ ǫn · B (1 + ǫ ) l − , ǫn · B (1 + ǫ ) l ) } . Note that for all k, k ′ ∈ N l , we have u k ≤ u k ′ (1 + ǫ ) . (18)We can show that for a set of vectors { q k : k ∈ N s ∩ N l } that lie in the same region and same utility-group, the greedy algorithm that processes the vectors in non-increasing order of length gives an O ( ǫ ) -optimalsolution. For simplicity we assume first that T = ∅ and q T = . Algorithm 4 G REEDY -C OVER ( u, { q k } k ∈ N , C ) Input:
A cost vector u ∈ R N + ; accuracy parameter ǫ ; vectors q k ∈ Q r + satisfying (18) and (17); a demand C ∈ R + ; Output: ǫ -optimal solution S to C OVER -L IN S ′ := min { u ( S ) | S ⊆ N, | S | ≤ ǫ , q S is feasible } Order the vectors q k , k ∈ N , such that k q k ≥ k q k ≥ · · · S ← ∅ ; k ← while k P k ∈ S q k k < C do k ← k + 1 S ← S ∪ { k } if u ( S ) ≤ u ( S ′ ) then return S else return S ′ Lemma 3.16.
Consider instance of problem C OVER -L IN described by a set of vectors { q k } k ∈ N satisfying(17) and (18). Then, for any sufficiently small constant ǫ > , Algorithm 4 outputs a solution S satisfying u ( S ) ≤ (1 + 9 ǫ ) O PT .Proof. Let S ∗ be an optimal solution. Since we consider every possible feasible solution of size at most ǫ instep 1, we may assume w.l.o.g. that | S ∗ | ≥ ǫ .We claim that | S | ≤ | S ∗ | − ǫ + 1 . To see this claim, let q S ∗ := P k ∈ S ∗ q k , and for k ∈ N , denote by b q k := Pj q S ∗ ( q k ) the projection of q k on q S ∗ . Let q w ∈ S be the last item added to S in step 6 of the algorithm.15hen it is clear that X k ∈ S \{ w } k b q k k < k q S ∗ k , (19)since otherwise S ′ := S \{ w } would satisfy the condition in line 4 (as k P k ∈ S ′ q k k ≥ P k ∈ S ′ k b q k k ≥k q S ∗ k ≥ C , by the feasibility of S ∗ ). By the angle condition (17): k b q k k ≥ (1 − ǫ ) k q k k for all k ≤ N. Itfollows by (19) that (1 − ǫ ) | S ′ | min k ∈ S ′ k q k k ≤ (1 − ǫ ) X k ∈ S ′ k q k k ≤ X k ∈ S ′ k b q k k < k q S ∗ k = X k ∈ S ∗ k b q k k ≤ X k ∈ S ∗ k q k k ≤ | S ∗ | max k ∈ S ∗ k q k k . The claim follows since max k ∈ S ∗ k q k k ≤ min k ∈ S ′ k q k k , by the greedy order in line 2.By the utility condition (18), u ( S ) ≤ | S | max k u k ≤ | S | (1 + ǫ ) min k u k ≤ (cid:18) | S ∗ | − ǫ + 1 (cid:19) (1 + ǫ ) min k u k ≤ (cid:18) − ǫ + ǫ (cid:19) (1 + ǫ ) u ( S ∗ ) . The lemma follows.Our QPTAS is given as Algorithm 5. For simplicity, we assume that the algorithm has already a correctguess of the bound B on the value of the optimal solution. After decomposing the instance into classesaccording to utility and region (steps 2 and 4), the algorithm enumerates over all possible selections of anonnegative integer n s,l associated to each region s and a utility class l (step 8); this number n s,l representsthe largest length vectors that are taken in the potential solution from the set N s ∩ N l . However, for technicalreasons, the algorithm does this only for pairs ( s, l ) for which the set N s ∩ N l contributes at least ǫ in theoptimal solution; the set of pairs that potentially do not satisfy this can be identified by enumeration (steps 6and 7). Algorithm 5 C OVER -L IN -QPTAS ( u, Q, C ) Input:
A cost vector u ∈ R N + ; accuracy parameter ǫ ; matrix Q ∈ CP ∗ n ; a demand C ∈ R + Output: An O ( √ ǫr ) -optimal solution S to C OVER -L IN ⊲ r = cp-rank ( Q ) Obtain sets T and V as explained above and set N ← [ n ] \ ( T ∪ V ) Decompose Q as Q = U U T , where U ∈ Q n × r + Let { q k } k ∈ N be the columns of U T corresponding to the indices in N Decompose the set N into space-classes N , . . . , N h and utility classes N , . . . , N ℓ S ← T ∪ N ⊲
Assume instance is feasible for each subset of pairs G ⊆ [ h ] × [ ℓ ] do for each possible selection ( T s,l ⊆ N s ∩ N l : | T s,l | ≤ ǫ , ( s, l ) ∈ G ) do for each possible selection ( n s,l ∈ { , . . . , | N s ∩ N l |} : s ∈ ([ h ] × [ ℓ ]) \ G ) do Let S s,l be the set of the n s,l vectors with largest length in N s ∩ N l S ′ ← T ∪ (cid:16)S ( s,l ) ∈G T s,l (cid:17) ∪ (cid:16)S s,l S s,l (cid:17) if k P k ∈ S ′ q k k ≥ C and u ( S ′ ) < u ( S ) then S ← S ′ return S Lemma 3.17.
For any sufficiently small ǫ > (for instance, ǫ < r +1) for r ≥ ), Algorithm 5 runs intime n O (( √ rǫ ) r +1 log n ) and outputs a solution S satisfying u ( S ) ≤ (1 + O ( √ ǫr )) O PT for any instance of C OVER -L IN . roof. The running time is obvious since it is dominated by the number of selections in steps 6, 7 and 8,which is at most (cid:0) nhℓ (cid:1) O ( hℓǫ ) = n O (( √ rǫ ) r +1 log n ) .To see the approximation ratio, fix λ := − p ǫ (2 − ǫ ) r − ǫ ! − . (20)Let S ∗ be an optimal solution (within N ), and for s ∈ [ h ] and l ∈ [ ℓ ] , let S ∗ s,l := S ∗ ∩ N s ∩ N l . Let G := { ( s, l ) ∈ [ h ] × [ ℓ ] : | S ∗ s,l | < ǫ } and T s,l := S ∗ ∩ N s ∩ N l , for ( s, l ) ∈ G .Let b S := S ∗ \ (cid:16)S ( s,l ) ∈G ( N s ∩ N l ) (cid:17) . For k ∈ N s ∩ N l , let b q k := Pj q S ∗ s,l ( q k ) be the projection of q k on q b S s,l := P k ∈ b S s,l q k . Define the set of pairs H := ( s, l ) ∈ [ h ] × [ ℓ ] : k X k ∈ N s ∩ N l b q k k ≥ λ k q b S s,l k and | S ∗ s,l | ≥ ǫ . Then according to Lemma 3.16 (or more precisely, its proof), for every ( s, l ) ∈ H there a choice n s,l ∈ { , . . . , | N s ∩ N l |} , such that if S s,l is the set of the n s,l vectors with largest length in N s ∩ N l ,then P k ∈ S s,l k b q k k ≥ λ k q b S s,l k and u ( S s,l ) ≤ (cid:16) λ − ǫ + ǫ (cid:17) (1 + ǫ ) u ( b S s,l ) . By this, we can define p k = τ · q k ,for τ := λ k q b Ss,l k P k ∈ Ss,l k b q k k ≤ , such that P k ∈ S s,l k b p k k = λ k q b S s,l k , where b p k := Pj q b Ss,l ( p k ) .Let us now apply Fact 3.15 with a = a s,l := q b S s,l , b = b s,l := p S s,l , and η = q b S to get that k Pj η ( b s,l ) k ≥ λ k Pj η ( a s,l ) k − p ǫ (2 − ǫ )1 − ǫ ! k a s,l k . (21)Summing the above inequalities for all ( s, l ) ∈ H , we get X ( s,l ) ∈H k Pj η ( b s,l ) k ≥ λ X ( s,l ) ∈H k Pj η ( a s,l ) k − p ǫ (2 − ǫ )1 − ǫ X ( s,l ) ∈H k a s,l k . (22)Note that P ( s,l ) ∈H a s,l = η by construction. By the Cauchy-Schwarz inequality and the nonnegativity of thevectors, k X ( s,l ) ∈H Pj η ( a s,l ) k = k η k = k X ( s,l ) ∈H a s,l k ≥ √ r k X ( s,l ) ∈H a s,l k = 1 √ r X ( s,l ) ∈H k a s,l k ≥ √ r X ( s,l ) ∈H k a s,l k . Using this in (22), we obtain X ( s,l ) ∈H k Pj η ( b s,l ) k ≥ λ − p ǫ (2 − ǫ ) r − ǫ ! X ( s,l ) ∈H k Pj η ( a s,l ) k ≥ X ( s,l ) ∈H k Pj η ( a s,l ) k , (23)by our choice of λ . From (22) and q S s,l ≥ b s,l , it follows by the feasibility of S ∗ that the solution defined by S = T ∪ (cid:16)S ( s,l ) ∈G T s,l (cid:17) ∪ (cid:16)S s,l S s,l (cid:17) is feasible. 17learly, one of the choices in each of the enumeration steps 6, 7, and 8 will capture the above choices G , T s,l for ( s, l ) ∈ G , and n s,l , for ( s, l ) ∈ ([ h ] × [ ℓ ]) \G . It follows the procedure returns a solution S with utility: u ( S ) ≤ u ( T ) + X ( s,l ) ∈G u ( T s,l ) + X ( s,l ) u ( S s,l ) ≤ u ( T ) + X ( s,l ) ∈G u ( T s,l ) + (cid:18) λ − ǫ + ǫ (cid:19) (1 + ǫ ) X ( s,l ) u ( b S s,l ) ≤ − ǫ − p ǫ (2 − ǫ ) r + ǫ ! (1 + ǫ ) u ( T ) + X ( s,l ) ∈G u ( T s,l ) + X ( s,l ) u ( b S s,l ) ≤ (1 + O ( √ ǫr )) u ( S ∗ ) . The lemma follows.
According to Corollary 2.3, given a matrix Q i of fixed cp-rank r i , a decomposition of the form Q i := U i U Ti − ∆ i , where U i ∈ Q r i × n , ∆ i ∈ R n × n + , and k ∆ i k ∞ ≤ δ , can be done in time poly( L , n O ( r i ) , log δ ) for any δ > . This leads to a technical issue: if one uses the approximate decomposition Q i := U i U Ti to solve theBQC problem, then the resulting solution can be either infeasible or far from optimal. When m = O (1) ,this issue can be dealt with at a small loss in the approximation ratio as follows. Let Q kji denote the ( k, j ) thentry of Q i . Note that, if the k th diagonal element Q kki of Q i is , then every entry in the k th row and the k thcolumn of Q i is , since Q i is positive semi-definite.We consider the case when f is linear. Let S ∗ be an optimal solution for P ACKING -BQC (resp., C
OVERING -BQC) and x ∗ be its characteristic vector. The idea is to show that there is a feasible solution e S ⊆ S ∗ (resp., e S ⊇ S ∗ ) such that f ( e S ) ≥ (1 − ǫ ) f ( S ∗ ) (resp., f ( e S ) ≤ (1 + ǫ ) f ( S ∗ ) ), and such that e x T Q i e x ≤ e C i (resp., e x T Q i e x ≥ e C i ), where e x is the characteristic vector of e S and e C i := C i − δn (resp., e C i := C i + δn ).Then it would be enough to solve a new instance with the same objective function but with the constraints x T U i U Ti x ≤ C i (resp., x T U i U Ti x ≥ e C i ). Note that e x is feasible for the new instance, since e x T U i U Ti e x = e x T Q i e x + e x T ∆ i e x ≤ e C i + n δ = C i (resp., e x T U i U Ti e x = e x T Q i e x + e x T ∆ i e x ≥ e C i ).Suppose that x is an α -approximate solution for the new instance of P ACKING -BQC, then x T Q i x = x T U i U Ti x − x T ∆ i x ≤ C i , since ∆ i ≥ , and hence x is an α (1 − ǫ ) -approximation for the original in-stance. Similarly, if x is an α -approximate solution for the new instance of C OVERING -BQC, then x T Q i x = x T U i U Ti x − x T ∆ i x ≥ e C i − δn = C i , and hence x is an α (1 + ǫ ) -approximation for the original instance.The ǫ -approximate solution e S can be obtained as follows. Let ǫ > be a fixed constant, and define λ = mǫ . Let U be the λ items of highest utility in the optimal solution S ∗ . This gives a partition of [ n ] intothree sets U ∪ V ∪ N , as described in Section 3.1. We can make the following assumptions w.l.o.g.:(A1) | S ∗ | ≥ mǫ ;(A2) for each i ∈ [ m ] the set κ i := { k ∈ S ∗ ∩ N : Q kki > } 6 = ∅ (resp., κ i := { k ∈ N \ S ∗ : Q kki > } 6 = ∅ ).(A1) can be assumed since otherwise the algorithm that enumerates over all possible sets of size at most mǫ isoptimal. (A2) can be assumed since if it is violated by i ∈ [ m ] , then after fixing the set U (resp, after fixingthe set U and fixing the variables in the set { k ∈ N : Q kki = 0 } to ), the i th constraint becomes redundant.Thus the algorithm that enumerates all subsets U of size at most mǫ and all subsets of constraints, removing asubset at a time (resp., enumerates all subsets U of size at most mǫ and all subsets of constraints, removing asubset at a time and fixing the variables in the set { k ∈ N : Q kki = 0 } to ), can assume an instance satisfying(A2). 18hen, assuming (A1) and (A2), for each i ∈ [ m ] , we pick an index k ( i ) ∈ κ i and set the correspondingvariable to (resp., to ). This yields an ǫ -optimal solution e S , whose characteristic vector x also satisfies x T Q i x ≤ ( x ∗ ) T Q i x ∗ − Q k ( i ) ,k ( i ) i ≤ e C i (resp., x T Q i x ≥ ( x ∗ ) T Q i x ∗ + Q k ( i ) ,k ( i ) i ≥ e C i ), if we choose δ := min i,k Q kki n . In this section we study optimization problem with multiple objectives. In such a problem it is concernedwith optimizing more than one objective function simultaneously. In fact, we consider the following binarymulti-objective quadratic problem (BMOQ):(BMOQ) { max f ( x ) , min g ( x ) } subject to x ∈ X ⊆ { , } n where f ( x ) = { f i ( x ) } i ∈I and g ( x ) = { g j ( x ) } j ∈J , and X = { x ∈ { , } n | x T Q i x + q Ti x ≤ C i ; Ax ≤ b ; i ∈ [ m ] } , where Q i ∈ CP ∗ n , q i ∈ R n + for all i ∈ [ m ] , A ∈ R d × n + , b ∈ R d + , and m, d are constant.Typically, a feasible solution that simultaneously optimizes each objective may not exist due to the trade-off between the different objectives. This issue can be captured by the notion of Pareto-optimal frontier (or
Pareto curve ), which is a set of all feasible solutions whose vector of the objectives is not dominated by anyother one. Unfortunately, the size of such a Pareto-optimal frontier set is often exponential in size of theinput. Hence, a natural goal in this setting is, given an instance of a (BMOQ) problem, to efficiently achievean ǫ -approximation of the Pareto optimality that consists of a number of solutions that is polynomial in sizeof the instance, for every constant ǫ > . Formally, let I be an instance of the multi-objective optimizationproblem (BMOQ), consider the following definitions taken from [MS13b]: Definition 4.1.
A Pareto-optimal frontier P of I is a subset of X , such that for any solution S ∈ P there isno solution S ′ ∈ X such that f i ( S ′ ) ≥ f i ( S ) for all i ∈ I and g j ( S ′ ) ≤ g j ( S ) for all j ∈ J , with strictinequality for at least one of them. Definition 4.2.
For ǫ > , an ǫ -approximate Pareto-optimal frontier of I , denoted by P ǫ , is a set of solutions,such that for all S ∈ X , there exists a solution S ′ ∈ P ǫ such that f j ( S ′ ) ≥ (1 − ǫ ) f j ( S ) and g j ( S ′ ) ≤ g j ( S ) / (1 − ǫ ) for all i ∈ I , j ∈ J . Computing an (approximate) Pareto-optimal frontier has been known to be a very efficient approach indesigning approximation schemes for combinatorial optimization problems with multiple objectives [EKP02,PY00, MS13b, MS13a]. Erlebach et al. [EKP02] give a PTAS for multi-objective multi-dimensional knapsackproblem. In this paper, we will extend their method to the case with quadratic constraints. Let I be an instanceof the BMOQ problem, and denote f j ( x ) = a Tj x, g j ( x ) = b Tj x for all j ∈ [ p ] . Let ǫ > be a fixed constant.The Algorithm P ARETO -O PT below will produce a set P ǫ , which is an ǫ -approximate Pareto-optimal frontierto the instance I , in polynomial time. Moreover, the size of P ǫ is guaranteed to be also polynomial in size ofthe input for every fixed ǫ > .We first compute ǫ -approximate solutions S j , T j to the maximization problem with single objective func-tions f j = a Tj x, g j = b Tj x subject to x ∈ X , respectively, by using an extension of Algorithm 1. As a result,for any feasible solution S , we have f j ( S ) ≤ (1 − ǫ ) − f j ( S j ) and g j ( S ) ≤ (1 − ǫ ) − g j ( T j ) for all j ∈ [ p ] .For each j ∈ [ p ] , we consider v j + 2 (lower) bounds { α jλ j , ≤ λ j ≤ v j + 1 } for the objective f j and w j + 3 (upper) bounds { β jη j , ≤ η j ≤ w j + 2 } for the objective g j . We define the set Λ × Γ ⊆ R p − , which19ontains all the tuples ( λ , . . . , λ p − , η , . . . , η p ) . Note that the size of Λ × Γ is bounded by a polynomial insize of the input since p is constant and v j , w j are also bounded by a polynomial in size of the input. The mainidea of the algorithm is that, for each tuple of the bounds, we try to maximize the last objective f p subjectto f j ≥ α jλ j , for all j ∈ [ p − , and g j ≤ β jη j , for all j ∈ [ p ] . To do that, we consider all possibility ofchoosing a subset of [ n ] of at most (a constant number) λ of items. Denote Σ as the set of all such subsets andlet X = Σ p . For each tuple ( X , . . . , X p ) ∈ X , where X j is a set of items of highest utility with respect to f i , we construct a tuple ( Y , . . . , Y p ) such that Y j is the set of all items that will not be put into the knapsackgiven that items in X j are the ones with highest utility with respect to f j in the solution. Now for each tuple ( X , . . . , X p , Y , . . . , Y p ) such that S j ∈ [ p ] ( X j ∩ Y j ) = ∅ , we solve the convex program (CP2 [ U ] ) and obtain(if exists) an ǫ -optimal solution x ∗ :(CP2 [ U ] ) max a Tp x subject to x T Q i x + q Ti x ≤ C i , i ∈ [ m ] ,Ax ≤ b,a Tj x ≥ α jλ j , j ∈ [ p − ,b Tj x ≤ β jη j , j ∈ [ p ] ,x k = 1 , for k ∈ U , x k = 0 , for k ∈ V ,x k ∈ [0 , , for k ∈ N. Define t i := U Ti [ ∗ ; N ] x ∗ N , t ′ i := T U Q i [ U ; N ] x ∗ N , q Ti [ ∗ ; N ] x ∗ N = v i for all i ∈ [ m ] , A [ ∗ ; N ] x ∗ N = σ , a Tj [ ∗ ; N ] x ∗ N = θ j for j ∈ [ p − , and b Tj [ ∗ ; N ] x ∗ N = θ ′ j for j ∈ [ p ] . We define a polytope P ( U ) ⊆ [0 , N asfollows: P ( U ) = y ∈ [0 , N (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) U Ti [ ∗ ; N ] y ≤ t i , T U Q i [ U ; N ] y ≤ t ′ i , q Ti [ ∗ ; N ] y ≤ v i , for i ∈ [ m ] A [ ∗ ; N ] y ≤ σa Tj [ ∗ ; N ] y ≥ θ j for j ∈ [ p − , b Tj [ ∗ ; N ] y ≤ θ ′ j for j ∈ [ p ] We can find (if exists) a (BFS) y in this polytope such that a Tp y ≥ a Tp x ∗ N . By rounding down this fractionalsolution y and setting x k ∈ { , } according to the assumption k ∈ U ∪ V , we obtain an integral solution x ∈ { , } n .There are totally n pλ tuples of the form ( X , . . . , X p , Y , . . . , Y p ) . For each such tuple such that thecondition in Step 18 is satisfied and (CP1 [ U ] ) is feasible, the algorithm outputs an integral solution x . Allpossible integral solutions obtained in this way are collected in the set P ǫ . This completes the description ofthe algorithm. Theorem 4.3.
Let I be an instance of the BMOQ problem. For every fixed ǫ > , the Algorithm 6 runs inpolynomial time in size of the input and produces an ǫ -approximate Pareto-optimal frontier P ǫ for I .Proof. It is not difficult to see that the running time of Algorithm 6 is polynomial in size of the input as theparameters ǫ, ǫ ′ , µ, λ are constants. Hence, we need only to show the approximation factor of the algorithm.More precisely, let S be an arbitrary feasible solution, we need to prove that there exists a solution S ′ ∈ P ǫ such that f j ( S ′ ) ≥ (1 − ǫ ) f j ( S ) and g j ( S ′ ) ≤ g j ( S ) / (1 − ǫ ) , for all j ∈ [ p ] . Let us define λ j = max { λ | α jλ ≤ f j ( S ) } and η j = min { η | β jη ≥ g j ( S ) } . It follows that α jλ j ≤ f j ( S ) ≤ α jλ j / (1 − ǫ ) and β jη j ≥ g j ( S ) ≥ (1 − ǫ ) β jη j Let δ = min { λ, | S |} . For each j ∈ [ p ] , let X j be the set that contains δ items in S with largest util-ity with respect to the objective f j . One can see that S is a feasible solution to (CP2 [ U ] ) with the tuple20 lgorithm 6 P ARETO -O PT ( { a Tj , b Tj , Q i , C i } j ∈ [ p ] , i ∈ [ m ] , ǫ ) Input:
Utilities, demand vectors and capacities ( { a Tj , b Tj , Q i , C i } j ∈ [ p ] , i ∈ [ m ] ) ; accuracy parameter ǫ Output: (1 − ǫ ) -approximate Pareto optimality P ǫ P ǫ ← ∅ Decompose Q i into U Ti U i , where U i has r i rows, for all i ∈ [ m ] ⊲ r i = cp-rank ( Q i ) for j ∈ [ p ] do S j ← PTAS ( a Tj , Q i , C i , i ∈ [ m ] , ǫ ) ⊲ Find an ǫ -approximate solution to f j T j ← PTAS ( b Tj , Q i , C i , i ∈ [ m ] , ǫ ) ⊲ Find an ǫ -approximate solution to g j v j ← l log (1 − ǫ ′ ) − f j ( S j ) m ; w j ← l log (1 − ǫ ′ ) − g j ( T j ) m α j ← ; α jλ j ← (1 − ǫ ) − λ j for all λ j ∈ [ v j + 1] β j ← ; β jη j ← (1 − ǫ ) − η j for all η j ∈ [ w j + 2] Λ ← [0 , v ] × . . . × [0 , v p − ] ; Γ ← [0 , w ] × . . . × [0 , w p ] φ ← P mi =1 ( r i + 1) ; λ ← (cid:24) ( φ + 2 p + d + m − ǫ ) ǫ (1 − ǫ ) (cid:25) Y ← ∅ ; X ← { ( X , . . . , X p ) | X j ⊆ [ n ] , | X j | ≤ λ, j ∈ [ p ] } for each ( X , . . . , X p ) ∈ X do for j ∈ [ p ] do Y j ← { k ∈ [ n ] \ X j | a jk > min { a jk ′ | k ′ ∈ X j }} Y ← Y ∪ ( Y , . . . , Y p ) for each tuple ( λ , . . . , λ p − , η , . . . , η p ) ∈ Λ × Γ do for each tuple ( X , . . . , X p , Y , . . . , Y p ) ∈ X × Y do if U ∩ V = ∅ then Find (if exists) an ǫ -optimal solution x ∗ to the convex program (CP2 [ U ] ) Define a polytope P ( U ) ⊆ [0 , N Find an BFS y of P ( U ) such that a Tp y ≥ a Tp x ∗ N x ← { ( x , . . . , x n ) | x k = ⌊ y k ⌋ for k ∈ N, x k = 1 , for k ∈ U , x k = 0 , for k ∈ V} ⊲ rounding down solution y P ε ← P ε ∪ x return P ε ( λ , . . . , λ p − , η , . . . , η p ) and the tuple ( X , . . . , X p ) . Let y be a BFS to the polytope P ( U ) and let x be thecorresponding rounded solution. By the same arguments provided as in the section 26, one can easily proofthat x is feasible. Denote S ′ by set of items corresponding to the integral solution x . We have: g j ( S ′ ) = b Tj x ≤ β jη j ≤ g j ( S ) / (1 − ǫ ) ≤ g j ( S ) / (1 − ǫ ) (24)Furthermore, for every j ∈ [ p − : a Tj x ∗ ≥ α jλ j ≥ (1 − ǫ ) f j ( S ) , and for j = p , we even have a Tp x ∗ ≥ (1 − ǫ ) opt ≥ (1 − ǫ ) f p ( S ) , because x ∗ is ǫ -optimal to (CP2 [ U ] ). Weconsider two cases below.Case I: If δ = | S | , (i.e. | S | ≤ λ ), the set S ′ will contain S as a subset. Hence, we have f j ( S ′ ) ≥ f j ( S ) for all j ∈ [ p ] .Case II: If δ = λ , we prove that f j ( S ′ ) ≥ (1 − ǫ ) f j ( S ) for all j ∈ [ p ] . Let a jk be the smallest utilityamong the λ items with highest utility in S with respect to objective f j , we have a Tj y + a Tj U ≥ a Tj x ∗ N + a Tj U = a Tj x ∗ ≥ (1 − ǫ ) O PT ≥ (1 − ǫ ) λa jk
21n the other hand, since x is obtained by rounding down y , and the fact that y has at most φ + 2 p + d + m − fractional components, we have f j ( S ′ ) = a Tj x = a Tj x N + a Tj U ≥ a Tj y − a jk ( φ + 2 p + d + m −
1) + a Tj U Therefore, f j ( S ′ ) ≥ a Tj x ∗ − φ + 2 p + d + m − λ λa jk ≥ a Tj x ∗ − φ + 2 p + d + m − λ (1 − ǫ ) a Tj x ∗ = (1 − φ + 2 p + d + m − λ (1 − ǫ ) ) a Tj x ∗ ≥ (1 − ǫ ǫ ) a Tj x ∗ = a Tj x ∗ ǫ Note that a Tj x ∗ ≥ (1 − ǫ ) f j ( S ) , we have f j ( S ′ ) ≥ a Tj x ∗ ǫ ≥ − ǫ ǫ f j ( S ) = (1 − ǫ ) f j ( S ) . (25)From (24) and (25), the proof is completed. We consider the binary quadratically constrained quadratic programming problem which takes the form:(BQCQP) max w ( x ) = x T Qx + q T x subject to x ∈ X ⊆ { , } n , where Q ∈ CP ∗ n , q ∈ R n + , and the feasible region X is defined as in the previous section. Theorem 4.4.
The problem (BQCQP) admits a PTAS.Proof.
Let I be an instance of (BQCQP) and let ǫ > be a fixed constant. Since p = cp-rank ( Q ) is fixed, onecan find in polynomial time a factorization Q = P ri =1 a j a Tj , a j ∈ R n + . Hence, the objective function can bewritten in the form P pj =1 ( a Tj x ) + a Tp +1 x , where a p +1 = q . Applying the Theorem 4.3 for the multi-objectiveinstance with |I| = p + 1 , |J | = 0 , we obtain P ǫ as an approximation of the Pareto-optimal frontier. Let x ∗ ∈ X be an optimal solution to the instance I , there exists a solution y ∈ P ǫ such that a Tj y ≥ (1 − ǫ ) a Tj x ∗ ,for all j ∈ [ p + 1] . Hence, w ( y ) = p X j =1 ( a Tj y ) + a Tp +1 y ≥ (1 − ǫ ) p X j =1 ( a Tj x ∗ ) + (1 − ǫ ) a Tp +1 x ∗ ≥ (1 − ǫ ) O PT . Let y ∗ = argmax { w ( x ) | x ∈ P ǫ } (note that the size of P ǫ is polynomial in size of the input), we have w ( y ∗ ) ≥ w ( y ) ≥ (1 − ǫ ) O PT . We consider the following problem where the objective is the product of a number of linear functions.(BMPP) max w ( x ) = ( a T x ) · ( a T x ) · . . . · ( a Tp x ) subject to x ∈ X ⊆ { , } n , where a j ∈ R n + for all j ∈ [ p ] . 22 heorem 4.5. If p is constant, the problem (BMPP) admits a PTAS.Proof. Again, applying the Theorem 4.3 for the multi-objective instance with |I| = p, |J | = 0 , we obtain anapproximation of the Pareto set P ǫ . Let x ∗ ∈ X be an optimal solution to the instance of (BMPP), there existsa solution y ∈ P ǫ such that a Tj y ≥ (1 − ǫ ) a Tj x ∗ , for all j ∈ [ p ] . Hence, w ( y ) = p Y j =1 a Tj y ≥ (1 − ǫ ) p p Y j =1 ( a Tj x ∗ ) = (1 − ǫ ) p O PT . Again, let y ∗ = argmax { w ( x ) | x ∈ P ǫ } , it follows w ( y ∗ ) ≥ (1 − ǫ ) p O PT . Note that p is constant, the resultfollows. We consider the binary quadratically constrained programming with a rational objective function.(S UM -R ATIO ) max w ( x ) = X pj =1 a Tj xb Tj x subject to x ∈ X ⊆ { , } n . This problem belongs to the class of fractional programming which has important applications in severalareas such as transportation, finance, engineering, statistics (see the survey paper by [SS03] and the referencestherein, for more applications).
Theorem 4.6.
Suppose that p is fixed. The problem (S UM -R ATIO ) admits a PTAS.Proof. applying the Theorem 4.3 for the instance with |I| = |J | = p , we obtain an approximation of thePareto set P ǫ of the instance of multi-objective problem. Let x ∗ ∈ X be an optimal solution to an instance of(S UM -R ATIO ) problem, there exists a solution y ∈ P ǫ such that a Tj y ≥ (1 − ǫ ) a Tj x ∗ and b Tj y ≤ b Tj x ∗ / (1 − ǫ ) ,for all j ∈ [ p ] . Hence, w ( y ) = p X j =1 a Tj yb Tj y ≥ (1 − ǫ ) p X j =1 a Tj x ∗ b Tj x ∗ = (1 − ǫ ) O PT . Again, let y ∗ = argmax { w ( x ) | x ∈ P ǫ } , the proof is completed. In this paper, we have investigated the approximability of binary quadratic programming problems when thecp-rank of the constraints’ matrices are bounded by a constant. It turns out that bounding the cp-rank of thesematrices makes several interesting variants of the quadratic programming problem tractable. In particular, ourresults hint that limiting the cp-rank makes the quadratic problem exhibit similar approximability behavior asthe linear case, assuming a constant number of packing or covering quadratic constrains. For the case with anynumber of quadratic constraints and linear objective function, one can easily obtain a O ( m √ r ) -approximationalgorithm for the packing problem, where r = max i ∈ [ m ] r i . The first interesting question is whether thereexists a (Q)PTAS for the covering problem with a constant number of quadratic constraints. Extending ourresults to the case of low nonnegative ranks, and finding non-trivial inapproximability bounds, parametrizedby the cp-rank of the constraints matrices are also another interesting open problems.23 eferences [AFLS01] K. Allemand, K. Fukuda, T. Liebling, and E. Steiner. A polynomial case of unconstrained zero-one quadratic optimization. Math. Program. , 91(1):49–52, 2001.[AGKM12] S. Arora, R. Ge, R. Kannan, and A. Moitra. Computing a nonnegative matrix factorization -provably. In
STOC , pages 145–162, 2012.[AKA94] B. Alidaee, G. Kochenberger, and A. Ahmadian. 0-1 quadratic programming approach for theoptimal solution of two scheduling problems.
International Journal of Systems Science , 25:401–408, 1994.[BPR96] S. Basu, R. Pollack, and M. Roy. On the combinatorial and algebraic complexity of quantifierelimination.
J. ACM , 43(6):1002–1045, 1996. Preliminary version in
FOCS
Completely Positive Matrices . Singapore: World Scientific,2003.[Bur09] S. Burer. On the copositive representation of binary and continuous nonconvex quadratic pro-grams.
Math. Program. , 120(2):479–495, 2009.[CEK14] C. Chau, K. Elbassioni, and M. Khonji. Truthful mechanisms for combinatorial ac electric powerallocation. In
AAMAS , pages 1005–1012, 2014.[Cha93] B. Chazelle. An optimal convex hull algorithm in any fixed dimension.
Discrete and Computa-tional Geometry , 10(1):377–409, 1993.[CHW76] A. Chandra, D. Hirschberg, and C. Wong. Approximate algorithms for some generalized knap-sack problems.
Theor. Comput. Sci. , 3(3):293–304, 1976.[DG14] P. Dickinson and L. Gijben. On the computational complexity of membership problems for thecompletely positive cone and its dual.
Computational Optimization and Applications , 57(2):403–415, 2014.[Dug09] S. Dughmi. Submodular functions: Extensions, distributions, and algorithms. a survey.
CoRR ,abs/0912.0322, 2009.[EKP02] T. Erlebach, H. Kellerer, and U. Pferschy. Approximating multiobjective knapsack problems.
Management Science , 48(12):1603–1612, 2002.[FC84] A. Frieze and M. Clarke. Approximation algorithms for the m -dimensional 0-1 knapsack prob-lem: worst-case and probabilistic analyses. European Journal of Operational Research , 15:100–109, 1984.[GHS80] G. Gallo, P. Hammer, and B. Simeone. Quadratic kanpsack problems.
Math. Program. , 12:132–149, 1980.[GKR11] V. Goyal, L. Kaya, and R. Ravi. An fptas for minimizing the product of two non-negative linearcost functions.
Math. Program. , 126(2):401–405, 2011.[GLS88] M. Gr¨otschel, L. Lov´asz, and A. Schrijver.
Geometric Algorithms and Combinatorial Optimiza-tion . Springer, New York, 1988.[GR13] V. Goyal and R. Ravi. An fptas for minimizing a class of low-rank quasi-concave functions overa convex set.
Oper. Res. Lett. , 41(2):191–196, 2013.24GV88] D. Grigoriev and N. Vorobjov. Solving systems of polynomial inequalities in subexponentialtime.
J. Symb. Comput. , 5(1/2):37–64, 1988.[Hal86] M. Hall.
Combinatorial Theory (2nd Ed.) . John Wiley & Sons, Inc., New York, NY, USA, 1986.[HXZ +
11] Z. He, S. Xie, R. Zdunek, G. Zhou, and A. Cichocki. Symmetric nonnegative matrix factor-ization: Algorithms and applications to probabilistic clustering.
IEEE Transactions on NeuralNetworks , 22(12):2117–2131, 2011.[IK75] O. Ibarra and C. Kim. Fast approximation algorithms for the knapsack and sum of subset prob-lems.
J. ACM , 22(4):463–468, 1975.[IP01] R. Impagliazzo and R. Paturi. On the complexity of k-sat.
J. Comput. Syst. Sci. , 62(2):367–375,2001.[JW02] D. Rader Jr. and G. Woeginger. The quadratic 0-1 knapsack problem with series-parallel support.
Oper. Res. Lett. , 30(3):159–166, 2002.[KN07] J. Kelner and E. Nikolova. On the hardness and smoothed complexity of quasi-concave mini-mization. In
FOCS , pages 472–482, 2007.[KS10] H. Kellerer and V. Strusevich. Fully polynomial approximation schemes for a symmetricquadratic knapsack problem and its scheduling applications.
Algorithmica , 57(4):769–795, 2010.[KS12] H. Kellerer and Vitaly A. Strusevich. The symmetric quadratic knapsack problem: approxima-tion and scheduling applications. , 10(2):111–161, 2012.[KST09] A. Kulik, H. Shachnai, and T. Tamir. Maximizing submodular set functions subject to multiplelinear constraints. In
SODA , pages 545–554, 2009.[KW07] W. Kern and G. Woeginger. Quadratic programming and combinatorial minimum weight productproblems.
Math. Program. , 110(3):641–649, 2007.[Lau70] D. Laughhunn. Quadratic binary programming with application to capital-budgeting problems.
Operations Research , 18:454–461, 1970.[LMNS09] J. Lee, V. Mirrokni, V. Nagarajan, and M. Sviridenko. Non-monotone submodular maximizationunder matroid and knapsack constraints. In
STOC , pages 323–332, 2009.[McM70] P. McMullen. The maximum numbers of faces of a convex polytope.
Mathematika , 17(02):179–184, 1970.[Moi13] A. Moitra. An almost optimal algorithm for computing nonnegative rank. In
SODA , pages1454–1464, 2013.[MS13a] S. Mittal and A. Schulz. An fptas for optimizing a class of low-rank functions over a polytope.
Math. Program. , 141(1-2):103–120, 2013.[MS13b] S. Mittal and A. Schulz. A general framework for designing approximation schemes for combi-natorial optimization problems with many objectives combined into one.
Operations Research ,61(2):386–397, 2013.[MY80] R. Mcbride and J. Yormark. An implicit enumeration algorithm for quadratic integer program-ming.
Manage. Sci. , 26(3):282296, 1980.[NT08] A. Nemirovski and M. Todd. Interior-point methods for optimization.
Acta Numerica , 17:191–234, 5 2008. 25Pis07] D. Pisinger. The quadratic knapsack problem - a survey.
Discrete Applied Mathematics ,155(5):623–648, 2007.[PS13] U. Pferschy and J. Schauer. Approximating the quadratic knapsack problem on special graphclasses. In
WAOA , pages 61–72, 2013.[PY00] C. Papadimitriou and M. Yannakakis. On the approximability of trade-offs and optimal accessof web sources. In
FOCS , pages 86–92, 2000.[Ren92] J. Renegar. On the computational complexity and geometry of the first-order theory of the reals.
J. Symb. Comput. , 13(3):255–352, 1992.[Sch86] A. Schrijver.
Theory of Linear and Integer Programming . Wiley, New York, 1986.[SS03] S. Schaible and J. Shi. Fractional programming: the sum-of-ratios case.
Optimization Methodsand Software , 18:219–229, 2003.[Svi04] M. Sviridenko. A note on maximizing a submodular set function subject to a knapsack constraint.
Oper. Res. Lett. , 32(1):41–43, 2004.[Vav09] S. Vavasis. On the complexity of nonnegative matrix factorization.
SIAM Journal on Optimiza-tion , 20(3):1364–1377, 2009.[VCZ11] J. Vondr´ak, C. Chekuri, and R. Zenklusen. Submodular function maximization via the multilinearrelaxation and contention resolution schemes. In
STOC , pages 783–792, 2011.[VHVD08] A. Vandendorpe, D. Ho, S. Vanduffel, and P. Dooren. On the parameterization of the creditrisk + model for estimating credit portfolio risk. Insurance: Mathematics and Economics , 42(2):736–745, 2008.[Von08] J. Vondr´ak. Optimal approximation for the submodular welfare problem in the value oraclemodel. In
STOC , pages 67–74, 2008.[Xu12] Z. Xu. A strongly polynomial fptas for the symmetric quadratic knapsack problem.
EuropeanJournal of Operational Research , 218(2):377–381, 2012.[YC13] L. Yu and C. Chau. Complex-demand knapsack problems and incentives in ac power systems.In
AAMAS , pages 973–980, 2013.[ZS05] R. Zass and A. Shashua. A unifying approach to hard and probabilistic clustering. In