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