Earth Mover's Distance (EMD)距离

原文： http://d.hatena.ne.jp/aidiary/20120804/1344058475
作者： sylvan5
翻译： Myautsai和他的朋友们(Google Translate、shuanger、qiu)

本文将讨论Earth Mover’s Distance (EMD)，和欧式距离一样，它们都是一种距离度量的定义、可以用来测量某两个分布之间的距离。EMD主要应用在图像处理和语音信号处理领域，在自然语言处理上很少有听说。

EMD 问题如下图所示

给定两个签名(或者叫分布、特征量集合)P和Q，P为m个特征量Pi和其权重wPi的集合，记作P={(P1,wP1),(P2,wP2),...(Pm,wPm)}，即上图左侧部分。同样的，还有一个分布Q，Q=(Q1,wQ1),(Q2,wQ2),...(Qn,wQn)，即上图右侧部分。在计算这个这两个签名的Earth Mover's Distance(EMD)前，我们要先定义好P、Q中任意取一个特征量( Pi and Qj )之间的距离dij（这个距离叫ground distance，两个签名之间EMD依赖于签名中特征量之间的ground distance）。当这两个特征量是向量时得到的dij是欧式距离，当这两个特征量是概率分布时得到的dij是相对熵(KL距离/Kullback–Leibler divergence)。现在，给定两个签名(P和Q)，只要计算好每两个特征量之间的距离，系统就能给出这两个签名之间的距离了。so easy！

不同情况下EMD使用方式也不一样，但还是有一些共通之处。比如权重都是指特征量的重要程度。例如，一个直方图对应一个签名的情况下，直方图中的每一根柱(bar)代表一个特征量，柱的高度就对应其权重。在之前的相似图像检索 (2009/10/3）一文中，我使用到了图像颜色分布直方图相交距离(Histogram Intersection )，也可以用在EMD中当作ground distance使用。最早提出EMD概念的论文中有提到，EMD最初就是用来做相似图片检索的。

运输问题概述

EMD 实际上是线性规划中运输问题的最优解。首先，简要描述下运输问题。我们假设这个例子是从多个工厂运输货物到多个仓库。在上图左侧，P从在P1 到 Pm代表m座工厂，工厂Pi有重量为wPi的货物。在上图右侧，Q从Q1到Qn代表n个仓库，仓库Qj最大容量为wQj。货物之间没有什么区别，都是同一类东西。每个仓库都希望装尽可能多的货物。如何尽可能高效把所有货物(实际上不一定是所有货物，部分也OK)从P运送到Q，就是运输问题的优化目标。在本例中，P、Q都是离散的，那么EMD可以用运输问题的Hungarian算法来计算它们之间的距离。挖个坑而已，这里不具体讨论。

这里定义，货物从工厂Pi运到仓库Qj，距离是dij，运送货物的重量为fij。这样一次运输需要的工作量为dij∗fij。显然，距离越远、或货物越重，工作量就越大。(注：运输可能是多对多的，即一个工厂运输货物到多个仓库，或者一个仓库接收多个工厂的货物。)货物从工厂运到仓库需要很多次这样的运输，经过一些计算和优化，这时我们得到了工作量总和的最小值W。

W=∑i=1m∑j=1ndijfij→min

距离 dij 是事先定义的，所以运输量 fij 是式中唯一的变量，对 fij 作如下4个约束：

(1) 运输过程从工厂P到仓库Q，不能反向。

fij≥0(1≤i≤m,1≤j≤n)

(2) 从工厂Pi一次次运出去的所有货物重量的和不可能超过该工厂中所有货物的总重量wPi

∑i=1nfij≤wPi(1≤i≤m)

(3) 仓库Qj接收的所有货物总重量不可能超过其总容量wQj。

∑i=1mfij≤wQj(1≤j≤n)

(4) 总运输量的上限是工厂中货物总重、仓库总容量中的最小值。

∑i=1m∑j=1nfij=min(∑i=1mwPi,∑j=1nwQj)

当仓库的总容量和工厂货物总重不一样时，我们才需要上述第4个限制条件。仓库总容量比货物总量大的话就可以全部运输了，所以这时候呢，运输量的上限就是货物总量。但如果货物总量比仓库总容量大的话，就不能全部运输了，这时候，运输量的上限就是仓库的总容量啦。方便起见，本例中当仓库的总容量和工厂货物总重是一样的。
运输问题的具体解答此处省略不讨论，假设这个时候我们已经得到了最优解f∗ij。为了使EMD不会随着总运输量的变化而变化，每一次的运输量还要除以总运输量，以达到归一化的目的。（in order to avoid favoring smaller signatures in the case of partial matching.）在后面的具体例子中会对它进行详细描述。

EMD(P,Q)=∑mi=1∑nj=1dijf∗ij∑mi=1∑nj=1f∗ij

很自然可以想到，给定两个签名，把一个变成另一个所需要的最小工作量，就是EMD对距离的定义，这里的「工作量」要基于用户对ground distance的定义，即特征量之间的距离的定义。然而，当特征量非常多的时候，由于要做一一匹配，其计算量是非常大的。因此，有人提出了一种将多个特征量组合起来做向量量化编码(Vector Quantization)后再组成签名的方法。

EMD的一些优点可见这里

举个栗子

这个例子是EMD提出者Rubner在其C语言实现的库中提到的。上图左侧给出的分布P是一些加权后的特征量，代表图像1的颜色直方图，因为R(ed)G(reen)B(lue)三原色，所以维度是3，某个颜色的多少用权重表示，比如左上角第一个特征量P1=[100,40,22]，其权重wP1=0.4，意思是R=100,G=40,B=22这个颜色在图片1中占的比重是0.4。同样地，分布Q代表图像2。现在我们尝试计算分布P与分布Q之间的EMD！

Rubner的C语言实现

首先我们尝试使用Rubner桑公开的C语言代码(example1.c)，编译依赖emd.c和emd.h。其中特征量类型feature_t在emd.h中定义如下：

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

具体实现代码见emd.c。对于上述例子的解答如下：

给定的dist()函数用于计算特征量间的距离，在此基础上emd()函数计算指定两个分布的签名之间距离。上面提到「本例中当仓库的总容量和工厂货物总重是一样的」，所以P、Q各自权重和都是1.0。运行结果：

emd = 160.542770

可得P、Q两个分布的签名之间的距离是160.54。只要归一化后保持P、Q两个总权重之间的比例，就算改变了权重值结果也是不变的。

/*分布P的权重 */
float w1[5] = { 0.4, 0.3, 0.2, 0.1 };
/*分布Q的权重 */
float w2[3] = { 0.5, 0.3, 0.2 };

R语言的实现

接下来，我们尝试使用R语言来解决同样的问题。R的lpSolve函数可以直接用来解决运输问题，但这个函数并不在标准库中，得先安装下：

install.packages("lpSolve")

创建一个如下内容的emd_sample.R文件

上述代码中，emd的计算没有直接使用两个签名，而是用两个签名之间的距离矩阵和权重。因为根据lp.transport()的定义中传入的参数需要是一个距离行列式。运行结果如下：

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

我用C语言实现的版本也得到了一样的结果，但是效率和R相比差，毕竟后者是通过距离矩阵计算得到的。如果你知道怎样写更好的话请告诉我。

还是用Python来实现吧！

Python的Scipy中似乎没有lp.transport()的对应函数，搜了下其他库（openopt和cvxopt）也没有特别合适的实现，只能蛋疼地在Python里用rpy2调用R的lp.transport()函数。rpy2是一个让你可从python里调用R的python包。看着有些复杂，因为调用了一个封装在R中的函数。

当运行

> python emd_sample.py
emd = 160.542762808

在C语言或R语言的实现中，我们也得到了同样的结果。

结束语

本文对与EMD的讨论力求准确，但是错误难免，敬请批评指正，同时请参考其他文献。

参考文献

Earth mover’s distance - Wikipedia <a href="http://en.wikipedia.org/wiki/Earth_mover" s_distance"="" style="border: 0px; margin: 0px; padding: 0px; vertical-align: baseline; color: rgb(116, 51, 153); background: transparent;">link
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言語実装link
Fast Earth Mover’s Distance (EMD) Code - EMDを高速計算する実装 link
柳本, 大松: Earth Mover’s Distanceを用いたテキスト分類、人工知能学会全国大会, 2007. - EMDの説明がわかりやすい。画像や音声の手法がテキストにも使えるんですね。
lpSolve - R言語のlpSolveのマニュアル。lp.transform()の詳しい仕様はここで。
Formal definition of EMD

this article is mainly based on the original text written by sylvan5 on aidiary.some additional contents are added by mckelvin.
本文主要基于sylvan5发表在aidiary的原文，在此基础上增加了一些内容。

from: http://blog.mckelv.in/articles/1474.html