Earth Mover's Distance (EMD)

Earth Mover's Distance (EMD) について調べたことを整理しておきます。EMDは、ユークリッド距離のような距離尺度の一つで、二つの分布の間の距離を測ることができます。言語処理ではあまり聞いたことなかったのですが、画像処理や音声処理では比較的有名な距離尺度のようです。

EMDが使える問題設定は下図のようになります。

EMDは特徴量と重みの集合（シグネチャと呼ぶ）で与えられる分布Pと分布Qの間の距離です。ここで、特徴量間では距離 が定義されているのが前提です。特徴量がベクトルのときはユークリッド距離、特徴量が確率分布のときはカルバック・ライブラー距離（情報量）などです。EMDは、特徴量の集合が2つ与えられたときに、1個1個の特徴量間の距離をもとに、特徴量集合間の距離を求められるんですね。これはすごい。

重みは具体的な応用によって使い方が変わりますが、その特徴量の重要度を表しています。たとえば、ヒストグラムだったら各棒が特徴量にあたり、棒の高さが重みにあたります。前に類似画像検索 (2009/10/3）で、画像の色のヒストグラムからHistogram Intersectionという距離を使いましたが、EMDを使って距離を求めることもできます。参考文献にあげたEMDの原論文は類似画像検索を対象にしています。

EMDなんて使わず、もっと単純に全特徴量のあらゆる組み合わせ間の距離の総和でもいいんじゃね？と思いましたけど、これだけだと重みを完全に無視していますね・・・重みが重要なんです！

考え方の基本は輸送問題

EMDの定義は、最適化問題の1つの輸送問題（Transportation Problem）の考え方に基づいています。なのでまずは輸送問題について簡単にまとめます。先の図において、Pの各場所P1, ... ,Pmには、重みの量だけ荷物が積まれているとします。そして、Qの各場所Q1, ... ,Qnには重みの量だけ格納できる倉庫があるとします。このとき、Pにある荷物をすべてQに運ぶ*1とき、どこからどこへどのくらい運ぶともっとも効率がよいかを求めるのが輸送問題です。

ここで、Pi から Qj へ輸送するコスト（距離）を とし、Pi から Qj へ輸送する荷物量（フロー）を と定義します。そして、Pi から Qj へ運ぶのに要する仕事量を と定義します。たとえば、距離が遠いところに大量の荷物を運ぶとそれだけ仕事量が増えるので直感とも一致します。このとき、総仕事量Wを下のように定義すると、W を最小化する を求めればもっとも効率のよい運び方だとわかります。

は与えられるのが前提なので、最適化する変数は輸送量 だけです。そして、輸送量 には下の4つの制約が加えられます。

(1) かならずPからQへ輸送する。逆方向はない。

(2) Piにある荷物以上は輸送できない

(3) Qjにある倉庫の容量以上は荷物を受け付けられない

(4) 輸送量の上限は、荷物の総量か倉庫の総容量の小さい方

最後の条件は荷物の総量と倉庫の容量が違うときに必要です。荷物総量より倉庫の総容量が大きかったら全部輸送できるので輸送量の上限は荷物の総量となりますが、荷物が倉庫の量より多かったら全部輸送できないので輸送量の上限は倉庫の総容量になります。今回、取り上げる例題は荷物の総量と倉庫の総容量は同じとしています。

輸送問題の解き方は省略しますが、解くと最適な が求まります。EMDはこの を用いて下のように定義されます。輸送量 の合計で割っているのは、輸送量によってEMDのスケールが変わらないように正規化しているためですね。これはあとで具体例で確認してみます。

EMDは輸送に必要な最適な仕事量が小さいほど二つのシグネチャの距離は近いという考え方なのでわりと自然な考え方だと思います。ただし、あらゆる特徴量間の組み合わせについて足し合わせが必要なので特徴量の数が多くなると計算量は非常に大きくなりそうです。そのため、特徴量の数が爆発しないようにベクトル量子化と組み合わせてシグネチャを作る手法が提案されています。これは、後に紹介予定です。

EMDの定義がわかったところで具体例を解いてみます

この例は、EMDの提案者のRubnerさんのライブラリに出てくる例題です。特徴量は3次元ベクトルで重みは浮動小数点数で与えられています。特徴量の値が0から255の3次元ベクトルなので、分布Pが画像1のカラーヒストグラム、分布Qが画像2のカラーヒストグラムを表しているようです。この分布Pと分布QのEMDを計算してみます！

RubnerのC言語実装

まずは、Rubnerさんが公開されているC言語のコードを使ってみます（example1.c）。実行には、emd.cとemd.hが必要です。また、emd.hのfeature_tの定義を問題に合わせて書き換える必要があります。今回は、特徴量が3次元ベクトルなので

typedef struct { int X,Y,Z; } feature_t;

と定義しています。emd.cのライブラリを使って上の例題を解くコードです。

#include <stdio.h>
#include <math.h>
#include "emd.h"

/* ユークリッド距離 */
float dist(feature_t *F1, feature_t *F2) {
    int dX = F1->X - F2->X;
    int dY = F1->Y - F2->Y;
    int dZ = F1->Z - F2->Z;
    return sqrt(dX*dX + dY*dY + dZ*dZ);
}

int main() {
    /* 分布Pの特徴ベクトル */
    feature_t f1[4] = { {100,40,22}, {211,20,2}, {32,190,150}, {2,100,100} };
    /* 分布Qの特徴ベクトル */
     feature_t f2[3] = { {0,0,0}, {50,100,80}, {255,255,255} };
    /* 分布Pの重み */
    float w1[5] = { 0.4, 0.3, 0.2, 0.1 };
    /* 分布Qの重み */
    float w2[3] = { 0.5, 0.3, 0.2 };
    /* 分布Pのシグネチャ */
    signature_t s1 = { 4, f1, w1 };
    /* 分布Qのシグネチャ */
     signature_t s2 = { 3, f2, w2};

    /* EMDを計算 */
    float e;
    e = emd(&s1, &s2, dist, 0, 0);
    printf("emd = %f\n", e);
    return 0;
}

emd()関数に2つの分布のシグネチャと特徴量間の距離を計算する関数を指定しています。この実装では、重みは浮動小数点になっています。合計するとどちらも1.0になるので荷物総量と倉庫の総容量は等しくなっています。実行すると、

emd = 160.542770

よって、分布Pと分布Qの距離は、160.54とわかります。先に書いたようにEMDは、輸送量で正規化しているため下のように比率を保ったまま重みを変えても結果は同じになります。

    /* 分布Pの重み */
    float w1[5] = { 4.0, 3.0, 2.0, 1.0 };
    /* 分布Qの重み */
    float w2[3] = { 5.0, 3.0, 2.0 };

R言語による実装

次は、R言語で同じ例題を解いてみます。R言語で輸送問題を解く関数は、lpSolveというライブラリに含まれていますが、標準では入っていないのでインストールします。

install.packages("lpSolve")

以下のemd_sample.Rファイルを作成します。

library(lpSolve)

# ユークリッド距離
euclid_dist <- function(f1, f2) {
    return(sqrt(sum((f1 - f2)^2)))
}

# EMDを計算
emd <- function(dist, w1, w2) {
    # lp.transport()を使うための準備
    costs <- dist
    row.signs <- rep("<", length(w1))
    row.rhs <- w1
    col.signs <- rep(">", length(w2))
    col.rhs <- w2
    # 輸送問題を解く
    t <- lp.transport(costs, "min", row.signs, row.rhs, col.signs, col.rhs)
    # 最適な輸送量を取得
    flow <- t$solution
    # 仕事量を計算
    work <- sum(flow * dist)
    # 正規化してEMDを計算
    e <- work / sum(flow)
    return(e)
}

# 特徴量
f1 = matrix(c(100, 40, 22, 211, 20, 2, 32, 190, 150, 2, 100, 100), 4, 3, byrow=T)
f2 = matrix(c(0, 0, 0, 50, 100, 80, 255, 255, 255), 3, 3, byrow=T)

# 重み（要整数！）
w1 = c(4, 3, 2, 1)
w2 = c(5, 3, 2)

n1 = length(f1[,1])
n2 = length(f2[,1])

# 距離行列を作成
dist = matrix(0, n1, n2)
for (i in 1:n1) {
    for (j in 1:n2) {
        dist[i, j] = euclid_dist(f1[i,], f2[j,])
    }
}

# 距離行列と重みからEMDを計算
e = emd(dist, w1, w2)
cat(sprintf("emd = %f\n", e))

この実装では、シグネチャを渡さずに、シグネチャから計算した距離行列と重みを渡しています。lp.transport()が距離行列を受け付けるのでそれに合わせましたが、まあ、どっちでもよいと思います。Rを起動して以下のように打つと実行できます。

> source("emd_sample.R")
emd = 160.542763

C言語版と同じ結果になりました！ただ、距離行列の計算がRっぽい書き方でないので効率悪いかも。もっといい書き方があったら教えてください。

それでもPythonで書きたいんだよ！

PythonのSciPyにはlp.transport()に対応する関数はないようです。他の最適化ライブラリ（openoptやcvxopt）もざっと探しましたが見つけられませんでした。自分で書いてもよかったのですが、せっかくRの関数があるのでrpy2を使ってPythonからRのlp.transport()を呼び出してみます。

rpy2は、Rの機能をPythonから呼び出せるようにするPythonライブラリです。使い方はやや複雑ですが、今回のようにPythonになくてRにあるアルゴリズムもたくさんあるので使えるとはかどります。まあ、Rをラッパーしているのでちょっと複雑になっちゃうけどね。

#coding:utf-8
import numpy as np
import rpy2.robjects as robjects

# Rのlp.transport()をインポート
robjects.r['library']('lpSolve')
transport = robjects.r['lp.transport']

def euclid_dist(feature1, feature2):
    """ユークリッド距離を計算"""
    if len(feature1) != len(feature2):
        print "ERROR: calc euclid_dist: %d <=> %d" % (len(feature1), len(feature2))
        return -1
    return np.sqrt(np.sum((feature1 - feature2) ** 2))

def emd(dist, w1, w2):
    """Rのtransport()関数を使ってEMDを計算"""
    # transport()の引数を用意
    costs = robjects.r['matrix'](robjects.FloatVector(dist),
                                 nrow=len(w1), ncol=len(w2),
                                 byrow=True)
    row_signs = ["<"] * len(w1)
    row_rhs = robjects.FloatVector(w1)
    col_signs = [">"] * len(w2)
    col_rhs = robjects.FloatVector(w2)

    t = transport(costs, "min", row_signs, row_rhs, col_signs, col_rhs)
    flow = t.rx2('solution')

    dist = dist.reshape(len(w1), len(w2))
    flow = np.array(flow)
    work = np.sum(flow * dist)
    emd = work / np.sum(flow)
    return emd

if __name__ == "__main__":
    f1 = np.array([ [100, 40, 22], [211, 20, 2], [32, 190, 150], [2, 100, 100] ])
    f2 = np.array([ [0, 0, 0], [50, 100, 80], [255, 255, 255] ])

    # 重みは自然数のみ
    w1 = np.array([4, 3, 2, 1])
    w2 = np.array([5, 3, 2])

    n1 = len(f1)
    n2 = len(f2)

    # 距離行列を作成
    dist = np.zeros(n1 * n2)
    for i in range(n1):
        for j in range(n2):
            dist[i * n2 + j] = euclid_dist(f1[i], f2[j])

    # 距離行列と重みからEMDを計算
    print "emd =", emd(dist, w1, w2)

実行すると、

> python emd_sample.py
emd = 160.542762808

となり、CやRと同じ結果が得られました。

おわりに

今回は、EMDについて調べたことをまとめてみました。できるだけ正確に書くようにしましたが、誤りもあるかもしれません。鵜呑みにしないで他の資料も当たってみてください。

あとで、EMDを類似楽曲検索という具体的な問題に応用してみる予定です。乞うご期待。

参考文献

• Earth mover's distance - Wikipedia
• Y. Rubner, C. Tomasi and L. J. Guibas: The earth mover's distance as a metric for image retrieval (PDF), International Journal of Computer Vision, 40(2), pp.99-121, 2000 - EMDの原論文。EMDを類似画像検索に適用しています。
• Code for the Earth Movers Distance (EMD) - Rubnerさんが公開されているC言語実装
• Fast Earth Mover's Distance (EMD) Code - EMDを高速計算する実装
• 柳本, 大松: Earth Mover's Distanceを用いたテキスト分類、人工知能学会全国大会, 2007. - EMDの説明がわかりやすい。画像や音声の手法がテキストにも使えるんですね。
• lpSolve - R言語のlpSolveのマニュアル。lp.transform()の詳しい仕様はここで。

*1:実際はすべてではなくQに入る分だけでOK