Calculation

Pythonで計算するEarth Mover's Distance

Mathematics Python SciPy

更新日

2020年12月11日

作成日

2020年12月12日

EMDとは

Earth mover's distance (EMD) とは統計学において二つの確率分布間の距離を測る指標です．数学では Wasserstein metric として知られているそうです．

特徴点の分布を $P=\{(p_1,w_{p_1}),\cdots,(p_m,w_{p_m})\}$ ， $Q=\{(q_1,w_{q_1}),\cdots,(q_n,w_{q_n})\}$ としたとき

\min \sum_{i=1}^{m} \sum_{j=1}^{n} f_{i,j} d_{i,j}

where $f_{i,j}$ is a flow from $p_i$ to $q_j$ , and $d_{i,j}$ is a (euclidean) distance.

EMD(P, Q) = \frac{\sum_{i=1}^{m} \sum_{j=1}^{n} f^{*}_{i,j} d_{i,j}}{\sum_{i=1}^{m} \sum_{j=1}^{n} f_{i,j}}

where $f^{*}_{i,j}$ is an optimal flow, and

\sum_{i=1}^{m} \sum_{j=1}^{n} f_{i,j} = \min \left( \sum_{i=1}^{m} w_{p_i}, \sum_{j=1}^{n} w_{q_j} \right)

と求めることができます（＋制約式）．要は，最適な経路を通るときのコストです．

参考

線形割当問題による解法

EMD は確率分布間の最小距離を測る手法で．原理的には複数の倉庫から複数の店舗へ商品を納品するときに，最も配送コストの低い配送ルートを選択するのと同じです．

この問題は，線形計画法において割当問題と呼ばれています．

与えられているのは距離とそれに対応するコストで，コストが等しければ最短経路になります．

線形割当問題は SciPy のパッケージを用いることで解くことができます．

sample.py

import numpy as np
from scipy.optimize import linear_sum_assignment
from scipy.spatial.distance import cdist

# 分布上の特徴点の数
n = 10
x1 = np.random.randint(0, 9, n)
x2 = np.random.randint(0, 9, n)

# ユークリッド距離
d = cdist(x1[:,np.newaxis], x2[:,np.newaxis])
# 線形割当問題の解
assignment = linear_sum_assignment(d)
# コスト
emd = d[assignment].sum() / n
print(emd)

cdist() はデフォルトでは二つの分布の特徴点間のそれぞれのユークリッド距離を計算します．

上記は各特徴点が1次元ですが，下記のように $l$ 次元に拡張することもできます．

sample.py

n = 10
l = 3
x1 = np.random.randint(0, 9, (n, l))
x2 = np.random.randint(0, 9, (n, l))

d = cdist(x1, x2)
assignment = linear_sum_assignment(d)
emd = d[assignment].sum() / n
print(emd)

参考

Wasserstein metricによる解法

一次元の分布間の距離であれば wasserstein_distance() でも計算できます．

sample.py

from scipy.stats import wasserstein_distance

wd = wasserstein_distance(x1, x2)
print(wd)

グラフにおける最小コストフローによる解法

二つの分布間ということは，グラフ理論における完全2部グラフで同じ問題を考えることができます．

NetworkXに用意されている min_cost_flow_cost() を使うとネットワークフローの最小化を解くことができるようです．