Pythonで計算するEarth Mover's Distance
EMDとは
Earth mover's distance (EMD) とは統計学において二つの確率分布間の距離を測る指標です．数学では Wasserstein metric として知られているそうです．
特徴点の分布を $P={(p_1,w_),\cdots,(p_m,w_)}$，$Q={(q_1,w_),\cdots,(q_n,w_)}$ としたとき
$$ \min \sum_^{m} \sum_^{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_^{m} \sum_^{n} f^{*}{i,j} d{i,j}}{\sum_^{m} \sum_^{n} f_{i,j}} $$
where $f^{*}_{i,j}$ is an optimal flow, and
$$ \sum_^{m} \sum_^{n} f_{i,j} = \min \left( \sum_^{m} w_, \sum_^{n} w_ \right) $$
と求めることができます（＋制約式）．要は，最適な経路を通るときのコストです．
参考
線形割当問題による解法
EMD は確率分布間の最小距離を測る手法で．原理的には複数の倉庫から複数の店舗へ商品を納品するときに，最も配送コストの低い配送ルートを選択するのと同じです．
この問題は，線形計画法において割当問題と呼ばれています．
与えられているのは距離とそれに対応するコストで，コストが等しければ最短経路になります．
線形割当問題は SciPy のパッケージを用いることで解くことができます．
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$次元に拡張することもできます．
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)
参考
- scipy.spatial.distance.cdist | SciPy v1.5.4 Reference Guide
- scipy.optimize.linear_sum_assignment | SciPy v1.5.4 Reference Guide
Wasserstein metricによる解法
一次元の分布間の距離であれば wasserstein_distance() でも計算できます．
from scipy.stats import wasserstein_distance
wd = wasserstein_distance(x1, x2)
print(wd)
グラフにおける最小コストフローによる解法
二つの分布間ということは，グラフ理論における完全2部グラフで同じ問題を考えることができます．
NetworkXに用意されている min_cost_flow_cost() を使うとネットワークフローの最小化を解くことができるようです．