Numerical Simulation of a Quadratic MIP with a highly rational term

Question

I am interested in solving the following minimization problem: $$ \begin{array}{cl} \displaystyle\min_{x,y}&\displaystyle\frac{1}{K}\sum_{i=1}^{K}\left(\frac{x_{i}}{y_{i}}-\frac{X}{Y}\right)^{2} \\[0.4cm] \mathrm{s.t.}&\sum_{i}x_{i}=X \\[0.2cm] &\sum_{i}y_{i}=Y \\[0.2cm] &x_{i}\leqslant y_{i},\quad i=1,2,\ldots,K \\[0.2cm] &x_{i},y_{i}\in\mathbb{Z}^{+},\quad i=1,2,\ldots, K \end{array} \tag{1} $$ with $X$ and $Y$ being fixed.

I attempted to simulate $(1)$ using $\textsf{cvxpy}$. However, it did not work as the quadratic term $\frac{x_{i}^{2}}{y_{i}^{2}}$ is neither convex, nor concave, and $\textsf{cvxpy}$ has no framework for such cases.

Question: How can we reformulate this problem so that it can be simulated by $\textsf{cvx/cvxpy}$, or what QMIP solver is able to handle an objective function of this form being highly rational and quadratic?

The solution to this minimization problem all boils down to $x_{i}=\lceil y_{i}\frac{X}{Y}\rceil$ or something close to that. For what its worth, I attempted to reformulate this problem by adding an envelop through $f_{i}:=\frac{x_{i}}{y_{i}}$. This leads to the following minimization problem: $$ \begin{array}{cl} \displaystyle\min_{f}&\displaystyle\frac{1}{K}\sum_{i=1}^{K}\left(f_{i}-\frac{X}{Y}\right)^{2} \\[0.4cm] \mathrm{s.t.}&y_{i}f_{i}-x_{i}=0 \\[0.2cm] &\sum_{i}x_{i}=X \\[0.2cm] &\sum_{i}y_{i}=Y \\[0.2cm] &f_{i}\leqslant 1,\quad i=1,2,\ldots,K \\[0.2cm] &x_{i},y_{i}\in\mathbb{Z}^{+},\quad i=1,2,\ldots, K \end{array} \tag{2} $$ However this formulation in $\textsf{cvxpy}$ was to no avail as there is no guided way to handle cases with bilinear constraints.

Answer 1

Getting an exact solution via brute force

I'll try to circle back later to formulate this as a MIP, but your problem as stated is small enough that you can just brute force the solution.

For instance, for K=3, X=941, Y=1223, the solution is

Best val = 8.22518e-12
xi =  307  317  317 
yi =  399  412  412 

Using 12 cores this takes 8.5 minutes to find using the code below.

This code finds said solution:

// Compile with: g++ -O3 42886-quadratic-mip.cpp -Wall -Wextra -fopenmp
#include <cassert>
#include <iostream>
#include <limits>
#include <iomanip>
#include <tuple>
#include <vector>

void generate_sum_sets_helper(
  const int number_of_vars,
  const int X,
  const bool allow_zeros,
  const int Xtotal,
  std::vector<int>& running_sol,
  const int running_sum,
  std::vector<std::vector<int>>& solutions
){
  if(X < 0){
    return;
  } else if(number_of_vars==0){
    if(running_sum==Xtotal){
      solutions.push_back(running_sol);
    }
  } else if(X==0){
    if(allow_zeros && running_sum==Xtotal){
      running_sol.push_back(0);
      generate_sum_sets_helper(number_of_vars-1, X, allow_zeros, Xtotal, running_sol, running_sum, solutions);
      running_sol.pop_back();
    }
  } else {
    for(int i=(allow_zeros?0:1);i<=X;i++){
      running_sol.push_back(i);
      generate_sum_sets_helper(number_of_vars-1,X-i,allow_zeros,Xtotal,running_sol,running_sum+i, solutions);
      running_sol.pop_back();
    }
  }
}

std::vector<std::vector<int>> generate_sum_sets(
  const int number_of_vars,
  const int X,
  const bool allow_zeros
){
  std::vector<int> running_sol;
  std::vector<std::vector<int>> solutions;
  generate_sum_sets_helper(number_of_vars, X, allow_zeros, X, running_sol, 0, solutions);
  return solutions;
}

double objective(
  const std::vector<int>& xvals,
  const std::vector<int>& yvals,
  const int X,
  const int Y
){
  assert(xvals.size()==yvals.size());
  double sum = 0;
  for(size_t i=0;i<xvals.size();i++){
    const auto temp = xvals[i] / static_cast<double>(yvals[i]) - X / static_cast<double>(Y);
    sum += temp * temp;
  }
  return sum / xvals.size();
}

bool all_less(
  const std::vector<int>& xvals,
  const std::vector<int>& yvals
){
  assert(xvals.size()==yvals.size());
  for(size_t i=0;i<xvals.size();i++){
    if(xvals[i]>=yvals[i]){
      return false;
    }
  }
  return true;
}

int main(){
  // Using prime numbers to prevent an exact solution
  constexpr auto Xmax =  941;
  constexpr auto Ymax = 1223;
  constexpr auto K = 3;
  static_assert(Xmax < Ymax);
  std::cerr<<"Generating xvals..."<<std::endl;
  const auto xvals = generate_sum_sets(K, Xmax, true);
  std::cerr<<"Generated "<< xvals.size() << " xvals..."<<std::endl;
  std::cerr<<"Generating yvals..."<<std::endl;
  const auto yvals = generate_sum_sets(K, Ymax, false);
  std::cerr<<"Generated "<< yvals.size() << " yvals..."<<std::endl;
  std::tuple<double, std::vector<int>, std::vector<int>> bestval = {std::numeric_limits<double>::infinity(), {}, {}};
  auto best_overall = bestval;
  std::cerr<<"Searching for solutions..."<<std::endl;
  #pragma omp parallel firstprivate(bestval)
  {
    #pragma omp for collapse(2)
    for(size_t xi=0;xi<xvals.size();xi++){
      for(size_t yi=0;yi<yvals.size();yi++){
        if(!all_less(xvals.at(xi),yvals.at(yi))){
          continue;
        }
        const auto oval = objective(xvals.at(xi),yvals.at(yi),Xmax, Ymax);
        if(oval < std::get<0>(bestval)){
          std::cerr<<"Best objective = "<<std::setprecision(10)<<oval<<std::endl;
          bestval = std::make_tuple(oval, xvals.at(xi), yvals.at(yi));
        }
      }
    }

    #pragma omp critical
    {
      if(std::get<0>(bestval) < std::get<0>(best_overall)){
        best_overall = bestval;
      }
    }
  }
  std::cout<<"Best val = "<<std::get<0>(best_overall)<<std::endl;
  std::cout<<"xi = ";
  for(const auto &x: std::get<1>(best_overall)){
    std::cout<<std::setw(4)<<x<<" ";
  }
  std::cout<<std::endl;
  std::cout<<"yi = ";
  for(const auto &x: std::get<2>(best_overall)){
    std::cout<<std::setw(4)<<x<<" ";
  }
  std::cout<<std::endl;

  return 0;
}

cvxpy doesn't work

In cvxpy you would write your envelope formulation as:

import cvxpy as cp

# Create variables
K = 3
X = 941
Y = 1223

# Add variables
xi = cp.Variable(K, integer=True)
yi = cp.Variable(K, integer=True)
fi = cp.Variable(K)

# Add constraints
constraints = [
    sum(xi) == X,
    sum(yi) == Y,
    xi >= 0,
    yi >= 1,
    xi <= X,
    yi <= Y,
    fi >= 0,
    fi <= 1,
]
for i in range(K):
    constraints.append(xi[i] <= yi[i]-1)
for i in range(K):
    constraints.append(yi[i]*fi[i] - xi[i] == 0)

# Add objective
objective_expr = cp.Minimize(sum([(fi[i] - X / Y)**2 for i in range(K)]) / K)

# Optimize the model
model = cp.Problem(objective_expr, constraints)
optval = model.solve()
print(optval)

but this doesn't work because the problem isn't convex:

cvxpy.error.DCPError: Problem does not follow DCP rules. Specifically:
The following constraints are not DCP:
var2[0] @ var3[0] + -var1[0] == 0.0 , because the following subexpressions are not:
|--  var2[0] @ var3[0]
var2[1] @ var3[1] + -var1[1] == 0.0 , because the following subexpressions are not:
|--  var2[1] @ var3[1]
var2[2] @ var3[2] + -var1[2] == 0.0 , because the following subexpressions are not:
|--  var2[2] @ var3[2]

Getting a maybe-good solution using MIPs in Gurobi

So, instead, we switch to using Gurobi, which is free with an academic license and otherwise will cost a blood sacrifice beneath a full moon. Note that you could try using alternatives like XPRESS.

In Gurobi the problem is formulated as follows:

import gurobipy as gp
from gurobipy import GRB

# Create a new model
model = gp.Model("qp")

#Set parameters
model.params.NonConvex = 2
model.setParam('FeasibilityTol', 1e-9)
model.setParam('OptimalityTol', 1e-9)
model.setParam('MIPGap', 1e-13)
model.setParam('IntFeasTol', 1e-9)
model.setParam('NumericFocus', 3) # VERY VERY IMPORTANT

# Create variables
K = 3
X = 941
Y = 1223

# Add variables
xi = model.addVars(K, lb=0, ub=X, vtype=GRB.INTEGER, name="x")
yi = model.addVars(K, lb=1, ub=Y, vtype=GRB.INTEGER, name="y")
fi = model.addVars(K, lb=0, ub=1, vtype=GRB.CONTINUOUS, name="f")

# Add constraints
model.addConstr(sum(xi[i] for i in range(K)) == X, name="constraint_sum_x")
model.addConstr(sum(yi[i] for i in range(K)) == Y, name="constraint_sum_y")
for i in range(K):
    model.addConstr(xi[i] <= yi[i]-1, name=f"constraint_xgy_{i+1}")
for i in range(K):
    model.addConstr(yi[i]*fi[i] - xi[i] == 0, name=f"constraint_f{i}")

# Add objective
objective_expr = sum((fi[i] - X / Y)**2 for i in range(K)) / K
model.setObjective(objective_expr, GRB.MINIMIZE)

# Optimize the model
model.optimize()

# Get the solution
if model.status == GRB.OPTIMAL:
    print(f"min_value = {model.objVal}")
    print(f"x_values = {[xi[i].x for i in range(K)]}")
    print(f"y_values = {[yi[i].x for i in range(K)]}")
    # Double check that we got the objective right
    print(f"Objective: ", sum([(xi[i].x/yi[i].x - X/Y)**2 for i in range(K)])/K)
else:
    print("Could not solve!")

model.printStats()
model.printQuality()

Note that there's a line I've labeled as VERY VERY IMPORTANT.

If we remove that line after 1.62s Gurobi confidently gives us:

Optimal solution found (tolerance 1.00e-13)
Best objective 5.844680295297e-11, best bound 5.844680295297e-11, gap 0.0000%
min_value = 5.844680295297167e-11
x_values = [327.0, 297.0, 317.0]
y_values = [425.0, 386.0, 412.0]
Objective:  5.844692423821425e-11

Solution quality statistics for model qp :
  Maximum violation:
    Bound       : 0.00000000e+00
    Constraint  : 0.00000000e+00
    Integrality : 0.00000000e+00

This is wrong! We know this because we brute forced the solution on a small problem.

We need to tighten every tolerance to its tightest value and set numeric focus to get the right answer:

model.setParam('FeasibilityTol', 1e-9)
model.setParam('OptimalityTol', 1e-9)
model.setParam('MIPGap', 1e-13)
model.setParam('IntFeasTol', 1e-9)
model.setParam('NumericFocus', 3) # VERY VERY IMPORTANT

and after 0.39s Gurobi hands that to us:

Optimal solution found (tolerance 1.00e-13)
Best objective 8.225198300238e-12, best bound 8.225198300238e-12, gap 0.0000%
min_value = 8.225198300237935e-12
x_values = [317.0, 307.0, 317.0]
y_values = [412.0, 399.0, 412.0]
Objective:  8.225179646366615e-12

Solution quality statistics for model qp :
  Maximum violation:
    Bound       : 0.00000000e+00
    Constraint  : 0.00000000e+00
    Integrality : 0.00000000e+00

What have we learned?

What we've learned is that Gurobi is very fast in comparison to the brute force solution because it's able to use the problem's struture and its advanced heuristics to ignore enormous parts of the search space. We've also learned that it will confidently hand us a suboptimal answer because the problem you've formulated has many answers that are very close to optimality and Gurobi has difficulty distinguishing between them unless we really push it to be careful.

With large values of K I'd certainly become worried that even with the settings used above I might not be getting the very best answer. Of course, with large values of K a brute force solution will also take much longer.