python - Python – Simplex为原始和对偶返回不同的值(Python – Simplex returning different values for the primal and the dual)

Question

I'm currently working on an implementation of the optimal transport problem using the simplex algorithm in Python, for the primal and the dual.

(我目前正在使用Python中的单纯形算法针对原始和对偶实现最佳运输问题。)

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.

(我认为问题可能来自于对偶问题，因为我也使用Sinkhorn的算法解决了原始问题，并且它返回的值与基于单纯形的原始解决方案相同。)

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)

  ask by MKpls translate from so

Answer 1

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

(运输问题（相等版本）的原始和对偶看起来像：)

When using bounds=None in the call to linprog, we tell the solver to use non-negative variables.

(在对linprog的调用中使用bounds=None时，我们告诉求解器使用非负变量。)

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.

(对于双重问题，您需要使用bounds=(None,None)表示变量应该是自由的。)

Now you should see that the optimal objective values for the primal and the dual problem are identical.

(现在您应该看到，原始问题和对偶问题的最佳目标值是相同的。)