Comparison to Gurobi: OT (emd) not giving optimal plan

[floyebolu]

Hi,

I'm not sure if this is a "discussion" or "issue"...

This may be a bit naive but I thought the emd solver provided an exact solution. I have a series of transportation problems I need to solve over the course of a simulation run and I went to POT based on suggestion from a colleague and because it has a nice/easy interface.

I get a much worse answer using ot.emd than if I create the ILP and solve it in Gurobi.

Please see my code below.

import numpy as np
import gurobipy as gp
import ot

# Load data
first_day = np.load('first_day_arrays.npz')

cost_matrix = first_day['cost_matrix'].T  # cost matrix with shape (3001, 20)
units_hist = first_day['units_hist']  # supply vector with shape (3001,)
reqs_hist = first_day['reqs_hist']  # demand vector with shape (20,)

# Python Optimal Transport
ot_plan = ot.emd(units_hist, reqs_hist, cost_matrix)

# Gurobi
pre_one_first = np.ones(cost_matrix.shape[0])
post_one_first = np.ones(cost_matrix.shape[1])[:, None]

gurobi_model = gp.Model("first_day")
gp_plan = gurobi_model.addMVar(shape=cost_matrix.shape, lb=0, vtype=gp.GRB.INTEGER, name="x")
gurobi_model.setObjective(np.kron(pre_one_first, post_one_first.T) @ (gp_plan.reshape(-1) * cost_matrix.reshape(-1)), sense=gp.GRB.MINIMIZE)
gurobi_model.addConstr(gp_plan @ post_one_first == units_hist[:, None])
gurobi_model.addConstr(pre_one_first @ gp_plan == reqs_hist)
gurobi_model.optimize()

# Compare transport costs
print('\n\n')
print(f"Python Optimal Transport cost: {np.sum(ot_plan * cost_matrix)}")
# Result: 4.737712175590463
print(f"Gurobi cost: {gurobi_model.objVal}")  # np.sum(gp_plan.X * cost_matrix) also works
# Result: 0.7318552570921428

# Compare plans
print('\n\n')
print(f'Shape of Python Optimal Transport plan: {ot_plan.shape}') # (3001, 20)
print(f'Shape of Gurobi plan: {gp_plan.X.shape}') # (3001, 20)
print(f'Are the plans equal? {np.allclose(ot_plan, gp_plan.X)}')
# Result: False

Is there something I'm missing about POT (e.g., optimality not guaranteed) or is something not working as it should?
I even changed the numItermax parameter in the ot.emd function from its default to 10_000_000 and up to 150_000_000, but the result did not change.

scratch.zip
Attached within scratch.zip:

• first_day_arrays.npz
• gp_plan.X.npz (in case one is not able to access/use Gurobi. Retrieve like so: gp_plan_X = np.load('gp_plan.X.npz')['gp_plan_X'])

Thanks for your help!

P.S.: Gurobi is free for academics, and if you pip install gurobipy one is able to solve "small" optimisation problems.

[floyebolu]

Based on a suggestion by @rflamary, I converted all the input data to np.float64 but I got the same result(s) :/

[cedricvincentcuaz]

Hello @floyebolu,
Thank you for pointing this out. Could you let us know whether or not the transport plan outputed by gurobi is admissible (satisfy well the marginal constraints)?
And if gurobi_model.objVal coincides with the operation used to compute the OT loss using numpy?

[floyebolu]

Hi @cedricvincentcuaz

The transport plan from Gurobi meets the constraints.
For your second question, I assume you mean that does gurobi_model.objVal equal the product of its plan and the cost matrix. Which the answer to is yes also.

supplied = gp_plan.X.sum(axis=1)
received = gp_plan.X.sum(axis=0)
print(supplied.shape) # (3001,)
print(received.shape) # (20,)

np.allclose(supplied, units_hist)
# Result: True

np.allclose(received, reqs_hist)
# Result: True


np.allclose(gurobi_model.ObjVal, np.sum(gp_plan.X * cost_matrix))
# Result: True

Thanks for all your help :)

[rflamary]

I had a look at the cost matrix and it appears that you have very large values and very small ones. I guess that gurobi has tricks to handle this but the OT solver in ot.emd uses float64 that can handle a dynamic of only around 1e16 and provides the best solution at this precision.

I just cliped the max values à 1e10 with cost_matrix2 = np.clip(cost_matrix, a_min=-np.inf, a_max=1e10) and the resulting OT plan and cost is now very close to gurobi at 0.7318552570921427 when computing the loss with the original matrix.

[floyebolu]

Thank you for answering, @rflamary

After consulting with a colleague and trying another solver, I am using another way to resolve this.

I have different large values signifying varying degrees of 'bad' matches in my transportation problem. E.g., 1e32 is to show that it is impossible to move between two points, and 1e16 is to allow for shortages (a unit is moved from a dummy supply point to demand point). So clipping at a maximum value is not suitable for my problem as I need to distinguish between these allocations in more challenging problem instances.

What I have done instead is to transform the cost matrix by multiplying the 'smaller' values by a constant. Then I've raised the 'largest' values to a fractional power to reduce them whilst maintaining relative size differences between these large cost/penalty values. This way the range of values has been squeezed: the maximum cost value is bounded at ~ 6e12.

The following code snippet should clarify what I mean.

# previous code setup and imports
# ...

cost_matrix = first_day['cost_matrix'].T  # cost matrix with shape (3001, 20)

# Scaling the matrix
orig_cost_matrix = cost_matrix * 1
large_indx = np.abs(cost_matrix) > 1e12
small_indx = np.abs(cost_matrix) <= 1
cost_matrix[large_indx] = np.power(cost_matrix[large_indx], 2/5)
cost_matrix[small_indx] *= 1000

With this I get a resulting cost of 0.7318552570921428.

Thanks for all the help :)