0. はじめに
この記事では輸送問題の解を「微分」する方法に関して説明します。とくにSinkhornアルゴリズムという手法について詳しく解説し、実際に輸送問題の解の「微分」を計算するPyTorchのサンプルコードを載せました。
解説中に出てくる定理や命題の証明はAppendixにまとめまてあります。AppendixおよびShinkhornアルゴリズムの収束に関する章は数学的な内容なので、興味のない方は飛ばしても問題ありません。
輸送問題と微分
輸送問題はその名の通り、複数の工場から複数の店舗への最適な商品の輸送の仕方を決める問題です。
各工場と店舗間の配送には配送料に応じたコストがかかり、総配送コストを最小化するように各輸送量を決めます。
輸送問題やその他類似の数理最適化に関する入門的な内容は例えば
数理計画法（数理最適化）
が参考になります。
この記事でいう輸送問題の解の「微分」とは、最適な輸送量の各要素を輸送コストの各要素（や各工場の商品数等）で微分した値をさします。
輸送問題の微分の応用
輸送問題の解が微分できると何がうれしいのでしょうか。
一番の利点はNeural Networkの出力を利用して輸送問題を解くタイプのタスクがend to end で学習できることです。
例えば
SuperGlue: Learning Feature Matching with Graph Neural Networks
は、2つの画像間のkey pointのマッチングタスクを扱っています。key pointの局所特徴量とkey point間の輸送コストをNeural Networkで計算し、それらを利用した輸送問題を解くことで最適な画像間の対応を求めています。
Differentiable Ranks and Sorting using Optimal Transport
では、実はソートがある種の輸送問題とみなせることを利用して、
- MNISTの画像のソート問題
- leaet quantile regression
などをend to end で解いています。
上記の論文はいずれも本記事で紹介するSinkhornアルゴリズムを利用しています。
（解の微分を考慮しなくてもShinkhornアルゴリズムはそれ単体で近似解を求めるための有用なアルゴリズムです。）
記号
- $R^n_+$
$\{ x|x\in R_n,x_i \geq 0,i=1,\dots ,n \}$ - $R^n_{+,*}$
$\{ x|x\in R_n,x_i > 0,i=1,\dots ,n \}$ - $1_n$
$(1,...,1) \in R^n$ - $\Sigma_n$
$\{x|x \in R^n_+, \sum x_i=1 \}$ - $U(a,b)$
$a \in R^n_+, b \in R^m_+, U(a,b)= $ $\{ P|P\in R^{n\times m}_+, P1_m=a, P^T1 _n=b \} $ - $\langle P,Q\rangle$
$P\in R^{n\times m},Q\in R^{n\times m},\langle P,Q\rangle=\sum_{i,j}P_{i,j}Q_{i,j}$ - $diag(a)$
$a \in R^n , diag(a):$ 対角成分が$a$である$n\times n$の対角行列
1. 問題の定式化
輸送問題
輸送問題は、$a\in R^n_+,b\in R^m_+, C\in R^{n\times m}_+$が与えられたとき 、$\langle P,C\rangle$を最小化する$P\in U(a,b)$を求めることです。
この記事では
$L_C(a,b)=min_{P\in U(a,b)}\langle P,C\rangle$
と書きます。
$a$が工場の供給量、$b$が店舗の需要量、$C_{i,j}$が工場$i$店舗$j$間の単位量当たりの輸送コスト、$P_{i,j}$が工場$i$店舗$j$間の輸送量、$\langle P,C\rangle$が輸送の総コストですね。最適な輸送を行った場合の総コストが$L_C(a,b)$です。
輸送問題の解の微分？
この記事で輸送問題の解の「微分」と言ったとき、
$P_*\in argmin_{P\in U(a,b)}\langle P,C\rangle$
として$\frac{\partial P_{*k.l}}{\partial C_{i,j}}$や$\frac{\partial P_{*k,l}}{\partial a_i}$のことを指します。
さて、そもそも最適解$P_*$は微分可能なのでしょうか。御想像の通りこのままでは微分可能ではないです。
さらに、何らかの工夫をして「微分っぽい何か」が定義できたとても、特定の例ではその値が常に0になってしまう気もします。例えば次のような輸送問題を考えてみます。
| 工場id | 供給量 | 店舗id | 需要量 | |
|---|---|---|---|---|
| 1 | 1 | a | 1 | |
| 2 | 1 | b | 1 | |
| 3 | 1 | c | 1 |
輸送コスト
| 工場\店舗 | a | b | c |
|---|---|---|---|
| 1 | 2 | 1 | 0 |
| 2 | 1 | 0 | 2 |
| 3 | 0 | 2 | 1 |
最適な輸送量は次のようになるはずです。
| 工場\店舗 | a | b | c |
|---|---|---|---|
| 1 | 0 | 0 | 1 |
| 2 | 0 | 1 | 0 |
| 3 | 1 | 0 | 0 |
まず、工場1と店舗aの輸送コストが仮に0.00001円安くなったとします。しかし、特定の工場店舗間で輸送コストが少し安くなった程度では(最も輸送コストの小さい店舗が変わらない限り)最適な輸送量は全く変わらないと予想できます。したがって仮に「微分っぽい何か」が定義できたとしても常に0になりそうです。
また、これまで輸送問題の「解」と言って来ましたが、解が一つになるとは限りません。
これらの問題をどのように解決すればよいのでしょうか。
正則化項つき輸送問題
上記の問題を解決するために次のような"正則化項"のついた輸送問題を考えます。
つまり、輸送量のエントロピーを
$H(P)=-\sum_{i,j}P_{i,j}(log(P_{i,j})-1)$
として、もとの問題の代わりに
$L_C^{\epsilon}(a,b)=min_{P\in U(a,b)}\langle P,C\rangle - \epsilon H(P)$ ★
を考えます。
$H(P)$が狭義の凸関数なのでこれは唯一の解を持ちますが、この解と、もとの輸送問題の解は次のような関係にあります。
命題 1.
$P_{\epsilon}$を輸送問題$L_C^{\epsilon}(a,b)=min_{P\in U(a,b)}\langle P,C\rangle - \epsilon H(P)$の解とする。このとき次が成り立つ。
$\lim_{\epsilon \to 0} P_{\epsilon} = argmin_{P}\{ -H(P):P \in U(a,b), \langle P,C\rangle=L_C(a,b)\}$
$\lim_{\epsilon \to \infty} P_{\epsilon} = argmin_{P}\{ -H(P):P \in U(a,b) \}$
つまり、正則化項つきの輸送問題★の解は、$\epsilon \to 0$の極限で、もとの問題の最適解のうち、もっともエントロピーが大きいものに収束します。
また、この形の問題はさまざまなよい性質を持ち、次の章で説明するSinkhornアルゴリズムにより、微分可能な形で近似解を求めることができます。
2. Sinkhornアルゴリズム
Shinkhornアルゴリズム
Shinkhornアルゴリズムは次の命題に基づいて★の解を逐次的に求めるアルゴリズムです。
命題 2.
ある$(u,v)\in R^n_+\times R^m_+$が存在して、★の解$P_*$は次のようにかける。
$P_*=diag(u)Kdiag(v)$
ここで$K$は$K_{i,j}=exp(-C_{i,j}/\epsilon)$なる行列。
ここで$u,v$には定数倍の不定性があることに注意してください。実際、$(u,v)\to (zu,v/z),z>0$としても結果は変わりません。
さて、$P_*$も$U(a,b)$の元なので、$u,v$は次のような制約を満たす必要があります。
$diag(u)Kdiag(v)1_m=a$ ①
$diag(v)K^Tdiag(u)1_n=b$ ②
ここで、$diag(u)Kdiag(v)1_m=u*(Kv)$（$*$はベクトルの要素ごとの掛け算、$Kv$は通常の行列とベクトルの掛け算）と書けることを利用して、$u,v$の初期値を適当に決めて
$u^{l+1}=a/(Kv^l)$ ($u$を更新して①を満たすようにする)
$v^{l+1}=b/(K^Tu^{l+1})$ ($v$を更新して②を満たすようにする)
と更新していくのがSinkhornアルゴリズムです。ここで$a/(Kv^l)$などはベクトルの要素ごとの割り算をさします。
はじめの更新式は$P^l=diag(u^l)Kdiag(v^l)$の各行の和を並べると$a$になるように、2番目の更新式は各列の和を並べると$b$になるようにしてます。したがってShinkhornアルゴリズムは、$P^0$の行と列を交互に正規化していくアルゴリズムといえます。
Shinkhornアルゴリズム
init $u=u^0,v=v^0,l=0$, calc K;
while $l$ < MAX_ITER:
$\ \ \ \ u=a/(Kv)$
$\ \ \ \ v=b/(K^Tu)$
$\ \ \ \ l++$
$P=diag(u)Kdiag(v)$
return $u,v,P$
※実際の計算では数値的な安定性のため、logを取ってから計算したり、スケール変換を行ったりします。
Shinkhornアルゴリズムの出力の微分
さて、Shinkhornアルゴリズムは掛け算および$K$の計算に必要な$exp$からなります。いずれも微分可能です！
実際、PyTorchなどの自動微分フレームワークを使用すれば、容易に微分が計算可能です。
例えば、PyTorchを利用してShinkhornアルゴリズムで輸送問題を解く"layer"は次のように実装できます。
import torch
from torch import nn
class OTLayer(nn.Module):
def __init__(self, epsilon):
super(OTLayer,self).__init__()
self.epsilon = epsilon
def forward(self, C, a, b, L):
K = torch.exp(-C/self.epsilon)
u = torch.ones_like(a)
v = torch.ones_like(b)
l = 0
while l < L:
u = a / torch.mv(K,v)
v = b / torch.mv(torch.t(K),u)
l += 1
return u, v, u.view(-1,1)*(K * v.view(1,-1))
実際に簡単な輸送問題を解かせて微分を計算してみます。
C = torch.tensor([
[2., 1., 0.],
[1., 0., 2.],
[0., 2., 1.]
], requires_grad=True)
a = torch.ones(3)/3.
b = torch.ones(3)/3.
ans = torch.Tensor([
[0, 0, 1],
[0, 1, 0],
[1, 0, 0]
])
ans=ans/3.
(輸送問題を解く)
otl = OTLayer(0.1)
u,v,P = otl(C,a,b,10)
print(P)
tensor([[6.8702e-10, 1.5133e-05, 3.3332e-01],
[1.5133e-05, 3.3332e-01, 6.8702e-10],
[3.3332e-01, 6.8702e-10, 1.5133e-05]], grad_fn=<MulBackward0>)
（解の微分）
P[0,0].backward()
print(C.grad)
tensor([[-6.8702e-09, 3.1170e-12, 6.8671e-09],
[ 3.1170e-12, -3.1169e-12, -1.2736e-16],
[ 5.9219e-12, -1.4151e-16, -5.9221e-12]])
はじめに紹介したように、輸送問題を解くOTLayerが作れることで、Neural Networkの出力を利用して輸送問題を解くタイプのタスクがend to end で学習できるようになります。
例えば、
C,a,b = SomeNeuralNetwork(input1, input2)
otl = OTLayer(0.1)
u,v,P = otl(C,a,b,10)
loss = nn.MSELoss()
output = loss(P, ans)
output.backward()
とすると$C,a,b$を計算するSomeNeuralNetworkまで誤差伝播することができます。
Shinkhornアルゴリズムの収束について
（本章とAppendixは数学的な内容になるので飛ばしても問題ありません。）
さて、前章の数値実験ではSinkhornアルゴリズムは正しい解に収束してそうな様子です。実際にそうであることを証明できるのですが、その前に$u,v$などの定数倍の不定性がある対象に対して収束性を議論するために少し準備を行います。
次のようにHilbert metricを定義すると、これは定数倍を同一視した空間(命題 3.の商空間)において距離を定めます。
定義 1.
任意の$u,u'\in R_{+,*}^n$ に対してHilbert metric $d_H$を
$d_H (u,u')=log (max_{i,j} \frac {u_i u_j'} {u_j u_i'})$
と定義する。
命題 3.
$R^n_{+,*}$ 上の同値関係を$u\sim u' \Leftrightarrow \exists z>0, u=zu'$としたとき、$d_H$ は商空間 $ R^n _{+,*} / \sim $ 上の距離を定める。
つまり、任意の $ u,u',u''\in R^n _{+,*} $ に対して次を満たす。
$d_H(u,u')=d_H(u',u)$
$d_H(u,u)\geq 0$
$d_H(u,u')=0 \Leftrightarrow u \sim u'$
$d_H(u,u'')\leq d_H(u,u')+d_H(u',u'')$
このHilbert metricを使って、Sinkhornアルゴリズムが★の最適解に収束することを示す次の定理が記述できます。
定理 1.
★の最適解を$P_*=diag(u_*)Kdiag(v_*)$とする。このとき$u_*,v_*$は$R_{+,*}^n/\sim$の元を一意に定める。
つまり $P_*=diag(u_*)Kdiag(v_*)=diag(u' _*)Kdiag(v' _*)$ とすると
$u _* \sim u' _*, v _* \sim v' _*$
これら$R_{+,*}^n/\sim$の元を同じ記号$u_*,v_*$であらわす。
Shinkhornアルゴリズムで得られるベクトル列${u^l} {v^l}$は$R_{+,*}^n/\sim$上$u_*,v_*$に収束し、特に
$d_H(u^l,u_*)=O(\lambda (K)^{2l})$
$d_H(v^l,v_*)=O(\lambda (K)^{2l})$
となる。ここで、
$\lambda(K)=\frac {\sqrt{\eta(K)} -1}{\sqrt{\eta(K)} +1} < 1$
$\eta (K) = max_{ijkl}\frac{K_{i,k}K_{j,l}}{K_{j,k}K_{i,l}}$
である。
また、
$P^l=diag(u^l)Kdiag(v^l)$
とすると
$| log(P^l)-log(P _* )| _{\infty} \leq d _H (u^l,u _*)+d _H (v^l,v _*)$
を満たす。つまり${P^l}$は$R^{m\times n}$上$P_*$に収束する。
$u^0$の値によって$u^l$は$R^n$上は振動したり異なる値に収束するかもしれませんが、その「向き」および$P^l$は収束し最適解$P_*$に近づくということですね。
3. まとめ
この記事では、
- 輸送問題の微分があるとうれしい事例を紹介しました。
- 正則化項つきの輸送問題を考えるとShinkhornアルゴリズムにより微分可能な形で近似解が得られることを説明しました。
- 簡単なPyTorchのコードでShinkhornアルゴリズムの動作を確認し、解の微分が計算できることを観察しました。
はじめに紹介した事例以外でも、NeuralNetworkの出力を利用して輸送問題を解きたいと事例はいろいろ想像できます。それらのタスクの重要さが確認されれば、Shinkhornアルゴリズムのnn.LayerがPytorchに実装されることになるかもしれません。
4. Appendix
本編で紹介した各命題、定理の証明を行います。
基本的に
を参考にして文献で省略されている部分を埋めました。
命題 1. の証明
命題 1.
$P_{\epsilon}$を輸送問題$L_C^{\epsilon}(a,b)=min_{P\in U(a,b)}\langle P,C\rangle - \epsilon H(P)$の解とする。このとき次が成り立つ。
$\lim_{\epsilon \to 0} P_{\epsilon} = argmin_{P}\{ -H(P):P \in U(a,b), \langle P,C\rangle=L_C(a,b)\}$
$\lim_{\epsilon \to \infty} P_{\epsilon} = argmin_{P}\{ -H(P):P \in U(a,b) \}$
証明
$P^0_*=argmin_{P}\{ -H(P):P \in U(a,b), \langle P,C\rangle=L_C(a,b)\}$とおく。
${ \epsilon_l }$を$\epsilon_l \to 0, \epsilon_l > 0$なる数列とする。また、
$P_l=argmin_{P\in U(a,b)}\langle P,C\rangle - \epsilon_l H(P)$
とする。
$U(a,b)$は有界閉集合なので${ P_l }$は収束する部分列$\{ P _{l _i} \} _i$をもつ。その収束先を$P _*$とおく。
任意の$P \in \{ P:P \in U(a,b), \langle P,C\rangle=L_C(a,b)\}$に対して次が成り立つ。
0 \leq \langle P_{l_i}, C \rangle - \langle P, C \rangle\\
=(\langle P_{l_i}, C \rangle-\epsilon_{l_i} H(P_{l_i}))-(\langle P, C \rangle-\epsilon_l H(P_{l_i})) \\
\leq (\langle P, C \rangle-\epsilon_{l_i} H(P))-(\langle P, C \rangle-\epsilon_{l_i} H(P_{l_i})) \\
=\epsilon_{l_i} (H(P_{l_i})-H(P))
つまり、
0 \leq \langle P_{l_i}, C \rangle - \langle P, C \rangle\leq \epsilon_{l_i} (H(P_{l_i})-H(P)) \\
0\leq H(P_{l_i})-H(P)
ここで$i\to \infty$とすると$\langle P _*, C \rangle = \langle P, C \rangle$および$H(P)\leq H(P _*)$を得る。したがって、$P _*=P^0 _*$。つまり、有界閉集合上の点列$\{P _l \}$の任意の収束する部分列が$P^0 _*$に収束するので、$\{P _l \}$も$P^0
_*$に収束する。これではじめの式が示せた。
2番目の式も同様の方法で示せる。
$P^{\infty} _*=argmin _{P}\{ -H(P) \}$とおく。
$\{ z _l \}$を$z _l \to \infty, z _l > 0$なる数列とする。また、
$P' _l=argmin _{P\in U(a,b)}\langle P,C\rangle - z _l H(P)$
とする。
$U(a,b)$は有界閉集合なので$\{ P' _l \}$は収束する部分列$\{ P' _{l _i} \} _i$をもつ。その収束先を$P' _*$とおく。ここで
0\leq -(H(P' _{l _i}))-(-H(P _*^{\infty})) \\
=(\frac{1}{z _{l _i}}\langle P' _{l _i}, C\rangle - H(P' _{l _i})) - (\frac{1}{z _{l _i}}\langle P'_{l _i}, C\rangle - H(P _*^{\infty})) \\
\leq (\frac{1}{z _{l _i}}\langle P _*^{\infty}, C\rangle - H(P _*^{\infty})) - (\frac{1}{z _{l _i}}\langle P' _{l _i}, C\rangle - H(P _*^{\infty})) \\
= \frac{1}{z _{l _i}}(\langle P _*^{\infty}, C\rangle - \langle P' _{l _i}, C\rangle)
が成り立つ。ここで$i\to \infty$とすると$H(P' _*)=H(P _*^{\infty})$つまり$P' _* = P _*^{\infty}$。はじめの式の場合と同様の議論で$\{ P' _l \}$も$P _*^{\infty}$に収束する。■
命題 2. の証明
命題 2.
ある$(u,v)\in R^n_+\times R^m_+$が存在して、★の解$P_*$は次のようにかける。
$P_*=diag(u)Kdiag(v)$
ここで$K$は$K_{i,j}=exp(-C_{i,j}/\epsilon)$なる行列。
証明
ラグランジュの未定乗数法を用いる。$P\in U(a,b)$という制約は$P1_m=a, P^T 1_n=b$という、それぞれn個m個の条件式を与えるが、それぞれの条件式に対する未定乗数を並べたものを$f\in R^n, g\in R^m$とする。
$$
L(P,f,g)=\langle P, C \rangle - \epsilon H(P) - \langle f, P1_m-a \rangle - \langle g, P^T 1_n-b \rangle
$$
の$P$による微分$=0$とおくと
$$
\frac{\partial L}{\partial P_{i,j}} = C_{i,j}+\epsilon log(P_{i,j})-f_i -g_j = 0
$$
より$P_{i,j}=exp(f_i / \epsilon) exp(-C_{i,j}/ \epsilon)exp(g_j / \epsilon)$。$u_i=exp(f_i / \epsilon), v_j=exp(g_j / \epsilon)$とすれば命題の式を得る。■
定理 1. の証明
定理 1.
★の最適解を$P_*=diag(u_*)Kdiag(v_*)$とする。このとき$u_*,v_*$は$R_{+,*}^n/\sim$の元を一意に定める。
つまり $P_*=diag(u_*)Kdiag(v_*)=diag(u' _*)Kdiag(v' _*)$ とすると
$u _* \sim u' _*, v _* \sim v' _*$
これら$R_{+,*}^n/\sim$の元を同じ記号$u_*,v_*$であらわす。
Shinkhornアルゴリズムで得られるベクトル列${u^l} {v^l}$は$R_{+,*}^n/\sim$上$u_*,v_*$に収束し、特に
$d_H(u^l,u_*)=O(\lambda (K)^{2l})$
$d_H(v^l,v_*)=O(\lambda (K)^{2l})$
となる。ここで、
$\lambda(K)=\frac {\sqrt{\eta(K)} -1}{\sqrt{\eta(K)} +1} < 1$
$\eta (K) = max_{ijkl}\frac{K_{i,k}K_{j,l}}{K_{j,k}K_{i,l}}$
である。
また、
$P^l=diag(u^l)Kdiag(v^l)$
とすると
$| log(P^l)-log(P _* )| _{\infty} \leq d _H (u^l,u _*)+d _H (v^l,v _*)$
を満たす。つまり${P^l}$は$R^{m\times n}$上$P_*$に収束する。
証明
★の最適解を$P _*=diag(u _*)Kdiag(v _*)=diag(u' _*)Kdiag(v' _*)$とする。任意の$i,j$について$u _{* i}K _{i,j}v _{*j} =u' _{*i} K _{i,j} v' _{*\ j}$、つまり$u _{*i} v _{*j}=u' _{*i}v' _{*j}$。したがって$j$もしくは$i$を固定するれば$u _* \sim u' _*, v _* \sim v' _*$を得る。
$u^l,v^l$が$R _{+,*}^n/\sim$上$u _*,v _*$に収束することの証明には次の定理を利用する。
定理 2.
任意の$K \in R _{+,*}^{n \times m}, v,v'\in R _{+,*}^m$に対して
$$
d _H(Kv,Kv')\leq \lambda(K) d _H(v,v')
$$
(定理 2. の証明は省略させてください)
ここで、任意の$v,v'\in R _{+,*}^n/\sim$に対して
d _H(v,v')=d_H(v/v', 1 _n)=d _H(1 _n /v, 1 _n /v')
となること、また、$P _* \in U(a,b)$なのでShinkhornアルゴリズムの章と同様の議論で
u _*=\frac{a}{Kv _*},v_*=\frac{b}{Ku _*}
が成り立つことに注意する。すると、
d_H(u^{l+1},u_*)=d_H(\frac{a}{Kv^l}, \frac{a}{Kv_*}) \\
= d_H(Kv^l, Kv_*) \leq \lambda (K)d_H(v^l, v_*) \\
=\lambda (K) d_H(\frac{b}{Ku^l}, \frac{b}{Ku_*}) \\
= \lambda (K) d_H(Ku^l, Ku_*) \leq \lambda (K)^2d_H(u^l, u_*) \\
\leq \dots \leq \lambda (K)^{2l+1}d_H(u^0, u_*) \\
$\lambda (K) < 1$なので$\{ u^l \}$は$u _*$に収束し、$u _*$との距離は$O(\lambda(K)^{2l})$のオーダーで小さくなる。
$\{v^l\}$に関しても同様。
次に${P^l}$が$R^{n\times m}$上$P_*$に収束することを示す。
まず、
diag(u^l / u_*)P_* diag(v^l / v_*) = P^l
だが、つねに
log(u^l_i/u_{*,i}) \leq d_H(u^l/u_*, 1_m)= d_H(u^l, u_*)\leq d_H(u^l, u_*) + d_H(v^l, v_*), i=1,\dots n.
となることに注意する。
また、必要であれば$u_*$を定数倍して、
1 \leq u^l_i/u_{*,i}, i=1,\dots n.
が成り立つようにする。
つまり
1 \leq u^l_i/u_{*,i}\leq exp(d_H(u^l, u_*) + d_H(v^l, v_*)), \ i=1,\dots n. ●
となる。この不等式をまとめて
1_n \leq diag(u^l/u_*)1_n \leq exp(d_H(u^l, u_*) + d_H(v^l, v_*))1_n
と書くことにする。逆数をとって
exp(-(d_H(u^l, u_*) + d_H(v^l, v_*)))1_n \leq diag(u_*/u^l)1_n \leq 1_n
つぎに各辺に$P^{l,T}$をかける。$P^{l,T}1_n=b$なので
exp(-(d_H(u^l, u_*) + d_H(v^l, v_*)))b \leq P^{l,T} diag(u_*/u^l)1_n \leq b
ここで$diag(u^l / u _* )P _* diag(v^l / v _*) = P^l$より両辺の転置をとって右側から$diag(u _* /u^l)$をかけると
diag(v^l / v_*)P_*^T = P^{l,T}diag(u_*/u^l)
さらに$P_*^T1_n=b$なので
P^{l,T} diag(u_*/u^l)1_n = diag(v^l/v_*)b
したがって上の不等式は
exp(-(d_H(u^l, u_*) + d_H(v^l, v_*)))b \leq diag(v^l/v_*)b \leq b
つまり
exp(-(d_H(u^l, u_*) + d_H(v^l, v_*))) \leq v^l_j/v_{*,j} \leq 1,\ j=1,\dots m.
●と各辺をかけ合わせると
j=1,\dots m, i=1,\dots n \\
exp(-(d_H(u^l, u_*) + d_H(v^l, v_*))) \leq v^l_ju^l_i/v_{*,j}u_{*,i} \leq exp(d_H(u^l, u_*) + d_H(v^l, v_*)),
を得る。
P^l_{i,j}/P_{*,i,j}=v^l_ju^l_i/v_{*,j}u_{*,i}
なのでlog+maxをとると
\|log(P^l)-log(P_*)\|_{\infty} \leq d_H(u^l,u_*)+d_H(v^l,v_*)
となる。■