Directly calling gurobipy API causes substantially longer runtime than calling cvxpy

Question

Background

I am trying to implement the reranking algorithm proposed in the paper (Eqn 3). The algorithm is cast into an integer linear program which is trying to find permutation matrix $\mathbf{X} \in \mathbb{R}^{n\times n}$, where $\mathbf{X}_{ij} =1$ means ranking item $i$ to position $j$.

Question

Initially, I implemented this algorithm using cvxpy and setting the solver to be Gurobi. The runtime is not satisfying. So I decided to directly call Gurobi Python API (i.e. gurobipy) rather than using the modeling interface (i.e. cvxpy). I expected the runtime to be improved as I skip the abstraction from cvxpy. However, the runtime is surprisingly much worse (see below).

# gurobipy runtime in seconds
174.727278
# cvxpy runtime in seconds
54.925721

The only two explanations I could think are

• My implementation based on gurobipy is not efficient.

• cvxpy somehow optimizes the formulation by using some transformation. This transformation makes it easier for Gurobi to discover the structure in the problem. However, when directly code the formulation using gurobipy, this structure is not well exploited by Gurobi.

My questions are

• Does any of my explanations make sense? Why this odd behavior might happen.
• What could I do to optimize my cvxpy code or gurobipy code to improve the runtime?

Code to reproduce my result

import numpy as np
import cvxpy as cp
import gurobipy as gp

from datetime import datetime
from gurobipy import GRB

def runGUROBIImpl(n, theta=0.5, nIter=10):
    for i in range(nIter):
        with gp.Env(empty=True) as env:
            env.setParam("OutputFlag", 0)
            env.start()
            with gp.Model(env=env) as model:
                C = np.random.rand(n, n)
                DCG = np.random.rand(n, n)
                IDCG = np.random.rand()

                P = model.addVars(n, n, vtype=GRB.BINARY)
                obj = sum([C[i, j] * P[i, j] for i in range(n) for j in range(n)])
                constrLHS = sum([DCG[i, j] * P[i, j] for i in range(n) for j in range(n)])

                model.setObjective(obj, GRB.MINIMIZE)
                model.addConstr(constrLHS >= theta * IDCG)
                model.addConstrs((P.sum("*", i) == 1 for i in range(n)))
                model.addConstrs((P.sum(i, "*") == 1 for i in range(n)))

                model.optimize()

                if not (model.status == GRB.OPTIMAL): print("unsuccessful...")

def runCVXPYImpl(n, theta=0.5, nIter=10):
    C, DCG, IDCG = cp.Parameter((n, n)), cp.Parameter((n, n)), cp.Parameter()
    P = cp.Variable((n, n), boolean=True)

    objective = cp.Minimize(cp.sum(cp.multiply(C, P)))
    constraints = [cp.sum(cp.multiply(DCG, P)) >= theta * IDCG,
                   cp.matmul(np.ones((1, n)), P) == np.ones((1, n)),
                   cp.matmul(P, np.ones((n, 1))) == np.ones((n, 1))]
    problem = cp.Problem(objective, constraints)

    for i in range(nIter):
        C.value = np.random.rand(n, n)
        DCG.value = np.random.rand(n, n)
        IDCG.value = np.random.rand()
        # solve problem
        problem.solve(solver=cp.GUROBI, verbose=False)

        if not (problem.status == cp.OPTIMAL): print("unsuccessful...")

np.random.seed(0)
start = datetime.now()
runGUROBIImpl(n=100, theta=0.5, nIter=300)
end = datetime.now()
print((end - start).total_seconds())

np.random.seed(0)
start = datetime.now()
runCVXPYImpl(n=100, theta=0.5, nIter=300)
end = datetime.now()
print((end - start).total_seconds())

Logs

The algorithm proposed in the paper requires the linear program to be solved several hundreds times with updated parameter matrix each time to get the final result $\mathbf{X}$.

As it would be much too verbose to provide logs for several hundreds solving attempts, the following are the logs when I set nIter=1 for both runGUROBIImpl() and runCVXPYImpl().

Note when nIter=1 the situation for runtime is different from when nIter=300 as the cvxpy requires a lot of time to construct problem.

# gurobipy runtime in seconds
0.61062
# cvxpy runtime in seconds
15.098772

runGUROBIImpl() Logs

Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (linux64)
Thread count: 2 physical cores, 4 logical processors, using up to 4 threads
Optimize a model with 201 rows, 10000 columns and 30000 nonzeros
Model fingerprint: 0x4962970c
Variable types: 0 continuous, 10000 integer (10000 binary)
Coefficient statistics:
  Matrix range     [2e-04, 1e+00]
  Objective range  [7e-05, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [2e-01, 1e+00]
Found heuristic solution: objective 48.7825807
Presolve removed 1 rows and 0 columns
Presolve time: 0.03s
Presolved: 200 rows, 10000 columns, 20000 nonzeros
Variable types: 0 continuous, 10000 integer (10000 binary)

Root relaxation: objective 1.415358e+00, 210 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

*    0     0               0       1.4153579    1.41536  0.00%     -    0s

Explored 0 nodes (210 simplex iterations) in 0.03 seconds
Thread count was 4 (of 4 available processors)

Solution count 2: 1.41536 48.7826 

Optimal solution found (tolerance 1.00e-04)
Best objective 1.415357907906e+00, best bound 1.415357907906e+00, gap 0.0000%

runCVXPYImpl() Logs

Parameter OutputFlag unchanged
   Value: 1  Min: 0  Max: 1  Default: 1
Changed value of parameter QCPDual to 1
   Prev: 0  Min: 0  Max: 1  Default: 0
Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (linux64)
Thread count: 2 physical cores, 4 logical processors, using up to 4 threads
Optimize a model with 201 rows, 10000 columns and 30000 nonzeros
Model fingerprint: 0x8b5dd50a
Variable types: 0 continuous, 10000 integer (10000 binary)
Coefficient statistics:
  Matrix range     [2e-04, 1e+00]
  Objective range  [7e-05, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [2e-01, 1e+00]
Found heuristic solution: objective 47.5478288
Presolve removed 1 rows and 0 columns
Presolve time: 0.03s
Presolved: 200 rows, 10000 columns, 20000 nonzeros
Variable types: 0 continuous, 10000 integer (10000 binary)

Root relaxation: objective 1.415358e+00, 210 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

*    0     0               0       1.4153579    1.41536  0.00%     -    0s

Explored 0 nodes (210 simplex iterations) in 0.03 seconds
Thread count was 4 (of 4 available processors)

Solution count 2: 1.41536 47.5478 

Optimal solution found (tolerance 1.00e-04)
Best objective 1.415357907906e+00, best bound 1.415357907906e+00, gap 0.0000%

Answer 1

As others have pointed out already, you are not solving the same instances. When writing out the MPS files using a gurobi.env file containing GURO_PAR_DUMP=1, we can see that the instances differ (here n=2, nIters=1; left is cvxpy, right is Gurobi):

To get the signs in order, you could change this line in the Gurobi method

model.addConstr(constrLHS >= theta * IDCG)

to

model.addConstr(-constrLHS <= -theta * IDCG)

Apparently, cvxpy changes this internally. This still leaves the variable order different in both models. I didn't manage to get this identical for both codes. I did verify that the only difference is in the ordering - the objective values coincide.

That Gurobi takes more time in some configurations (especially when the instance size is larger) is probably due to the fact that you are constructing a new model every time (you should reuse the environment, by the way). In cvxpy on the other hand, the model structure is only constructed once and then reused for all instances. This also doesn't require creating a new set of variables every time. I am not sure how this is handled internally, though, since at some point the model data needs to be handed over to Gurobi.

Answer 2

Note that the model fingerprint differs. I suspect that the variable or constraint orders are different.