我目前正在使用 Python 中的单纯形算法来实现最优传输问题，用于原始和对偶。但是，我没有得到相同的原始值和对偶值。
我认为问题可能来自对偶，因为我还使用 Sinkhorn 算法解决了原始问题，它返回的值与基于单纯形的原始解决方案相同。
但是，我真的找不到哪里出了问题，而且我对问题可能来自哪里已经没有想法了。这是我的代码，我希望它很清楚！
感谢您抽出宝贵时间阅读，并提前感谢您提供的任何帮助！！！
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)
