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)

参考

Wasserstein metricによる解法

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

from scipy.stats import wasserstein_distance

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

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

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

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

参考