Integer programming problem with simple quadratic objective function in Python

Question

I have $n$ objects that need to be divided among $k$ groups. Each group must receive at least $5$ objects. In addition, the percentage of objects in group $i$ should be as close as possible to $p_i$ where $\sum\limits_{i=1}^k p_i = 1$. I was thinking of formulating this as an integer programming problem with the number of objects in group $i$ denoted by $h_i$. The problem then becomes:

\begin{align}\min&\quad\sum\limits_{i=1}^k\left(\frac{h_i}{n}-p_i\right)^2\\\text{s.t.}&\quad h_i \geq 5\\&\quad\sum_i h_i = n\end{align}

And of course, all the $h_i$ are integers. I was looking at solving this in Python since the production environment I need to deploy to has this setup. I looked into integer programming for Python but can only find linear programs so far. Any pointers?

Answer 1

MOSEK can handle large-scale problems of a wide variety, including conic and Mixed-Integer programming problems. The Fusion API provided by MOSEK makes it straightforward to utilize the MOSEK solver in a structured manner.

In the problem that you state, one can rephrase the quadratic objective function by introducing an auxiliary scalar variable $H$, as follows: \begin{align}\text{min}&\quad H\\\text{s.t.}&\quad H \geq \sqrt{\sum_{i=1}^k \bigg( \frac{h_i}{n} - p_i\bigg)^{2}}\end{align}

while the other constraints are unchanged. I have expressed the objective in terms of the epigraph of the Euclidean norm of the vector $\vec{d} = \vec{(h/n)} - \vec{p}$ that can be implemented exactly in Fusion as a $k+1$ dimensional quadratic cone, i.e. $(H,\vec{d})\in Q^{k+1}$.

The Fusion implementation for your problem could be as follows:

import numpy as np
from mosek.fusion import *

#n is the number of objects and p is the vector of ratios (size k)
def quad_int_prog(n,p):
    #Number of groups
    k = len(p)
    #Making a Mosek model
    with Model('k_groups') as M:
        #h (vector of size k): Integer and h>=5
        h = M.variable('h',k,Domain.integral(Domain.greaterThan(5))) 
        #Auxiliary Variable (epigraph)
        H = M.variable('H',Domain.greaterThan(0))

        #Linear Constraint: sum(h_i)=n
        M.constraint('h_sum_constraint',Expr.sum(h),Domain.equalsTo(n))
        #Conic Constraint: (H,d) belongs to Quadratic cone domain of size k+1 
        H_d_vector = Expr.vstack(H,Expr.sub(Expr.mul(1/n,h),p))
        M.constraint('Percentage_constraint',H_d_vector,Domain.inQCone(k+1))

        #Objective: Minimize H 
        M.objective('Objective',ObjectiveSense.Minimize,H)
        #Solving... (To enable log output, use M.setLogHandler(sys.stdout))
        M.solve()
        #Optimal value for h
        h_optimal = h.level()
    return(h_optimal)

Being a large-scale solver, MOSEK can easily handle scaled up versions of the discussed model. For instance, on my desktop (with i7-6700HQ and 16 GB of RAM) I could solve a problem of the above-stated structure with $k=3000$ and $n=20000$ (with a randomly generated $p$) in about 4.23 seconds with a relative tolerance gap of less than $2 \%$. Remember, it can be crucial to set a termination criteria when dealing with Mixed-Integer problems, such as setting a tolerance gap or providing a maximum number of iterations. Moreover, I highly recommend enabling the log-output in MOSEK (using the command M.setLogHandler(sys.stdout)) because it not only tells you the progress the solver has made, but can also provide key insights into the model and the structure of your problem.

If you need a quick reference guide for mathematical optimization and general model building (like the rephrasing presented above using the Euclidean norm), check out the MOSEK modelling cookbook. For a comprehensive guide to using Fusion and/or other products that MOSEK offers, see here. It is possible to acquire a free trial license (valid for 30 days), and a free academic license (valid for 1 year) is also provided.

Answer 2

As the $h_i$ are integers, you can reformulate your problem as a Mixed Integer Linear Program. Here's how:

Let $c_i^m = \left(\frac{m}{n} - p_i\right)^2$ denote the penalty that is received when assigning $m$ items to group $i$. Now introduce binary decision variables $x_i^m$ that are equal to 1 if $m$ items are assigned to group $i$ and 0 otherwise. As you have a constraint stating that at least 5 items must be assigned to every group and there are $k$ groups, you need to define these variables for $5\leq m \leq n - 5(k-1)$ (EDIT: you can probably find much better bounds. For example, suppose that $\lfloor{p_in\rfloor}\geq 5$ for all $i$, then it is never optimal to assign less than $\lfloor{p_in\rfloor}$ objects to group $i$.)

Then, your problem is equivalent to

\begin{align} \min \quad &\sum_{i=1}^k \sum_{m=5}^{n-5(k-1)} c_i^mx_i^m, \\ \text{s.t.} \quad &\sum_{m=5}^{n-5(k-1)} x_i^m =1 \quad&&\forall i, \\ &\sum_{i=1}^k \sum_{m=5}^{n-5(k-1)} mx_i^m = n, \\ &x_i^m \in \{0,1\} \quad &&\forall i,m. \end{align}

The first constraint makes sure that for each group, only one number of items can be assigned. The second constraint guarantees that in total $n$ items are assigned. Given that your problem has such simple structure (I imagine a rounding approach would perform well), I'm pretty sure this is not the most efficient way to solve your problem, but at least you can implement it easily as a MIP in Python.

Answer 3

I finally found cvxpy, which seems to work for me. Not sure if it has integer programming, but for now I just used non-integer and then split the fractional parts back among the solution array. Not ideal, but better than nothing. Code below. If someone has a solution that involves actual integer programming using an open-source library, please post an answer.

import cvxpy as cp
import numpy as np

## Number of objects.
n=50

## Required proportions (array length should be n and sum to 1).
p = np.array([.4,.2,.15,.1,.05,.05,.025,.025])

x = cp.Variable(len(p))
constraints = [4 <= x, x <= n, sum(x) == n]
objective = cp.Minimize(cp.sum_squares(x/n - p))
problm = cp.Problem(objective, constraints)

res = problm.solve()

vals = x.value // 1
excess = int(sum(x.value % 1))

excess_vals_unif = np.ones(len(p))* excess//len(p)

excess_vals_nonunif = np.concatenate((np.ones(excess % len(p)),\
                np.zeros(len(p)-excess%len(p)))\
                ,axis=0)

final_vals = vals + excess_vals_unif + excess_vals_nonunif
print(final_vals)

Answer 4

Your problem is an Integer Quadratic Knapsack Problem. Your specific problem is considered by Bretthauer et al. The book by Kellerer et al. is an excellent source for information on knapsack problems in general.

Modeling

Bretthauer et al. consider the following problem: $$\begin{eqnarray} \min & & \frac{1}{2} x^\top D x - a^\top x\\ s.t. & & \sum_{i=1}^n p_i x_i = c\\ & & l_i \le x_i \le u_i & \forall i = 1 \dots, n\\ & & x_i \in \mathbb{N}_{\ge 0} & \forall i = 1 \dots, n \end{eqnarray}$$ with $D$ a diagonal matrix with positive components.

Note that $$\left(\frac{h_i}{n} - p_i\right)^2 = \frac{1}{2} \left(\frac{2}{n^2}\right) h_i^2 - \left(\frac{2 p_i}{n}\right) h_i + p_i^2$$ In the optimization, we can ignore the constant $p_ i^2$ and we can define $D_{ii} = \frac{2}{n^2}$ and $a_i = \frac{2 p_i}{n}$ to arrive at the formulation of Bretthauer et al.

Branch and bound

Bretthauer et al. solve this problem with branch and bound. This paper was published in 1995, and in the meantime many different algorithms have been published. This paper, and the book by Kellerer et al., however, could serve as a starting point.

Dynamic programming

There is a well-known dynamic programming formulation for knapsack problems. Dynamic programming is discussed in the book by Kellerer et al., and also on the knapsack Wikipedia page. Because the optimal $h_i$ and $h_j$ are only related through the constraint, I see no reason why this dynamic program could not be adapted to your problem. The only preprocessing that is necessary is that you set all $h_i = 5$ as the starting value, and adjust the remaining capacity accordingly.

This would be different if $D$ is not diagonal. In this case, we would have terms like $D_{ij} h_i h_j$ in the objective, and the choice of $h_i$ changes how $h_j$ influences the objective.

Implementing the dynamic program is relatively straightforward, but it may not be fast enough for your problem size. Or it might, I simply don't know!