Блог

Я пытался вычислить расстояние Вассерштейна между двумя одномерными распределениями Гаусса со средним значением 0,0 и 4,0, с отклонениями 9,0 и 16,0 соответственно. Я использовал модуль scipy.linprog.optimize и использовал метод «внутренней точки», как указано в следующей ссылке https://yetanothermathprogrammingconsultant.blogspot.com/2019/10/scipy-linear-programming-large-but-easy.html . Тем не менее, требуется более 17 часов, и все еще (мой код) выполняется для решения задач матрицы 300 x 300 LP (т. Е. 300 исходных узлов и 300 узлов назначения. Однако в документе говорится, что проблему можно решить с 1000 исходными узлами и 1000 узлами назначения. (т. е.) проблему LP можно решить с помощью 1 000 000 (одного миллиона) решающих переменных. Что не так с моим кодом? Почему это занимает так много времени? Нужна ли нам большая память (или кластеры) для решения таких проблем? мой код

 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 ограничений). Я проверил, что моя функция затрат симметрична и неотрицательна. и каждое значение суммируется с единицей, поскольку это просто вероятности.

В сообщении в блоге слово sparse используется много раз. Не без причины. Чрезвычайно важно хранить матрицу A как разреженную матрицу. В противном случае вы не сможете справиться с большими проблемами. В сообщении в блоге очень подробно обсуждается разница в требованиях к памяти транспортной матрицы LP, поэтому этот момент было трудно пропустить.

Вот несколько примеров кода о том, как настроить транспортную модель с 1000 исходными узлами и 1000 узлами назначения с помощью scipy.optimize.linprog. Опять же, матрица LP имеет 2000 строк и 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,10,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()
  

Похоже, вы полностью упустили суть сообщения в блоге.