1 引言

最近看到一篇特征匹配相關(guān)的論文，思想是將特征匹配問題轉(zhuǎn)化為最優(yōu)傳輸問題求解，于是我去學(xué)習(xí)了一下最優(yōu)傳輸問題。
本文主要是對博文 Notes on Optimal Transport 的學(xué)習(xí)做一個(gè)記錄總結(jié)，該博文寫的不錯(cuò)，推薦閱讀。

2 例子：分甜點(diǎn)

文章作者以一個(gè)簡單的甜點(diǎn)分配例子引入了最優(yōu)傳輸問題。
向量 $\mathbf{r}=[3,3,3,4,2,2,2,1{]}^{?}$ 表示 $n=8$ 個(gè)人需要的甜點(diǎn)數(shù)：

向量 $\mathbf{c}=[4,2,6,4,4{]}^{?}$ 表示 $m=5$ 種甜點(diǎn)的數(shù)量：

矩陣 $\mathbf{M}\in {\mathbb{R}}^{5\times 8}$ 表示每個(gè)人對各種甜點(diǎn)的偏好，尺度區(qū)間$[?2,2]$，-2表示非常不喜歡，2表示非常喜歡：

我們的目標(biāo)，就是要根據(jù)甜點(diǎn)的數(shù)量，同時(shí)考慮每個(gè)人的需求和偏好，將所有甜點(diǎn)合理地分配到每個(gè)人手中。

3 最優(yōu)傳輸問題

最優(yōu)運(yùn)輸問題的目標(biāo)就是以最小的成本將一個(gè)概率分布轉(zhuǎn)換為另一個(gè)概率分布。上面的分甜點(diǎn)的目標(biāo)，用最優(yōu)傳輸問題的定義來說，就是將概率分布 $\mathbf{c}$ 以最小的成本轉(zhuǎn)換到概率分布 $\mathbf{r}$ 。
這就需要我們求得一個(gè)分配方案，由矩陣 $P\in {\mathbb{R}}^{n\times m}$ 表示，存儲每個(gè)人分得的每個(gè)甜點(diǎn)的情況。

根據(jù)現(xiàn)實(shí)條件，這個(gè)分配矩陣 $P$ 顯然具有以下約束：

1. 分配的甜點(diǎn)數(shù)量不能為負(fù)數(shù)；
2. 每個(gè)人的需求都要滿足，即 $P$ 的行和服從分布 $\mathbf{r}$；
3. 每種甜點(diǎn)要全部分完，即 $P$ 的列和服從分布 $\mathbf{c}$。

于是在分布 $\mathbf{r}$、$\mathbf{c}$ 約束下， $P$ 的解空間可以做如下定義：
$$U(\mathbf{r},\mathbf{c})=\left\{P\in {\mathbb{R}}_{>0}^{n\times m}\mid P{1}_m=\mathbf{r},P^{?}{1}_n=\mathbf{c}\right\}\text{\hspace{1em}(1)}$$
PS：這是博文的原公式，這里我有個(gè)疑問，為什么 $P$ 的元素要求嚴(yán)格大于0，而不是大于等于0？希望有同學(xué)能夠解答我的疑惑（感謝）

如前面所述，我們希望最小化轉(zhuǎn)換成本，可以簡單地反轉(zhuǎn)偏好矩陣 $\mathbf{M}$ 的符號，就可以得到成本矩陣（cost matrix）。于是就有了最優(yōu)傳輸問題的公式化表示：
$$d_M(\mathbf{r},\mathbf{c})=\underset{P\in U(\mathbf{r},\mathbf{c})}{\min }\sum _{i,j}{P}_{ij}{M}_{ij}\text{\hspace{1em}(2)}$$

標(biāo)量 $d_M$ 也被稱為推土機(jī)距離（earth mover distance），因?yàn)樗梢越忉尀橹辽僖苿?dòng)多少“泥土”（成本）才能將一個(gè)土堆（分布）變成另一個(gè)土堆（分布）。

4 Sinkhorn算法

4.1 Sinkhorn距離

Sinkhorn距離是對推土機(jī)距離的一種改進(jìn)，在其基礎(chǔ)上引入了熵正則化項(xiàng)：
$$d_M^{\lambda }(\mathbf{r},\mathbf{c})=\underset{P\in U(\mathbf{r},\mathbf{c})}{\min }\sum _{i,j}{P}_{ij}{M}_{ij}?\frac{1}{\lambda }h(P)\text{\hspace{1em}(3)}$$
其中 $h(P)=?\sum P_{ij}\log P_{ij}$ 稱作 $P$ 的信息熵（information entropy），$P$ 分布越均勻，信息熵越大。

熵正則化參數(shù) $\lambda$ 負(fù)責(zé)調(diào)整信息熵的影響程度，$\lambda$ 越大，信息熵的影響越小，最終結(jié)果受成本矩陣的影響更大，即更多地考慮每個(gè)人的喜好；反之，最終結(jié)果則更傾向于均勻分配，每種甜點(diǎn)將平均分配給每個(gè)人。

4.2 算法流程

新增的熵正則化項(xiàng)似乎讓問題更加難以優(yōu)化，但Sinkhorn算法提供了一種簡單且有效的方法應(yīng)對這一問題，Sinkhorn算法認(rèn)為，最優(yōu)分配矩陣 $P_{\lambda }^{?}$ 的元素應(yīng)該具有如下形式：
$$(P_{\lambda }^{?}{)}_{ij}={\alpha }_i{\beta }_j{e}^{?\lambda M_{ij}}\text{\hspace{1em}(4)}$$
其中正是 ${\alpha }_1,...,{\alpha }_n$ 和 ${\beta }_1,...,{\beta }_n$ 使得$P^{?}$ 滿足分配矩陣的三個(gè)約束。如何推導(dǎo)出這一形式可以參考SuperGlue中的最優(yōu)傳輸算法詳解一文。

具體流程如下：

給定： 代價(jià)矩陣 $M$, 分布 $\mathbf{r}$, 分布 $\mathbf{c}$, 熵正則化參數(shù) $\lambda$
初始化： 分配矩陣 $P_{\lambda }=e^{?\lambda M}$
重復(fù)：

1. 縮放行，使得 $P$ 的行和逼近分布 $\mathbf{r}$
2. 縮放列，使得 $P$ 的列和逼近分布 $\mathbf{c}$

直到： 收斂

4.3 代碼實(shí)驗(yàn)

以下是Sinkhorn代碼實(shí)現(xiàn)：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt


r = np.array([3, 3, 3, 4, 2, 2, 2, 1])
c = np.array([4, 2, 6, 4, 4])
M = np.array(
    [[2, 2, 1, 0, 0], 
    [0, -2, -2, -2, -2], 
    [1, 2, 2, 2, -1], 
    [2, 1, 0, 1, -1],
    [0.5, 2, 2, 1, 0], 
    [0, 1, 1, 1, -1], 
    [-2, 2, 2, 1, 1], 
    [2, 1, 2, 1, -1]],
    dtype=float) 
M = -M # 將M變號，從偏好轉(zhuǎn)為代價(jià)

def compute_optimal_transport(M, r, c, lam, eplison=1e-8):
    """
    Computes the optimal transport matrix and Slinkhorn distance using the
    Sinkhorn-Knopp algorithm

    Inputs:
        - M : cost matrix (n x m)
        - r : vector of marginals (n, )
        - c : vector of marginals (m, )
        - lam : strength of the entropic regularization
        - epsilon : convergence parameter

    Outputs:
        - P : optimal transport matrix (n x m)
        - dist : Sinkhorn distance
    """
    n, m = M.shape  # 8, 5
    P = np.exp(-lam * M) # (8, 5)
    P /= P.sum()  # 歸一化
    u = np.zeros(n) # (8, )
    # normalize this matrix
    while np.max(np.abs(u - P.sum(1))) > eplison: # 這里是用行和判斷收斂
        # 對行和列進(jìn)行縮放，使用到了numpy的廣播機(jī)制，不了解廣播機(jī)制的同學(xué)可以去百度一下
        u = P.sum(1) # 行和 (8, )
        P *= (r / u).reshape((-1, 1)) # 縮放行元素，使行和逼近r
        v = P.sum(0) # 列和 (5, )
        P *= (c / v).reshape((1, -1)) # 縮放列元素，使列和逼近c(diǎn)
    return P, np.sum(P * M) # 返回分配矩陣和Sinkhorn距離

我們來看看在不同 $\lambda$ 下，得到的分配矩陣有什么特點(diǎn)：

lam = 0.1

P, d = compute_optimal_transport(M,
        r,
        c, lam=lam)

partition = pd.DataFrame(P, index=np.arange(1, 9), columns=np.arange(1, 6))
ax = partition.plot(kind='bar', stacked=True)
print('Sinkhorn distance: {}'.format(d))
ax.set_ylabel('portions')
ax.set_title('Optimal distribution ($\lambda={}$)'.format(lam))

可以看到每個(gè)人分配得到的甜點(diǎn)基本上都符合初始甜點(diǎn)的分布比例 $\mathbf{c}=[4,2,6,4,4{]}^{?}$ 。

試著調(diào)大 $\lambda$：

可以看到最終的分配向每個(gè)人的偏好靠攏了。文章來源地址http://www.zghlxwxcb.cn/news/detail-434851.html

到了這里，關(guān)于最優(yōu)傳輸問題與Sinkhorn算法的文章就介紹完了。如果您還想了解更多內(nèi)容，請?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！