最佳运输代码Scipy线性编程优化需要更长的时间

如何解决最佳运输代码Scipy线性编程优化需要更长的时间

我一直在尝试计算平均值为0.0和4.0且方差分别为9.0和16.0的两个一维高斯分布之间的Wasserstein距离。我使用了scipy.linprog.optimize模块，并使用了以下链接中所述的“ interior-point”方法 https://yetanothermathprogrammingconsultant.blogspot.com/2019/10/scipy-linear-programming-large-but-easy.html。但是，它花费了超过17个小时的时间，并且仍在运行（我的代码正在运行）以解决300 x 300 LP矩阵问题（即300个源节点和300个目标节点）。但是，该文件说有可能用1000个源节点和1000个目的节点来解决问题。（即，一个人可以用1,000,000（一百万）个决定性变量来解决LP问题。我的代码有什么问题？为什么要花这么长时间？我们是否需要大容量内存（或群集）来解决此类问题？我的代码

from datetime import datetime
start_time = datetime.Now()
from scipy.optimize import linprog
import scipy

#Initializing the LP matrix
Piprob=np.zeros(500*500).reshape(500,500)

def Piprobmin(Krv,rhoi,rhoj):
  r1=np.shape(Krv)[0]
  r2=np.shape(Krv)[1]
  print("r1,r2",r1,r2)
  #Computing the LP Matrix which has just two ones in each column
  pmat=np.zeros((r1+r2)*(r1*r2)).reshape((r1+r2),(r1*r2))
  for i in range(r1+r2):
      for j in range(r1*r2):
          if((i<r1) and (j<((i+1)*r2)) and (j>=(i*r2))):
              pmat[i][j]=1
          if(i>=r1):
              for k in range(r1*r2):
                  if j==(i-r1)+(k*r2):
                      pmat[i][j]=1
  #flattening the cost matrix into one dimensional array
  krvf=Krv.flatten()
  tempr=np.append(rhoi,rhoj)
  Xv=[] #Creating list for joint probability matrix elements
                          
  res = scipy.optimize.linprog(c=krvf,method='interior-point',A_eq=pmat,b_eq=tempr,options= 
  {'sparse':True,'disp':True})

  print("res=\n",res)
  wv=res.fun
  for  l1 in range(r1*r2):
      Xv.append(res.x[l1])
  Yv=np.array(Xv)
  Yv=Yv.reshape(r1,r2)
  #returning Yv-joint probability and,Wv-minimized wasserstein distance
  return Yv,wv

  Piprob,W=Piprobmin(K,result1,result2) #K-cost function matrix,result1 is the first 
                                        #marginal,result2 is the second marginal
  end_time = datetime.Now()
  print('Duration: {}'.format(end_time - start_time))

成本函数的大小为300 X 300，大小，每个边际都有300点（总共600个约束）。我确认我的成本函数是对称且非负的。并且每个边际之和都等于一个概率。

解决方法

在blog post中，单词 sparse 经常使用。不是没有原因的。将A矩阵存储为稀疏矩阵非常重要。否则，您将无法处理大问题。该博客文章详细讨论了运输LP矩阵在内存要求方面的差异，因此这一点应该不容错过。

这是一些示例代码，说明如何使用 scipy.optimize.linprog 建立具有1000个源节点和1000个目标节点的运输模型。同样，LP矩阵具有2,000行和1,000,000列，并且存储稀疏。

import numpy as np
import scipy as sp
import scipy.sparse as sparse
import scipy.optimize as opt
from memory_profiler import profile


def GenerateData(M,N):
   np.random.seed(123)

   # form objective function
   c = np.random.uniform(0,10,(M,N))

   # demand,supply
   s = np.random.uniform(0,15,M)
   d = np.random.uniform(0,N)
   assert np.sum(d) <= np.sum(s),"supply too small"

   #print('c',c)
   #print('s',s)
   #print('d',d)
   return {'c':c,'s':s,'d':d,'n':N,'m':M}

def FormLPData(data):

   rhs = np.append(data['s'],-data['d'])

   # form A
   # column (i,j)=n*i+j has two nonzeroes:
   #    1 at row i with rhs supply(i)
   #    1 at row N+j with rhs demand(j)
   N = data['n']
   M = data['m']
   NZ = 2*N*M
   irow = np.zeros(NZ,dtype=int)
   jcol = np.zeros(NZ,dtype=int)
   value = np.zeros(NZ)
   for i in range(N):
      for j in range(M):
         k = M*i+j
         k1 = 2*k
         k2 = k1+1
         irow[k1] = i
         jcol[k1] = k
         value[k1] = 1.0
         irow[k2] = N+j
         jcol[k2] = k
         value[k2] = -1.0

   A = sparse.coo_matrix((value,(irow,jcol)))

   #print('A',A)
   #print('rhs',rhs)

   return {'A':A,'rhs':rhs}



@profile
def run():
   # dimensions
   M = 1000  # sources
   N = 1000 # destinations
   data = GenerateData(M,N)
   lpdata = FormLPData(data)
   res = opt.linprog(c=np.reshape(data['c'],M*N),A_ub=lpdata['A'],b_ub=lpdata['rhs'],options={'sparse':True,'disp':True})


if __name__ == '__main__':
   run()

因此，您似乎完全错过了有关博客文章的全部内容。