Python Earth Mover Расстояние 2D массивов

Question

Я хотел бы рассчитать расстояние движителя Земли между двумя двумерными массивами (это не изображения).

Прямо сейчас я просматриваю две библиотеки: scipy ( https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.wasserstein_distance.html) и pyemd ( https://pypi.org/project/pyemd/).

#define a sampeling method
def sampeling2D(n, mu1, std1, mu2, std2):
   #sample from N(0, 1) in the 2D hyperspace
   x = np.random.randn(n, 2)

   #scale N(0, 1) -> N(mu, std)
   x[:,0] = (x[:,0]*std1) + mu1
   x[:,1] = (x[:,1]*std2) + mu2

   return x

#generate two sets
Y1 = sampeling2D(1000, 0, 1, 0, 1)
Y2 = sampeling2D(1000, -1, 1, -1, 1)

#compute the distance
distance = pyemd.emd_samples(Y1, Y2)

Хотя версия scipy не принимает 2D-массивы и возвращает ошибку, метод pyemd возвращает значение. Если вы видите из документации, он говорит, что он принимает только 1D массивы, поэтому я думаю, что вывод неправильный. Как я могу рассчитать это расстояние в этом случае?

Answer 1

Поэтому, если я вас правильно понимаю, вы пытаетесь перенести распределение выборки, то есть рассчитать расстояние для установки, где все кластеры имеют вес 1. В общем, вы можете рассматривать расчет EMD как пример минимального потока затрат, и в вашем случае это сводится к проблеме линейного присваивания: ваши два массива - это разбиения в двудольном графе, а веса между двумя вершинами - ваше расстояние выбора. Предполагая, что вы хотите использовать евклидову норму в качестве метрики, веса ребер, то есть расстояния до земли, можно получить с помощью scipy.spatial.distance.cdistи на самом деле SciPy предоставляет решатель для задачи линейного назначения суммы, а также в scipy.optimize.linear_sum_assignment (в котором недавно произошли огромные улучшения производительности, которые будут доступны в SciPy 1.4. Это может вас заинтересовать, если у вас возникнут проблемы с производительностью; реализация 1.3 немного медленная для 1000x1000 входов).

Другими словами, то, что вы хотите сделать, сводится к

from scipy.spatial.distance import cdist
from scipy.optimize import linear_sum_assignment

d = cdist(Y1, Y2)
assignment = linear_sum_assignment(d)
print(d[assignment].sum() / n)

Возможно, было бы полезно проверить, что результат этого расчета совпадает с тем, что вы получите от решателя с минимальными затратами; один такой решатель доступен в NetworkX, где мы можем построить график вручную:

import networkx as nx

G = nx.DiGraph()

# Represent elements in Y1 by 0, ..., 999, and elements in
# Y2 by 1000, ..., 1999.
for i in range(n):
    G.add_node(i, demand=-1)
    G.add_node(n + i, demand=1)

for i in range(n):
    for j in range(n):
        G.add_edge(i, n + j, capacity=1, weight=d[i, j])

На этом этапе мы можем убедиться, что описанный выше подход согласуется с минимальным потоком затрат:

In [16]: d[assignment].sum() == nx.algorithms.min_cost_flow_cost(G)
Out[16]: True

Точно так же полезно видеть, что результат согласуется с scipy.stats.wasserstein_distance для одномерных входов:

from scipy.stats import wasserstein_distance

np.random.seed(0)
n = 100

Y1 = np.random.randn(n)
Y2 = np.random.randn(n) - 2
d =  np.abs(Y1 - Y2.reshape((n, 1)))

assignment = linear_sum_assignment(d)
print(d[assignment].sum() / n)       # 1.9777950447866477
print(wasserstein_distance(Y1, Y2))  # 1.977795044786648