Search code examples
pythonscipylinear-programmingsimplex

Python – Simplex returning different values for the primal and the dual


I’m currently working on an implementation of the optimal transport problem using the simplex algorithm in Python, for the primal and the dual. However, I don't get the same values for the primal and the dual.

I think the problem might come from the dual, because I have also solved the primal problem using Sinkhorn's algorithm, and it returned the same value as the simplex-based primal solution.

However, I really can't find what is wrong and I am running out of ideas as to where the problem might come from. Here is my code, I hope it is clear !

Thank you for taking the time to read, and thank you in advance for any help that might come !!!

import numpy as np
from random import random
import scipy
from scipy import optimize

#descriptions of the primal (resp dual) problem line 42 (resp 48 and 50)

n = 10 #number of points in distribution x_a
m = 20 #number of points in distribution x_b

x_a = [random() for i in range(n)] #n points chosen randomly between 0 and 1
x_a.sort() #sorted list
x_b = [2*random() for j in range(m)] #m points chosen randomly between 0 and 2
x_b.sort()

a = [1/n for i in range(n)] #distribution vector for x_a (uniform)
b = [1/m for i in range(m)] #distribution vector for x_b (uniform)
B = a+b #concatenation of a and b

#cost function (quadratic cost)
def cost(x,y) :
    n = len(x)
    p = len(y)
    tab = [[0 for j in range(p)] for i in range(n)]
    for i in range(n):
        for j in range(p):
            tab[i][j] = (x[i]-y[j])**2
    return tab

#definition of the matrix A
a0 = np.kron([1 for i in range(m)], np.eye(n, n))
a1 = np.kron(np.eye(m, m), [1 for i in range(n)])
A = np.concatenate((a0,a1), axis = 0)

#matrix of cost which we put columnwise into a 1D vector
cost_matrix = cost(x_a,x_b)
cost_column = []
for j in range(0,m):
    for i in range(n):
        cost_column += [cost_matrix[i][j]]

#Primal problem : Minimize cost_column*p under the constraint A*p=B (p decision variable)
proba = scipy.optimize.linprog(cost_column, A_ub=None, b_ub=None, A_eq=A, b_eq=B, bounds=None, method='simplex', callback=None, options=None, x0=None)
objective_result_primal = proba.fun
print(objective_result_primal)

A_transpose = np.transpose(A)
#Dual problem : Maximize B*h under the constraint A_transpose*h<=cost_column
B2 = [-B[i] for i in range(m+n)]#we use the opposite of B to turn the maximization problem into a minimization one (so that we can use the simplex)
#The dual problem becomes : Minimize B2*h under the constraint A_transpose*h<=cost_column
proba = scipy.optimize.linprog(B2, A_ub=A_transpose, b_ub=cost_column, A_eq=None, b_eq=None, bounds=None, method='simplex', callback=None, options=None, x0=None)
objective_result_dual = - proba.fun
print(objective_result_dual)

Solution

  • The primal and dual for the (equality version of the) transportation problem look like:

    enter image description here

    When using bounds=None in the call to linprog, we tell the solver to use non-negative variables. This is correct for the primal, but not for the dual problem. For the dual problem you need to use bounds=(None,None) which indicates the variables should be free. Now you should see that the optimal objective values for the primal and the dual problem are identical.