Upgrade to Pro — share decks privately, control downloads, hide ads and more …

最適輸送の解き方

3.8k

 最適輸送の解き方

最適輸送問題(Wasserstein 距離)を解く方法についてのさまざまなアプローチ・アルゴリズムを紹介します。
線形計画を使った定式化の基礎からはじめて、以下の五つのアルゴリズムを紹介します。
1. ネットワークシンプレックス法
2. ハンガリアン法
3. Sinkhorn アルゴリズム
4. ニューラルネットワークによる推定
5. スライス法
このスライドは第三回 0x-seminar https://sites.google.com/view/uda-0x-seminar/home/0x03 で使用したものです。自己完結するよう心がけたのでセミナーに参加していない人にも役立つスライドになっています。

『最適輸送の理論とアルゴリズム』好評発売中! https://www.amazon.co.jp/dp/4065305144

佐藤竜馬 (Ryoma Sato)

January 13, 2023
Tweet

Transcript

  1. 2 / 270 KYOTO UNIVERSITY 最適輸送アルゴリズムの様々なアプローチを紹介する  最適輸送のアルゴリズムには様々なアプローチがあり それぞれ異なる長所・短所を持っている 

    今日のゴール: 解きたい問題に合わせて適切なアルゴリズムを選べるようになる  最適輸送の主・双対問題は豊かな構造を持っている  それらを操れるようになると、アルゴリズムを選ぶだけでなく自分で 設計できるようになる  今日の発展的ゴール: 最適輸送の主・双対問題を自由自在に使いこなせるようになる 意欲のある方向け
  2. 3 / 270 KYOTO UNIVERSITY もくじ  p.7 今日の講演の構成 

    p.17 線形計画による定式化  p.45 ネットワークシンプレックス法  p.128 双対問題について  p.149 ハンガリアン法  p.198 Sinkhorn アルゴリズム  p.222 ニューラルネットワークによる推定  p.245 スライス法  p.262 まとめ  p.267 参考文献
  3. 5 / 270 KYOTO UNIVERSITY 自己紹介: 京大で最適輸送や GNN を研究してます 

    名前: 佐藤竜馬  京大博士一年(鹿島・山田研究室) 研究分野: グラフニューラルネットワーク (GNN) * Ryoma Sato, Makoto Yamada, Hisashi Kashima. Approximation Ratios of Graph Neural Networks for Combinatorial Problems. NeurIPS 2019. * Ryoma Sato, Makoto Yamada, Hisashi Kashima. Random Features Strengthen Graph Neural Networks. SDM 2021. 最適輸送  今日のトピック * Ryoma Sato, Makoto Yamada, Hisashi Kashima. Fast Unbalanced Optimal Transport on a Tree. NeurIPS 2020. * Yuki Takezawa, Ryoma Sato, Makoto Yamada. Supervised Tree-Wasserstein Distance. ICML 2021.  講演者近影 「人工知能」七月号に GNN の 日本語解説記事書きました
  4. 6 / 270 KYOTO UNIVERSITY 自己紹介: 競技プログラミングしてました  高校生〜学部生の頃は競技プログラミングをしていました 

    主な戦績(自慢): 日本情報オリンピック 2014 金賞 国際情報オリンピック 2014 日本代表、銅メダル ACM ICPC 2018 アジア地区予選ジャカルタ大会 2位 ACM ICPC 2019 世界大会出場  アルゴリズムの素養は主にここからきています  機械学習 + アルゴリズムの交わりを得意にしています
  5. 8 / 270 KYOTO UNIVERSITY 今日は最適輸送の色んな解き方をみていく  以下の五つのアルゴリズムの長所・短所をみていく  その1:

    線形計画の定式化とシンプレックス法  その2: ハンガリアン法  その3: Sinkhorn アルゴリズム  その4: ニューラルネットワークによる推定  その5: スライス法
  6. 9 / 270 KYOTO UNIVERSITY 線形計画の定式化は全ての基本  その1: 線形計画の定式化はその後の解き方の基本になる 

    ネットワークシンプレックス法はこれを主問題で解くアルゴリズム  長所 : ソルバーがよく整理されている CPU 上だと早い 厳密に解ける  短所 : 実行時間を見積もるのが難しい
  7. 10 / 270 KYOTO UNIVERSITY ハンガリアン法は一対一の割り当て問題によく使われる  その2: ハンガリアン法は線形計画の双対を解く 

    ユースケースは限られるが、一対一の割当問題を解くためには よく使われている  長所 : CPU 上だと早い シンプルに実装できる  短所 : 一対一マッチングにしか使えない (対応するように拡張はできる) 応募者 仕事
  8. 11 / 270 KYOTO UNIVERSITY Sinkhorn は微分可能で GPU フレンドリー 

    その3: Sinkhorn アルゴリズムは目的関数に正則化を加え、 行列演算で最適輸送を解くアルゴリズム  Sinkhorn が [Cuturi+ 2013] で紹介されたことで 機械学習界隈で最適輸送が使われるようになるキッカケになった  長所 :コストに関して微分可能 GPU 上だと早い 複数の最適輸送を並列に解ける 最適解が唯一に定まる  短所 : 厳密な最適輸送の解ではない オーバーフローや解の丸めに気をつける必要あり Marco Cuturi. Sinkhorn Distances: Lightspeed Computation of Optimal Transport. NeurIPS 2013. Gabriel Peyré, Marco Cuturi. Computational Optimal Transport. 2019.
  9. 12 / 270 KYOTO UNIVERSITY ニューラルネットワークによる推定は GAN の基本技術  その4:

    ニューラルネットワークによる推定は双対問題の特殊形を ニューラルネットワークの最適化に帰着する  Wasserstein GAN という名前で GAN(生成モデル)の 生成分布とデータ分布の距離をはかるのに使われている  長所 : コストについて微分可能 分布からのサンプルをミニバッチで繰り返し解くのに 向いている  短所 : 単一の入力には向いていない コストが距離で定義されていないと 使えない Tero Karras, Timo Aila, Samuli Laine, Jaakko Lehtinen. Progressive Growing of GANs for Improved Quality, Stability, and Variation. ICLR 2018.  GAN による顔写真の自動生成
  10. 13 / 270 KYOTO UNIVERSITY スライス法は簡単かつ高速に計算ができる  その5: スライス法は一次元上の最適輸送が簡単に解けることを 利用する

     非常に簡単であることから、とりあえず最適輸送っぽいものを得たい ときには有用  長所 : 実装が簡単 非常に高速  短所 : 元の最適輸送問題の良い近似に なっているとは限らない https://www.programmersought.com/article/67174999352/
  11. 14 / 270 KYOTO UNIVERSITY 今日は最適輸送の色んな解き方をみていく  その1: 線形計画の定式化  

     その2: ハンガリアン法  その3: Sinkhorn アルゴリズム  その4: ニューラルネットワークによる推定  その5: スライス法  上の話題とは独立している それぞれ違う アプローチで 双対問題を解く  一番大事! 時間をかけて丁寧に追う
  12. 15 / 270 KYOTO UNIVERSITY 今日の講演の聴き方について  スライドの中では証明がたびたび登場します  厳密性よりは直感を提供することを優先して説明します

     重要な洞察は証明後にあらためてフォローするなどして、 証明については理解せずに聴き飛ばしても追えるような構成には しているので安心してください  各章はほぼ独立していますが、たとえばニューラルネットワークの推定 を読む前に以前の章で双対問題の直感を養っておくとよりスムーズに 理解できるはずです  このスライドはウェブ上で公開しているので、講演が終わったあとに 証明部分を読み直して理解を深めるなどして活用してください
  13. 16 / 270 KYOTO UNIVERSITY 凡例  重要な用語は赤色で強調している  重要なメッセージは青色で強調している

    フキダシは補足情報 地文で理解できないときの 助けを目的としている
  14. 19 / 270 KYOTO UNIVERSITY 最適輸送は最小の輸送コストを求める問題  最適輸送問題(ざっくり): 工場会社が倉庫を n

    個、工場を m 個所有している 倉庫 i には小麦粉が a i > 0 グラムあり工場 j には小麦粉が b j > 0 グラム必要である。 倉庫 i から工場 j に 1 グラム輸送するのにコスト C ij かかるとき、 要件を満たすのに必要な最小コストは?
  15. 20 / 270 KYOTO UNIVERSITY 例題  例: a 1

    = 2 a 2 = 1 a 3 = 1 b 1 = 1 b 2 = 2 b 3 = 1 C 1 2 3 1 1 2 2 2 2 1 2 3 2 2 1
  16. 21 / 270 KYOTO UNIVERSITY 例題答え  例 答え: 1

    + 2 + 1 + 1 = 5 a 1 = 2 a 2 = 1 a 3 = 1 b 1 = 1 b 2 = 2 b 3 = 1 C 1 2 3 1 1 2 2 2 2 1 2 3 2 2 1 1*コスト1 1*コスト2 1*コスト1 1*コスト1
  17. 22 / 270 KYOTO UNIVERSITY 最適輸送のもう少しフォーマルな定義  最適輸送問題: 入力: ベクトル

    a ∈ Rn, ベクトル b ∈ Rm, 行列 C ∈ Rn✕m 出力: {a 1 , ..., a n } を {b 1 , ..., b m } に輸送する最小コスト  ただし a 1 + ... + a n = b 1 + ... + b m = 1 を仮定する つまり「小麦粉」は不足することも余ることもない 合計が 1 なのは割合だけ考えていると思えばよい a i から b j に P ij だけ輸送するのには C ij ✕ P ij コストかかる 輸送コストは輸送量に線形増加 R, は実数集合 平文では黒板太字を使うのが 面倒なので R を使用します
  18. 23 / 270 KYOTO UNIVERSITY 最適輸送の線形計画による定式化  最適輸送問題: 入力: ベクトル

    a ∈ Rn, ベクトル b ∈ Rm, 行列 C ∈ Rn✕m 最適化変数 P ∈ Rn✕m : P ij は i から j に輸送する量  総コスト  輸送量は非負  余りなし  不足なし あわせて質量制約とよぶ (mass constraint) C が距離行列 D と正数 p ≧ 1 を用いて C ij = D ij p を表せると き、最適値の 1/p 乗の値を Wasserstein 距離という。 最適輸送と Wasserstein 距離 を区別せずに用いることもある
  19. 25 / 270 KYOTO UNIVERSITY これからは問題の性質を知り高速化することを目指す  線形計画に落とし込めたので、線形計画ソルバ (CPLEX など)

    に投げれば解ける  めでたしめでたし  ... とは限らない!  第一に、汎用ソルバは専用ソルバより一般に遅い  第二に、最適輸送の答えが小さくなるような配置を見つけたい、 のように別の最適化の中に最適輸送を組み込むときには 最適輸送の最適化の性質を知っておく必要がある  ただし、小さいインスタンスについて最適値を知りたいだけなら本当に ここでおしまいです。お疲れさまでした。
  20. 26 / 270 KYOTO UNIVERSITY 双対は最適値が同じで相補性が成り立つ  線形計画の双対: 最小化の線形計画には対となる最大化の線形計画問題がある もとの線形計画を主問題、対となる線形計画を双対問題とよぶ

    ざっくり言うと、主問題の制約一つにつき双対問題に変数が一つ、 主問題の変数一つにつき双対問題に制約が一つできる  重要な性質(いまから詳しく見ていく): 主問題と双対問題の最適値は一致する(強双対性) 主問題の最適解の変数が正のとき双対問題の対応する制約は 等号成立する。逆にこれが成り立つとき最適解となる(相補性)  最適輸送では主問題と双対問題両方を使いこなすことが重要 この講演の目標は主・双対問題を自由自在に操れるようになること
  21. 27 / 270 KYOTO UNIVERSITY 最適輸送の双対問題  最適輸送の双対問題は右のように表せる 導出は省略。普通の線形計画と同様にできるので試してみるとよい。 主問題

    双対問題 f i : a i の余りなし制約に対応する変数 g j : b j の不足なし制約に対応する変数 f i + g j ≦ C ij : 主問題の変数 P ij に対応する制約
  22. 29 / 270 KYOTO UNIVERSITY 物流業者登場。双対問題に対応する  トラックで小麦粉を輸送するのをやめて、外部の物流業者を利用する  物流業者は倉庫

    i の地域ではグラム f i 円払うと引き取りを行い 工場 j の地域ではグラム g j 円で販売を行う  倉庫の小麦粉を全て引き取ってもらい、工場で必要な分を 買い取ると実質輸送が完了したことになる
  23. 32 / 270 KYOTO UNIVERSITY 物流業者が双対問題を解く  双対問題の主人公、双対問題を「解く」人は物流業者  物流業者は倉庫

    i で a i の小麦を単価 f i で回収するので a i f i 円 工場 j で b j の小麦を単価 g j で売るので b j g j 円売り上げる  物流業者の売上は これを最大化したがっている 双対問題の目的関数!
  24. 33 / 270 KYOTO UNIVERSITY 明らかな工場会社の得を維持しつつ利益を最大化する  しかし、あまりに引取価格や販売価格を高くしすぎると工場会社に 利用してもらえない 

    f i + g j ≦ C ij であれば、工場会社は倉庫 i から工場 j にトラックで 移すより引取 + 購入を行った方が得する  この条件があらゆる (i, j) について成り立っていれば、工場会社は 安心して物流業者を利用できる 少なくとも自社でトラックを動かすより物流業者に頼ったほうがコストが 低いことが保証される  物流業者は工場会社が安心できる範囲で売上を最大化すること をめざす → 双対問題 双対問題の制約!
  25. 35 / 270 KYOTO UNIVERSITY 最適輸送の双対問題は変数が少ないが制約が多い  主問題は各場所対の輸送量を定める問題 双対問題は各場所の販売・回収コストを定める問題 

    主問題は変数が nm 個、制約が n + m 個 双対問題は変数が n + m 個、制約が nm 個  どちらが解きやすいかはケースバイケース 大 小 大 小
  26. 36 / 270 KYOTO UNIVERSITY 双対問題の最適解には自由度あり  双対問題の最適解には定数ぶんだけ自由度がある  が最適解なら定数

    について も最適解  理由: なら また、  販売費を安くするぶん回収費を高くして帳尻を合わせられる  平均ゼロや特定変数の値が 0 になる解を正規の解とすることがある 実行可能 目的関数の値は変わらない
  27. 37 / 270 KYOTO UNIVERSITY 弱双対定理: 工場会社は物流業者に任せると損しない  弱双対定理: 任意の主問題の実行可能解

    P と双対問題の実行可能解 (f, g) について、  つまり、双対問題の目的関数の値は主問題の目的関数の値より 常に小さい  直感的には先程の議論から明らか (f, g) をどのような輸送プランに対しても明らかに損しないように選ぶ のが実行可能条件
  28. 38 / 270 KYOTO UNIVERSITY 弱双対定理の証明  弱双対定理の証明 P の実行可能性より

    P ij ≧ 0 (f, g) の実行可能性より C ij ≧ f i + g j P の実行可能性より ∑ j P ij = a i , ∑ i P ij = b j
  29. 39 / 270 KYOTO UNIVERSITY 強双対定理: 物流業者が頑張るとプラマイゼロまでいける  強双対定理: を満たす主問題の

    実行可能解 P と双対問題の実行可能解 (f, g) が存在する  つまり、双対問題と主問題の最適値は一致する  小麦粉の例でいうと、物流業者が値段設定を頑張れば、 工場会社が自社で輸送する最適コストまでは売上を伸ばせる (それ以上は弱双対定理より不可能)  証明は少し煩雑なので省略 線形計画や最適化の教科書を参照 主問題 最小化 双対問題 最大化 ちょうど一致
  30. 40 / 270 KYOTO UNIVERSITY 相補性: 主問題で正の箇所は等号成立  定義 相補性: 主問題の実行可能解

    P について、双対問題の解 (f, g) は以下の 条件を満たすとき相補的であるという P ij > 0 ならば f i + g j = C ij  つまり、輸送が発生する倉庫と工場の間では不等式制約が 等号で成立するとき相補的  等号が成立するということは、その間の輸送を流通業者に任せた としても(もちろん損することはないが)得せずコストは変わらない  f i + g j = C ij が成り立つ組 (i, j) をタイトと呼ぶ f i + g j < C ij が成り立つ組 (i, j) をルーズと呼ぶ 不等式制約がギリギリ 縁まで来ているイメージ
  31. 41 / 270 KYOTO UNIVERSITY 相補的ならば両方最適解  命題 相補性を使った最適性の十分条件: 主問題の実行可能解 P

    について双対問題の実行可能解 (f, g) が相補的であるならばこれらは両方それぞれの問題の最適解  証明(弱双対と似ている): 相補性より正の項については C ij = f i + g j それ以外の項は何であっても消えるので無問題 弱双対定理より、 主問題より目的関数は大きく なり得ないのでこれが最適
  32. 42 / 270 KYOTO UNIVERSITY 相補的でないならどちらかは最適解でない  命題 相補性を使った最適性の必要条件: 主問題の実行可能解 P

    について双対問題の実行可能解 (f, g) が相補的でないならばこれらのどちらかは最適解でない  証明(弱双対と似ている): 非相補性より P ij > 0, C ij > f i + g j なる項があるので等号不成立 強双対定理より、 主問題か双対問題のどちらかは さらに改善できる
  33. 43 / 270 KYOTO UNIVERSITY 最適化をすることは相補的な組を見つけることと同じ  つまり、実行可能解 P と

    (f, g) の両方が最適解である 必要十分条件は相補的であること  主問題や双対問題を最適化する代わりに、相補的な組を見つける ことをめざしてもよい  これで線形計画としての基本的な性質は終了 ここからはこの性質を使ってどのように解くかをみていく
  34. 44 / 270 KYOTO UNIVERSITY この章のゴールの確認  最適輸送の主問題と双対問題の線形計画による定式化を知る 主問題は工場会社が全ての小麦を輸送する最小コスト 双対問題は物流業者が全ての小麦を回収・販売する最大売上

     主問題と双対問題の関係を知る 双対問題の値は主問題の値より小さい(弱双対性) 双対問題の最適値は主問題の最適値と同じ(強双対性) 主・双対実行可能解が相補的であるときかつそのときのみ最適
  35. 46 / 270 KYOTO UNIVERSITY 線形計画の定式化は全ての基本  その1: 線形計画の定式化はその後の解き方の基本になる 

    ネットワークシンプレックス法はこれを主問題で解くアルゴリズム  長所 : ソルバーがよく整理されている CPU 上だと早い 厳密に解ける  短所 : 実行時間を見積もるのが難しい
  36. 47 / 270 KYOTO UNIVERSITY シンプレックス法は一般の線形計画用の効率的手法  シンプレックス法: 線形計画の実行可能領域が凸多面体であること・最適値が 頂点にあることを利用し、多面体の頂点を解が良くなる方向に

    走査していくアルゴリズム  最悪ケースでは指数時間かかるが、 多くのケースでは非常に高速である  シンプレックス法は一般の線形計画用のアルゴリズム  ここではシンプレックス法を最適輸送専用に高速化・単純化した ネットワークシンプレックス法を紹介する  どのように頂点を辿るかが効率の上で重要
  37. 49 / 270 KYOTO UNIVERSITY 別の解の中点になっていない解を頂点という P Q R P

    Q R P Q R 頂点以外だとそのような Q, R を取れるが 頂点だと取ろうとすると はみ出てしまう  まず、最適輸送問題の実行可能領域の頂点はどのような解に対応 しているかを考える  定義 頂点: 実行可能解 P が頂点であるとは、別の実行可能解 Q, R を用いて P = (Q + R) / 2 と表せないことである
  38. 51 / 270 KYOTO UNIVERSITY 使用グラフ: 輸送が発生する倉庫と工場に辺がある P 1 2

    3 1 1 1 0 2 0 1 0 3 0 0 2 1 2 3 1 2 3 使用グラフは この資料独自の用語 小麦問題の類推から 左側のノードを倉庫 右側のノードを工場とよぶ  定義 使用グラフ: 左側に n 個、右側に m 個ノードを置き、左側の i と右側の j に 辺を張った(二部)グラフ G を考える。 主問題の解 P について P ij > 0 であれば (i, j) が G の辺として 含まれるとき、G は P に付随する使用グラフという  実際に輸送が発生する倉庫と工場間に辺が張られているようなグラフ
  39. 53 / 270 KYOTO UNIVERSITY 頂点解は使用グラフとして木をもつ  命題 頂点の必要条件: 主問題の実行可能解 P

    が頂点ならば、木を使用グラフとしてもつ  証明(ざっくり): 対偶、木を使用グラフとしてもたなければ頂点でないを示す このとき P ij > 0 となる (i, j) だけから閉路 H を構成できる H を適当な倉庫から辿っていくと、対応する頂点は倉庫、工場、 倉庫、工場、倉庫、... というように倉庫と工場が交互に現れる ちょうど一周した時点でストップすることにする このとき通過した辺に対応する P の値を順番に P 1 , ..., P L とおく つづく 倉庫どうし、 工場どうしは 結ばれることが ないため L は閉路の長さ P 2 P 1 P 3 P 4 P 3 P 5 P 6
  40. 54 / 270 KYOTO UNIVERSITY 頂点解は使用グラフとして木をもつ:証明つづき  証明(続き): 偶数番目の要素を ε

    > 0 大きくし、奇数番目の要素を ε 小さく して P から新しい解 P’ を作る(閉路の要素以外は変えない) P i は正辺なので値は正 ε を十分小さくとると解は正のまま H に含まれる各工場について倉庫 → 工場側の辺が ε 小さくなり、 工場 → 倉庫側の辺が ε 大きくなるので不足なし制約を満たす 倉庫についても同様に余りなし制約を満たすので、P’ は実行可能 ε’ = -ε で同じことすると 実行可能解 P’’ が得られ P = (P’ + P’’) / 2 なので P は頂点でない Gabriel Peyré, Marco Cuturi. Computational Optimal Transport. 2019.
  41. 55 / 270 KYOTO UNIVERSITY 木が付随する解はただ一つ  補題 木が付随する解の唯一性: 木を任意の取ったときこの質量制約を満たす主問題の解のうち この木が付随するものはただ一つ

     証明(ざっくり): 木が定まれば主問題の解の値が次々に定まることを示す グラフが木なので次数が 1 のノード、すなわち一本しか使用グラフの 辺が繋がっていないノードが存在する このノードを i とよび、簡単のため倉庫と仮定する ノード i に接続する辺 (i, j) に相当する P ij の値は余りなし制約 より P ij = a i と一意に定まる つづく i から出る他の辺は 使用グラフに含まれないので使えない j に全て送るしかない P が負かもしれないので 実行可能とは限らない
  42. 56 / 270 KYOTO UNIVERSITY 木が付随する解はただ一つ: 証明つづき  証明(続き): P

    ij が確定したことにより、j 側の要求が P ij だけ満たされたので b j から P ij を引く P ij が確定したのでグラフから頂点 i と辺 (i, j) を取り除く こうしてもグラフは木のまま あとは再帰的に上記のステップを繰り返すと全ての辺の値が一意に 定まる  詳しくは次のスライドからはじまる例を参照
  43. 58 / 270 KYOTO UNIVERSITY 木の使用グラフが定まると解は一意に定まる: 例  例: 0.1

    0.2 0.5 0.1 0.1 0.6 0.2 0.2 次数 1 なので、相手から必ず 0.1 輸送されないといけない
  44. 59 / 270 KYOTO UNIVERSITY 木の使用グラフが定まると解は一意に定まる: 例  例: 0.0

    0.2 0.5 0.1 0.1 0.5 0.2 0.2 0.1 輸送したので残りが減る 0.1 この辺の輸送量が 0.1 に決まる
  45. 62 / 270 KYOTO UNIVERSITY 木の使用グラフが定まると解は一意に定まる: 例  例: 0.0

    0.0 0.5 0.1 0.1 0.3 0.2 0.2 0.1 0.2 最初次数 1 でなかったが 輸送コストが確定していって 次数が 1 になった 残った小麦は残りの辺を使って 全て輸送するしかない
  46. 67 / 270 KYOTO UNIVERSITY 木の使用グラフが定まると解は一意に定まる: 例  例: 0.0

    00 0.0 0.1 0.1 0.0 0.0 0.2 0.1 0.2 0.3 0.2 0.0 0.0 輸送されるというパターンもある
  47. 71 / 270 KYOTO UNIVERSITY 木の使用グラフが定まると解は一意に定まる: 例  例: 0.0

    0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.2 0.3 0.2 0.0 0.1 0.1 ∑ a i = ∑ b j の仮定より 必ず最後は辻褄があう
  48. 81 / 270 KYOTO UNIVERSITY 木の使用グラフが定まると解は一意に定まる: 例  別の例:無理やり続行 0.0

    0.0 0.0 0.3 0.1 0.0 0.2 0.2 0.1 0.2 0.5 -0.2 負の量を送られると逆に増える
  49. 83 / 270 KYOTO UNIVERSITY 木の使用グラフに付随する解が実行不能のときもある  輸送量を定めるうちに在庫量・要求量が負になることもある  そのままプロセスを続行させることはできるが、輸送プランにも負値が

    現れ実行可能ではない場合がある  つまり、木 T を固定したとき、この木を使用グラフとして持つ解は ちょうど 1 つ ただし、実行可能でない場合はあるので実行可能解は 0 or 1 個
  50. 84 / 270 KYOTO UNIVERSITY 木を使用グラフとして持つなら頂点  命題 頂点の十分条件: 主実行可能解 P

    が木を使用グラフとしてもつなら P は頂点  証明(ざっくり): 実行可能解 Q, R を用いて P = (Q + R) / 2 と表せるとする Q と R の要素は非負なので P ij = 0 である要素は Q, R でもゼロ (仮に Q ij > 0 なら R ij ≧ 0 より Q ij + R ij > 0 となる) よって Q, R にも P と同じ木が使用グラフとして付随し、一意性より P = Q = R
  51. 85 / 270 KYOTO UNIVERSITY 頂点の条件より探索空間を劇的に絞れる  証明はややこしいが、  実行可能解

    P が頂点であることと使用グラフとして木をもつ ことは同値である  木にはちょうど一つの解に付随 を知っていればよい  線形計画の最適解は頂点にあるので、これにより解の候補を絞れる
  52. 86 / 270 KYOTO UNIVERSITY 頂点だけ探索すれば良い → 頂点を全て探索してみる  愚直なアルゴリズム:

    倉庫と工場を結ぶ木を全て考え 各木を使用グラフとしてもつ解(一意に定まる)を計算し、 コストが一番安いものを出力する  木のパターンは非常に多いので効率は悪い しかし、解集合が無限にあったところを、有限まで絞り込めているので 劇的な改善
  53. 87 / 270 KYOTO UNIVERSITY 使用グラフと相補性は相性がよい  使用グラフの定義は相補性と相性が良さそう  使用グラフ

    G の全ての辺 (i, j) で f i + g j = C ij となるとき (f, g) は G について相補的であるという。  主問題の実行可能解 P と双対問題の実行可能解 (f, g) を考える P に付随する任意の使用グラフ G をとる (G, f, g) が相補的である → (P, f, g) が相補的である 理由: P ij > 0 な辺 (i, j) は使用グラフ G に含まれるのでタイト  P ij > 0 でない辺もタイトなので、使用グラフについての相補性は P についての相補性より強い  しかし、特殊な例ではこれが十分条件になる
  54. 88 / 270 KYOTO UNIVERSITY 木と相補的な双対実行可能解が存在する  命題 木の使用グラフと相補性: 以下の性質をもつ木

    T が存在する 主問題のある実行可能解 P と双対問題の実行可能解 (f, g) が 存在し、P が T を使用グラフとしてもち、(f, g) は T について相補的  証明(ざっくり): P を実行可能領域の頂点にある最適解の中から適当にとる このとき P は頂点性よりある木 T を使用グラフとしてもつ P は最適解なので相補的な双対実行可能解 (f, g) が存在する f i + g j = C ij を満たす (i, j) ∈ T は森 F になる F が連結であれば F = T となり終了 つづく
  55. 89 / 270 KYOTO UNIVERSITY 木と相補的な双対実行可能解が存在する: 証明続き  証明(つづき): F

    の異なる連結成分を結ぶ組 (i, j) の中で s(i, j) = C ij - f i - g j が最も小さいものを取ってきて i*, j* とおく i* が繋がっている側の連結成分に含まれる倉庫側の変数を s(i*, j*) 増加させ、工場側の変数を s(i*, j*) 減少させる F の連結成分内の倉庫・工場対では f i + g j の値は変わらず、この 連結成分と他の連結成分の対では s(i, j) だけ減少 or 増加する 連結成分間では少なくとも s(i*, j*) だけ余裕があったので 依然として実行可能 f i* + g j* = C i*j* なので F に (i, j) を追加しこのプロセスを 繰り返すことで連結成分を 1 つづつ減らして F を木にできる
  56. 90 / 270 KYOTO UNIVERSITY 木と相補的な双対実行可能解が存在する 緑辺は正の輸送量とする このとき相補性よりタイト 灰辺は輸送量ゼロ このときタイトとは限らない

    緑辺と灰辺を合わせて使用グラフ 連結成分全体の値を 定数変化させてもその中 では制約は守られたまま (cf 双対問題の自由度) 限界まで変化させると 他の連結成分との間が タイトになり繋がる
  57. 91 / 270 KYOTO UNIVERSITY 木と相補的な双対解はほぼ一意に定まる  相補解の一意性: 木 T

    を定めると T と相補的な f 1 = 0 を満たす双対問題の解 がただ 1 つ存在する  証明(ざっくり): T 上で倉庫 i = 1 と繋がっている工場 j の値は相補性より g j = C 1j - f 1 = C ij と定まる 同様に j と繋がる倉庫も定まり芋づる式に全ての値が定まる 閉路が無いので相補性の条件が矛盾することはなく最後まで定まる  T に含まれない辺については f i + g j > C ij となってしまうことがある ので実行可能解とはなるとは限らないことに注意  固定するのは f 1 でなくてよいが双対解の自由度より何かは固定する 実行可能とは限らない
  58. 92 / 270 KYOTO UNIVERSITY 木と相補的な双対解はほぼ一意に定まる: 例  例: 使用グラフ(木)

    300 100 50 80 100 50 200 f 1 辺上の値は輸送コスト 輸送量ではない
  59. 94 / 270 KYOTO UNIVERSITY 木と相補的な双対解はほぼ一意に定まる: 例  例: 300

    100 50 80 100 50 200 0 200 50 100 200 - 0 = 200 50 - 0 = 50 100 - 0 = 100
  60. 95 / 270 KYOTO UNIVERSITY 木と相補的な双対解はほぼ一意に定まる: 例  例: 300

    100 50 80 100 50 200 0 200 50 100 -20 -50 80 - 100 = -20 50 - 100 = -50
  61. 96 / 270 KYOTO UNIVERSITY 木と相補的な双対解はほぼ一意に定まる: 例  例: 300

    100 50 80 100 50 200 0 200 50 100 -20 -50 150 350 100 - (-50) = 150 300 - (-50) = 350
  62. 98 / 270 KYOTO UNIVERSITY 相補性を使って途中で最適性を判断できる  ちょっと改善した愚直なアルゴリズム: 倉庫と工場を結ぶ木を全て考え 各木

    T について T を使用グラフとしてもつ主問題の解 P と T と相補的な双対問題の解 (f, g) を計算し、双方が実行可能で あれば (P, f, g) を最適解として出力する  そのような (P, f, g) は存在するので解は出力される  (T, f, g) が相補的 -> (P, f, g) が相補的 -> (P, f, g) は最適  以前の愚直アルゴリズムは全ての T について計算して最小を出力 したが、改善バージョンは実行可能解が見つかった時点でその後を 計算せずとも直ちに最適とわかり停止できる点で効率がよい
  63. 99 / 270 KYOTO UNIVERSITY NS法は主問題の頂点を辿りつつ相補的な双対解を探す  ネットワークシンプレックス法の基本アイデア: 主問題の実行可能解 P

    と P に付随する使用木 T の組 (P, T) を 変化させていって T と相補的な双対実行可能解 (f, g) を見つける つまり主問題の実行可能領域の頂点を辿っていって
  64. 100 / 270 KYOTO UNIVERSITY まず初期頂点解を見つける  ステップ 1: 何でもよいので頂点解を見つける

     North west corner rule というアルゴリズムで見つけられる: 倉庫 1 から、工場 1 に目一杯渡し、工場 2 に目一杯渡し ... を 繰り返し、全部尽きた時点で今度は尽きた場所の工場からはじめて 倉庫 2 から目一杯貰い、倉庫 3 から目一杯貰い ... を繰り返し、 全部満たされた時点で今度は満たされた場所の倉庫からはじめて 次の工場に目一杯渡し、その次の工場に目一杯渡し... を繰り返す  詳しくは次のスライドからはじまるアニメーションを参照
  65. 101 / 270 KYOTO UNIVERSITY 頂点解の見つけ方: 例  例: 0.1

    0.2 0.5 0.1 0.1 0.6 0.1 0.3 まずは工場 1 から
  66. 102 / 270 KYOTO UNIVERSITY 頂点解の見つけ方: 例  例: 0.1

    0.2 0.5 0.1 0.1 0.6 0.1 0.3 0.1 まで受け入れ可
  67. 103 / 270 KYOTO UNIVERSITY 頂点解の見つけ方: 例  例: 0.1

    0.2 0.5 0.1 0.0 0.6 0.1 0.2 0.1 -0.1 目一杯輸送する
  68. 106 / 270 KYOTO UNIVERSITY 頂点解の見つけ方: 例  例: 0.1

    0.2 0.5 0.0 0.0 0.6 0.1 0.1 0.1 0.1 0.5 まで受け入れ可 0.1 まで輸送可
  69. 107 / 270 KYOTO UNIVERSITY 頂点解の見つけ方: 例  例: 0.1

    0.2 0.4 0.0 0.0 0.6 0.1 0.0 0.1 0.1 0.1 全部輸送 まだ余裕あり
  70. 108 / 270 KYOTO UNIVERSITY 頂点解の見つけ方: 例  例: 0.1

    0.2 0.4 0.0 0.0 0.6 0.1 0.0 0.1 0.1 0.1 次はこちら側を中心に考える
  71. 109 / 270 KYOTO UNIVERSITY 頂点解の見つけ方: 例  例: 0.1

    0.2 0.3 0.0 0.0 0.6 0.0 0.0 0.1 0.1 0.1 0.1
  72. 110 / 270 KYOTO UNIVERSITY 頂点解の見つけ方: 例  例: 0.1

    0.2 0.3 0.0 0.0 0.6 0.0 0.0 0.1 0.1 0.1 0.1
  73. 111 / 270 KYOTO UNIVERSITY 頂点解の見つけ方: 例  例: 0.1

    0.2 0.0 0.0 0.0 0.3 0.0 0.0 0.1 0.1 0.1 0.1 0.3 容量を全部使い果たした 次はこちら側を中心に考える
  74. 112 / 270 KYOTO UNIVERSITY 頂点解の見つけ方: 例  例: 0.1

    0.2 0.0 0.0 0.0 0.3 0.0 0.0 0.1 0.1 0.1 0.1 0.3
  75. 113 / 270 KYOTO UNIVERSITY 頂点解の見つけ方: 例  例: 0.1

    0.0 0.0 0.0 0.0 0.1 0.0 0.0 0.1 0.1 0.1 0.1 0.3 0.2 -0.2
  76. 114 / 270 KYOTO UNIVERSITY 頂点解の見つけ方: 例  例: 0.1

    0.0 0.0 0.0 0.0 0.1 0.0 0.0 0.1 0.1 0.1 0.1 0.3 0.2
  77. 115 / 270 KYOTO UNIVERSITY 頂点解の見つけ方: 例  例: 0.0

    0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.1 0.1 0.1 0.3 0.2 0.1 ∑ a i = ∑ b j の仮定より 必ず最後は辻褄があう
  78. 116 / 270 KYOTO UNIVERSITY 頂点解の見つけ方: 例  例: 頂点解とそれに付随する木が見つかった

    これが探索の初期解 0.1 0.2 0.5 0.1 0.1 0.6 0.1 0.3 0.1 0.1 0.1 0.1 0.3 0.2 0.1
  79. 117 / 270 KYOTO UNIVERSITY 付随する双対解が実行可能なら終了、不能なら次へ  ステップ 2: 今の木に付随する

    f 1 = 0 な双対解 (f, g) を計算する  もしこれが実行可能ならば最適解なのでこれで終了  実行可能でないなら f i* + g j* > C i*j* なる違反 (i*, j*) が存在する  木辺上では f i + g j = C ij となっているので (i*, j*) は木の辺でない  (i*, j*) 上で双対問題の制約を満たすようにしたい  ステップ 3: 次に試す木は今の木とほとんど同じだが、(i*, j*) を 使用グラフの辺に含めることで f i* + g j* = C i*j* が満たされるようにする  どうすればそのような (P, T) の組を見つけられるか 木に付随する双対解は 芋づる式に計算できる
  80. 118 / 270 KYOTO UNIVERSITY 含めたい辺を木に含めるにはどうしたらよいか  (i*, j*) を

    T に含める  すると T には閉路が一つできる  倉庫 → 工場の辺 (i*, j*) からはじめて、工場 → 倉庫の辺、 倉庫 → 工場の辺、... と一周まわるまで辿る  偶数番目の辺(工場 → 倉庫の辺)に相当する P ij のうち最小値 を p とおき添字を (i’, j’) とする  偶数番目の辺に相当する P ij を一律 p 減らし、奇数番目の辺に 相当する P ij を一律 p 増やす そして (i’, j’) を T から削る → T が木になる p 一番輸送が少ない辺 (i*, j*)
  81. 120 / 270 KYOTO UNIVERSITY ネットワークシンプレックスの更新の仕方 0 +p  こう

     こうしてできた P は実行可能で T を使用グラフとしてもつ、なぜか? +p +p -p -p -p
  82. 121 / 270 KYOTO UNIVERSITY 新しい木と解のペアも頂点に対応している  実行可能性: 閉路の中で一番輸送量が少ない箇所と同じぶんだけ減らしたので P

    ij は非負のまま 閉路の偶数番目と奇数番目を同じ分だけ増減させたので 余りなし制約、不足なし制約は依然満たされている  付随性: 除かれた辺 (i’, j’) は もともとの輸送量 p だけ輸送量が 減らされたので変化後は 輸送量ゼロになっており大丈夫 Gabriel Peyré, Marco Cuturi. Computational Optimal Transport. 2019.
  83. 123 / 270 KYOTO UNIVERSITY NS法のステップで目的関数は広義単調減少  単調性: この操作で主問題の目的関数は広義単調に減少していく 

    証明: (i*, j*) を含めたときにできた閉路を (i*, j*) からはじめて、辿った ノードを順番に並べる: i 1 = i*, j 1 = j*, i 2 , j 2 , ..., i l , j l , i l+1 = i* 木上では f i + g j = C ij が成り立っていることに注意すると、 (奇数番目の辺のコスト合計)- (偶数番目の辺のコスト合計) = この量 × p (≧ 0) だけ目的関数が減る ほとんど全ての f i と g j  がキャンセル (i*, j*) は違反辺 p 一番輸送が少ない辺 (i*, j*) i 1 i 2 i 3 j 1 j 2 j 3
  84. 124 / 270 KYOTO UNIVERSITY 改善する方向に頂点を辿っていくのでシンプレックス法  木は頂点解に対応する 辺を一つ削除して別の辺を追加した木は隣接頂点となる 

    つまり、ネットワークシンプレックス法は目的関数が悪化しない方向に 隣接する頂点を辿っていく これがシンプレックス法という名前の由来  頂点を移動したとき目的関数が変化しないこともある 追加する辺を適当に選び続けると同じ解をぐるぐる回る可能性がある これは通常のシンプレックス法と同様 無限ループを避けるうまいピポット選択(カニンガムのルール)が必要 「組合せ最適化」 9.6 章 ネットワーク単体法や 「グラフ・ネットワーク・組合せ論」 2.3.2 章などを参照
  85. 125 / 270 KYOTO UNIVERSITY NS法は最悪計算量は悪いが実際は高速  ピボット選択をうまくしてもネットワークシンプレックス法の最悪計算量は 多項式時間で抑えられない これは通常のシンプレックス法と同様

     ネットワークシンプレックス法は多くのインスタンスについて非常に高速 これも通常のシンプレックス法と同様  ネットワークシンプレックス法は頂点の性質をうまく使って高速に頂点を 辿れるので通常のシンプレックス法よりもさらに高速  インスタンスに特殊な性質(例えば低次元性)がないかつ CPU 上 で解くのであればだいたい最速  例えばグラフアルゴリズムライブラリ Lemon や networkx 最適輸送ライブラリ POT で使える 多項式で抑えられる 変種もある
  86. 126 / 270 KYOTO UNIVERSITY NS法は主問題を解いた。次からは双対問題を紹介  ネットワークシンプレックス法は主問題の解を実行可能領域内で 辿りつつ、主問題の目的関数を減少させていく そして相補的な双対問題が実行可能になれば終了

     以上の点でネットワークシンプレックス法は主問題のアルゴリズム  以降は双対問題を更に掘り下げ、双対問題を直接解いていく アルゴリズムを紹介する
  87. 127 / 270 KYOTO UNIVERSITY この章のゴールの確認  最適輸送問題において多面体の頂点が表している解を知る 使用グラフが木であるときかつそのときのみ頂点にある 

    多面体の頂点を効率よく辿る方法を知る 違反している双対問題の制約を見つけ、その辺が使用グラフに 含まれるように変更すると「より良い」頂点に移ったことになる
  88. 132 / 270 KYOTO UNIVERSITY 双対問題の最適解は各制約の感度を表す  まず、通常の線形計画と同様に、双対問題は感度分析を提供する  主問題の制約に対応する双対変数の最適解はその制約が

    少し大きくなった時に最適解がどれほど大きくなるかを表している  たとえば、(f*, g*) を双対問題の最適解とすると i 番目の倉庫に入っ ている小麦粉の量が ε 増えたときに ε f* i くらい輸送コストが増え、 j 番目の工場が要求する小麦粉の量が ε 増えたときに ε g* j くらい 輸送コストが増えることを表す  直感的には、ちょっとくらい小麦粉の量が変わっても双対問題の解は 微小しか変わらない。簡単化して解が同じだとすると、双対問題の 目的関数は a i が ε 増えると f* i だけ増える
  89. 133 / 270 KYOTO UNIVERSITY 例: 僻地にある倉庫では足元見られて処分代が高くなる  例: C

    1 2 3 1 100 100 100 2 1 1 1 3 1 1 1 コスト大 f 1 + g i ≦ C 1j = 100 f 1 大 大きくできる = 流通業者は足元を見ることができる = ちょっとでも小麦粉が増えるとコスト激増 こちらは他の制約 f 2 + g j ≦ C 2j = 1 があるので大きくできない
  90. 135 / 270 KYOTO UNIVERSITY 双対変数は負号を付けて処分額を販売額にしてよい  問: g j

    の「小麦粉の販売額」は分かるが f i の「小麦粉の処分額」は 不自然では?普通小麦粉の処分にお金を払わないだろ。 倉庫では小麦粉を周りに売って、工場では小麦粉を周りから買う方が 問題として自然では?  答: f i は負の値も取る。 -f i 円で小麦粉を売っていると考えられる。  f’ i = -f i と変数変換して、流通業者が f’ i 円で小麦粉を買い取ると してもよい。このとき、流通業者の利益は これを最大化したい 仕入れ経費 売上
  91. 136 / 270 KYOTO UNIVERSITY 負号の付け方で様々な等価な双対問題が考えられる (1) f: 倉庫での処分コスト, g:

    工場での購入コスト (2) f’: 倉庫での下取り利益, g: 工場での購入コスト (3) f: 倉庫での処分コスト, g’: 工場での引取利益 (4) f’: 倉庫での下取り利益, g’: 工場での引取利益  全て等価(最適値は同じ、最適解も簡単に変換できる)なので 一番イメージしやすいものを念頭に置くのが良い  式としては (1) と (3) がこのスライドでも文献でもよく登場する
  92. 139 / 270 KYOTO UNIVERSITY 大小を上下と思うと f は上に g’ は下に最適化する

     各変数について、値が大きいほど「上」、小さいほど「下」にあるという イメージを持ってみる  a i > 0, b j > 0 より、f i は上にあればあるほど、g’ j は下に あればあるほどよい  f i を 1 センチ上に上げると a i ポイント貰える g’ j を 1 センチ下に下げると b j ポイント貰える
  93. 140 / 270 KYOTO UNIVERSITY f i と g j

    のペアを C ij 以上離しすぎると制約違反  ただし、全く自由に f i と g’ j を動かせる訳ではない f i g’ j < C ij OK f i OK f i g’ j = C ij OK f i g’ j > C ij NG 1 つでもこういう組が できるとアウト
  94. 141 / 270 KYOTO UNIVERSITY 双対問題は紐が付いたボールを上下に引っ張るイメージ  双対問題の物理的イメージ: f に対応する赤玉

    n 個と g’ に対応する青玉 m 個を用意する f i に対応する赤玉と g’ j に対応する青玉の間に長さ C ij の紐を張る 合計 nm 本の紐がある 紐が切れないように赤玉を上に引っ張り、青玉を下に引っ張り できるだけポイントを稼ぎたい f i を 1 センチ上に上げると a i ポイント貰える g’ j を 1 センチ下に下げると b j ポイント貰える 相補性より「紐がピンと張っている組」の間に輸送が生じる f i g’ j > C ij f i g’ が上にある場合はセーフなので正確には g’ が 上にある限りは切れない魔法の紐を考える g’ は下に引っ張るのでこういうケースは起きづらい OK
  95. 143 / 270 KYOTO UNIVERSITY C-transform は f を固定して g’

    を下に下げる変換  ここでも 3 番目の定式化 (f, g’) の最適化を考える  C-transform は今手元にある双対問題の解 (f, g’) から より良い新しい解を生成する方法  イメージとしては、f に対応する赤玉を完全に固定してしまって、 糸が弛んでる g’ に対応する青玉をできるだけ下に引っ張る  赤玉を固定すると青玉どうしは干渉しあわないので並列に引っ張っても (どの順番引っ張っても)結果は同じ  赤玉が下に動くだけなので、目的関数は改善される(か変わらない)
  96. 144 / 270 KYOTO UNIVERSITY C変換は簡単に計算でき、制約は守られる  双対問題の解 f について

    f の C-transform fC ∈ Rm は  f i と繋がっている紐は f i - C ij まで下げられる どれか一つでも切れるとアウトなので、一番保守的な制約を採用する  fC は g’ 側の変数を表していることに注意。(f, fC) が新たな解になる。  明らかに (f, fC) は実行可能  証明:
  97. 145 / 270 KYOTO UNIVERSITY C変換はベストな相方を見つける  (f, g’) が実行可能解のとき

    (f, fC) の目的関数値の方が悪くない  証明:  つまり、(f, fC) は (f, *) という形の実行可能解の中でベストな解  実は、双対問題は f, g’ の両方を自由変数する必要はなく、 f を探索空間として、その都度ベストな g’ を C-transform で計算し その時の値を目的関数とすればよい  が、実際は f, g’ の両方を自由変数とした方が簡単な場合も多い g’ の実行可能性より g’ j ≧ f i - C ij
  98. 146 / 270 KYOTO UNIVERSITY C-transform は g’ を固定して f

    を上に上げる変換  逆に g’ 側を固定して f を上に引っ張る変換も考えられる これを C-transform という  やはり (g’C, g’) は実行可能であり (*, g’) の形の解の中でベスト
  99. 147 / 270 KYOTO UNIVERSITY 交互に変換する戦略は失敗する  適当な f からはじめて、(f,

    fC) → (fCC, fC) → (fCC, fCCC) → ... と交互にベストな相方を見つけていって解を改善していくアルゴリズムが 考えられる  がこれは上手くいかない なぜなら、fC = fCCC となりわずか 2 ステップで更新はストップ 直感的には、fC は f に対応できる一番下であるのに f よりもいくつかの f i が上に動いた fCC を相手にしては fC からさらに下に動くことは不可能 しかも (fCC, fC) が双対問題の最適解とは限らない  C-transform をソフトにして交互に最適解まで変換できるように した(と見ることができる)のが後で紹介する Sinkhorn 法
  100. 148 / 270 KYOTO UNIVERSITY この章のゴールの確認  双対問題の最適解が表しているものを理解する 各小麦粉の量を増やした時に増えるコストの感度を表している 

    双対問題の四つの変種を知る 変数の符号によって計四種類正負の異なる定式化がある  双対問題の物理的イメージを知る 変数をボールに、制約を糸に例えると f を上にひっぱり g を 下にひっぱるようなイメージ  C-transform を知る 片方の変数を固定した時にもう片方を最大限引っ張る操作 これにより解が良くなる(悪くならない)
  101. 150 / 270 KYOTO UNIVERSITY ハンガリアン法は一対一の割り当て問題によく使われる  その2: ハンガリアン法は線形計画の双対を解く 

    ユースケースは限られるが、一対一の割当問題を解くためには よく使われている  長所 : CPU 上だと早い シンプルに実装できる  短所 : 一対一マッチングにしか使えない (対応するように拡張はできる) 応募者 仕事
  102. 152 / 270 KYOTO UNIVERSITY 倉庫と工場を一対一に対応させる割当問題を考える  この章では特殊だが重要な例、割当問題、のみを考える  具体的には

    n = m, a 1 = ... = a n , b 1 = ... = b m を仮定する  倉庫と工場が n 個ずつあり、一対一に対応させる 倉庫 i と工場 j を対応させるのはコスト C ij かかる 最小の合計コストは? 同数・同量なので n, a, b が入力から消えた 問題がシンプルになった
  103. 153 / 270 KYOTO UNIVERSITY 割当問題はさまざまな応用がある  問題を単純化したが、輸送問題以外にも応用先は多い  例えば機械学習ではクラスタリングの精度を計るのに使う

    クラスタ 1 クラスタ 2 クラスタ 3 クラス A クラス B クラス C C 1B = 4: クラスタ 1 が クラス B を表していると したときの誤分類の個数 クラスタリングアルゴリズムは サンプルをまとめるだけで 何のクラスを表しているかは 出力しないので精度をはかる ためには対応関係を明らかに する必要がある 最適割当は ベストシナリオでの 精度に対応する
  104. 154 / 270 KYOTO UNIVERSITY ハンガリアン法が割り当て問題の最も有名なアルゴリズム  割当問題を解く最も有名なアルゴリズムはハンガリアン法  古典的なアルゴリズムだが今でも色んな場面で頻繁に使われる

    Object-Centric Learning with Slot Attention (NeurIPS 2020) 集合予測の問題において出力順番が不同なのでハンガリアン法でマッチング On Testing for Biases in Peer Review (NeurIPS 2019) 査読論文割当システムにハンガリアン法を使う Partially View-aligned Clustering (NeurIPS 2020) 複数ビューのクラスタの対応にハンガリアン法を使う(提案手法はこれの改良版)
  105. 155 / 270 KYOTO UNIVERSITY 割当問題の定式化 主問題 双対問題 今回は a

    i = 1 としておく a i = 1/n としても等価 といった解も実行可能だが、一対一対応を表す順列行列 が最適値を取り、ハンガリアン法もそのような解を見つける  問題の入力: 割り当てコスト行列 C ∈ Rn×m 最適輸送では入力であった a, b は今回定数
  106. 156 / 270 KYOTO UNIVERSITY ハンガリアン法は双対変数を最適化していく  ハンガリアン法は双対変数を実行可能領域に収めながら 実行可能でない主問題の相補的な解を実行可能に近づけていく ↔

    ネットワークシンプレックス法は主変数は実行可能に収めながら 実行可能でない双対問題の相補解を実行可能に近づけていく  ハンガリアン法のループ中、常に成り立つ条件: (1) は実行可能 (2) ならば はタイト ( ) (3) (4) (1), (2) より P が実行可能になれば直ちに最適解  P ij = 1 となっているとき (i, j) がマッチしているという
  107. 157 / 270 KYOTO UNIVERSITY ハンガリアン法は順番にマッチする相手を決めていく  まず f i

    = 0, g’ j =0, P ij = 0 と初期化する  i = 1, 2, ..., n と順番にマッチする相手を決めていく  ノード i から、左から右へのタイトな辺と右から左のマッチしている辺 だけを辿ってたどり着くノード集合を X とする  X の状態によって二種類に場合分け (a) 未マッチな右側の頂点が X に含まれる場合 (b) 未マッチな右側の頂点が X に含まれない場合
  108. 158 / 270 KYOTO UNIVERSITY ハンガリアン法の場合分け 1  (a) 未マッチな右側の頂点が

    X に含まれる場合 このパスの辺は全てタイトである。このパスの P ij を全て反転させる 双対変数は変化せず以下は全て満たされたまま (2) ならば はタイト (3) (4) i がこれによりマッチするので i ← i+1 として次のイテレーションへ i 1 = i j 1 未 済 i 2 j 2 未 済 i 3 j 3 未 i 1 = i j 1 済 未 i 2 j 2 済 未 i 3 j 3 済
  109. 159 / 270 KYOTO UNIVERSITY ハンガリアン法の場合分け 2  (b) 未マッチな右側の頂点が

    X に含まれない場合  X に含まれる左側のノード i から X に含まれない右側のノード j への 余裕度 C ij - (f i - g j ) が一番小さい値を ε > 0 とする  X に含まれるノードに対応する双対変数を全て ε だけ増加させる  i, j が両方 X に含まれる or 両方含まれないときは f i - g’ j 変化なし j だけ X に含まれるなら P ij = 0 かつ f i - g’ j は減る方向に変化 i だけ X に含まれるなら f i - g’ j は増えるが ε のとり方より制約は守る  この変化により余裕度が一番小さかった辺がタイトになった その辺を含めてもう一度 X を計算しなおす(前々スライドへ)
  110. 160 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 0 0

    0 0 0 0 0  例: C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60
  111. 161 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 0 0

    0 0 0 0 0  例: i = 1 まずはこのノードの 相手を探す C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60
  112. 162 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 0 0

    0 0 0 0 0  例: X: 未マッチ右ノードなし → ケース (b) へ C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60
  113. 163 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 0 0

    0 0 0 0 0  例: ケース (b) C ij - f i + g j = 30 最も余裕なし C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60
  114. 164 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 0 0

    0 0 0 0 30  例: ケース (b) タイトに +30 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60
  115. 165 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 0 0

    0 0 0 0 30  例: X: 未マッチ右ノードあり → ケース (a) へ 未 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60
  116. 166 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 0 0

    0 0 0 0 30  例: ケース (a) マッチ C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60
  117. 167 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 0 0

    0 0 0 0 30  例: i = 2 次はこのノードの 相手を探す C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60
  118. 168 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 0 0

    0 0 0 0 30  例: X: 未マッチ右ノードなし → ケース (b) へ C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60
  119. 169 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 0 0

    0 0 0 0 30  例: ケース (b) C ij - f i + g j = 20 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60
  120. 170 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 0 0

    0 0 0 20 30  例: ケース (b) タイトに +20 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60
  121. 171 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 0 0

    0 0 0 20 30  例: X: 未マッチ右ノードなし → ケース (b) へ C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60
  122. 172 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 0 0

    0 0 0 20 30  例: ケース (b) 余裕= 10 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60
  123. 173 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 0 0

    10 0 0 30 40  例: ケース (b) タイトに +10 +10 +10 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60
  124. 174 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 0 0

    10 0 0 30 40  例: X: 未マッチ右ノードあり → ケース (a) へ C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60 未
  125. 175 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 0 0

    10 0 0 30 40  例: ケース (a) 反転 反転 反転 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60
  126. 176 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 0 0

    10 0 0 30 40  例: i = 3 次はこのノードの 相手を探す C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60
  127. 177 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 0 0

    10 0 0 30 40  例: X: 未マッチ右ノードなし → ケース (b) へ C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60
  128. 178 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 0 0

    10 0 0 30 40 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60  例: ケース (b) 余裕= 30
  129. 179 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 0 0

    10 0 30 30 40 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60  例: ケース (b) +30 タイトに
  130. 180 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 0 0

    10 0 30 30 40 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60  例: X: 未マッチ右ノードなし → ケース (b) へ
  131. 181 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 0 0

    10 0 30 30 40 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60  例: ケース (b) 余裕= 10
  132. 182 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 10 0

    20 0 40 40 50 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60  例: ケース (b) +10 +10 +10 +10 +10 タイトに
  133. 183 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 10 0

    20 0 40 40 50 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60  例: X: 未マッチ右ノードあり → ケース (a) へ 未
  134. 184 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 10 0

    20 0 40 40 50 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60  例: ケース (a) 反転
  135. 185 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 10 0

    20 0 40 40 50 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60  例: i = 4 最後はこのノード
  136. 186 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 10 0

    20 0 40 40 50 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60  例: X: 未マッチ右ノードなし → ケース (b) へ
  137. 187 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 10 0

    20 0 40 40 50 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60  例: ケース (b) 余裕= 20
  138. 188 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 10 0

    20 20 40 40 50 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60  例: ケース (b) +20 タイトに
  139. 189 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 10 0

    20 20 40 40 50 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60  例: X: 未マッチ右ノードなし → ケース (b) へ この辺は使えない
  140. 190 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 10 0

    20 20 40 40 50 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60  例: ケース (b)
  141. 191 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 20 0

    30 30 40 50 60 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60  例: ケース (b) +10 +10 +10 +10 +10 タイトに ルーズに
  142. 192 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 20 0

    30 30 40 50 60 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60  例: X: 未マッチ右ノードあり → ケース (a) へ 未
  143. 193 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 20 0

    30 30 40 50 60 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60  例: ケース (a) 反転 反転 反転 反転 反転
  144. 194 / 270 KYOTO UNIVERSITY ハンガリアン法: 例 0 20 0

    30 30 40 50 60 C 1 2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60  例: 主・双対最適解が求まった
  145. 195 / 270 KYOTO UNIVERSITY ハンガリアン法は最悪ケース三乗時間  (b) が起こるたびに X

    に含まれる右ノードが 1 つは増えていくので n 回 (b) が起こるまでには未マッチ右ノードに到達できるようになる  (a) が起こるとマッチ組は 1 組増える  最悪ケースでも n2 ステップ以内に最適解を見つけて停止できる  (a), (b) を効率よく実装すれば最悪計算量 O(n3) 時間で解ける ちゃんと効率良く実装するのはそれなりに難しい  たとえば古いバージョンの scipy などで実装されている 今は同じ関数でも違うアルゴリズムが使われている
  146. 196 / 270 KYOTO UNIVERSITY ハンガリアン法は一般の最適輸送にも拡張できる  今回は簡単のため割り当て問題を考えたが、同様のアルゴリズムは 一般の最適輸送問題にも適用できる 

    たとえば最短パス反復アルゴリズム、主双対アルゴリズム  「組合せ最適化」 9.4 章 最短パス反復アルゴリズム 「グラフ・ネットワーク・組合せ論」 2.3.2 章など参照
  147. 197 / 270 KYOTO UNIVERSITY この章のゴールの確認  最適輸送の特殊ケースである割当問題を知る a と

    b が全て同じな一対一を求める問題を割当問題という クラスタリングの評価などに使える  割当問題を解くハンガリアン法を理解する 実行可能な双対解を動かしながら相補的な主解を 一つずつマッチさせて実行可能にしていく
  148. 199 / 270 KYOTO UNIVERSITY Sinkhorn は微分可能で GPU フレンドリー 

    その3: Sinkhorn アルゴリズムは目的関数に正則化を加え、 行列演算で最適輸送を解くアルゴリズム  Sinkhorn が [Cuturi+ 2013] で紹介されたことで 機械学習界隈で最適輸送が使われるようになるキッカケになった  長所 :コストに関して微分可能 GPU 上だと早い 複数の最適輸送を並列に解ける 最適解が唯一に定まる  短所 : 厳密な最適輸送の解ではない オーバーフローや解の丸めに気をつける必要あり Marco Cuturi. Sinkhorn Distances: Lightspeed Computation of Optimal Transport. NeurIPS 2013. Gabriel Peyré, Marco Cuturi. Computational Optimal Transport. 2019.
  149. 202 / 270 KYOTO UNIVERSITY エントロピー正則化つき問題はより簡単に解ける  つまり、輸送コストを下げつつ、エントロピーが大きい解を求める  ε

    を十分小さく取ると、輸送コストを下げるの項が支配的になり 元の最適輸送問題とほぼ同じ このときエントロピーは解のタイブレークに使われる  これまでは輸送行列 P がスパースな解を求めてきたが、 エントロピー正則化つき問題では最適値が P が密行列になる  目的関数が線形でなくなり、解も密になるので難しい問題になってそう だが、実はエントロピー正則化つき問題はより簡単に解ける
  150. 203 / 270 KYOTO UNIVERSITY エントロピー正則化つき問題の最適解は一意  嬉しさ1 : エントロピー正則化つき問題の目的関数は強凸

     制約が線形制約だけということは実行可能領域は凸なので 最適解は一意に定まる ↔ 線形計画は面全体で最適をとることがあり、一意に定まらない 線形 エントロピーは強凹 → 全体で強凸
  151. 204 / 270 KYOTO UNIVERSITY エントロピー正則化つき問題はコストについて微分可能  嬉しさ2 : コストに対して微分可能

     正則化なしの問題では解がスパースであり、たとえば遠くの全然関係 ない倉庫と工場のコストを少し下げたところで最適解には全く影響を 与えず、有用な勾配情報が得られなかった  エントロピー正則化つきでは、遠くの組であっても最適値で少量は 輸送が生じるので、遠くの組のコストの感度がわかる  これにより、ランダムなコストからはじめて、 最適輸送を解く → 勾配を見て最適値が下がるようにコスト更新 → 最適輸送を解く → コスト更新 → ... というように良いコストの設計を最適化できる もちろん制約がないと全コスト 0 が最適なので、コストに制約は入れる
  152. 205 / 270 KYOTO UNIVERSITY 双対問題を考える。まずはラグランジュ乗数を導入  エントロピー正則化つき問題でも双対を考える  ラグランジュ乗数

    f, g を導入するとラグランジュ関数は以下になる これを P ij について微分すると イコールゼロとおくと L が最小となる P は以下の通り P ij ≧ 0 についての乗数は あとでどうせ消えるので省略
  153. 206 / 270 KYOTO UNIVERSITY 双対問題は正則化なし版の制約をソフトにしたもの  前ページの最適 P* をラグランジュ関数に代入すると双対問題は

     エントロピー正則化つき最適輸送の双対問題は制約なし最大化問題  より最終項は f i + g j ≦ C ij のソフト版 ε → 0 で厳密制約になる 正則化なしと同じ ソフト版双対制約 嬉しい
  154. 207 / 270 KYOTO UNIVERSITY 双対問題の目的関数の勾配は f 内 g 内で絡みなし

     双対問題の目的関数を D をおく  D の f と g についての勾配は以下の通り  f i の勾配の中に f k (k ≠ i), g j の勾配の中に g k (k ≠ j) は 出てこない 嬉しい
  155. 208 / 270 KYOTO UNIVERSITY f ごと、g ごとに最適値が厳密に求まる  Sinkhorn

    アルゴリズムの基本的な考えは座標向上法  f を固定したときの g の最適値は勾配イコールゼロとおいて g を固定したときの f の最適値も同様に  Sinkhorn アルゴリズムは f 固定での g の厳密最適化 → g 固定での f の厳密最適化を交互に繰り返す
  156. 209 / 270 KYOTO UNIVERSITY Sinkhorn は f と g

    を交互に最適化する  対数領域での Sinkhorn アルゴリズム  ステップ 1: f(1) を適当に初期化(例えばゼロベクトル), t = 1 に  ステップ 2:  ステップ 3:  ステップ 4: f と g が収束するまで t ← t + 1 でステップ 2 へ  目的関数が微分可能で各ステップ唯一解なので大域最適に収束
  157. 210 / 270 KYOTO UNIVERSITY Sinkhorn はソフト版 C-transform  f

    と g の一方を固定してもう一方を最適化するのは C-transform に 似ている  一般に x 1 , ..., x n ∈ R について log sum exp は  Sinkhorn のイテレーションはソフト版 C-transform と考えられる  交互 C-transform は大域収束しなかったが Sinkhorn は強凸 + 制約なしより大域収束 f i - C ij が一番大きい = C ij - f i が一番小さい = 一番制約が厳しい ほぼ max 関数 ε → 0 で厳密に C-transform
  158. 212 / 270 KYOTO UNIVERSITY ギブスカーネルの定義  さらに K ∈

    Rn×m を以下で定める  K をギブス (Gibbs) カーネルという  i と j の類似度を表している  例えば C がユークリッド距離の自乗のときガウスカーネルになる コストに負号 → 類似度 exp は単調 Sinkhorn のココに登場
  159. 213 / 270 KYOTO UNIVERSITY Sinkhorn のイテレーションは行列積でかける  ギブスカーネル行列を使うとさらにシンプルに 

    ただし KT は転置行列、割り算は要素ごと、分母は行列ベクトル積  非常に単純に実装できる  行列ベクトル積がメインなので GPU で高速計算可能 ふつう Sinkhorn というと この形の更新アルゴリズムを指す
  160. 214 / 270 KYOTO UNIVERSITY Sinkhorn 変数から元の変数への戻し方  Sinkhorn により

    (u, v) が求まると、変数変換の式より元の 双対変数は  主問題の変数はラグランジュ関数から双対関数を求めた式より f, g が完璧に収束しない限り、求めた P ij は厳密には実行可能とは 限らないことに注意  この P ij から違反分をいい感じに分配して実行可能解を計算する アルゴリズムも提案されている [Altschuler+ 2017] Jason Altschuler, Jonathan Weed, Philippe Rigollet. Near-linear time approximation algorithms for optimal transport via Sinkhorn iteration. NeurIPS 2017.
  161. 215 / 270 KYOTO UNIVERSITY Sinkhorn は正則化が小さいとき数値的に不安定  注意: 変数変換により計算を指数領域で行っている

     特に、ε が小さいとこれらの値は非常に大きくなり、オーバーフローを 起こす可能性がある  数値的に不安定な時は、log sum exp などの計算の分遅いが 対数領域で計算を行う必要がある
  162. 216 / 270 KYOTO UNIVERSITY Sinkhorn アルゴリズム: 例 C 1

    2 3 4 1 30 80 40 90 2 20 50 90 80 3 80 70 30 40 4 70 50 10 60  例: ハンガリアン法で使った同じ例を使う 以下はハンガリアン法で求めた最適解
  163. 217 / 270 KYOTO UNIVERSITY Sinkhorn アルゴリズム: 例  例:

    入力を入れてもわずか 25 行 メインロジックはわずか 3 行 マッチングの可視化 1 -- 1 2 -- 2 3 -- 4 4 -- 3 となりほぼ正則化なしと同じ 正則化のため他のマスも少し 正の値をとる
  164. 218 / 270 KYOTO UNIVERSITY 微分の使いどころ: 配置問題  微分の使いどころの例: あなたは工場会社の社長で

    m 個の工場を持っている j 番目の工場の位置は (x j , y j ) で小麦粉を b j 要求する n 個の倉庫を建設する。i 番目の倉庫は小麦粉を a i 保存できる 小麦粉を 1 グラム 1 メートル輸送するのに 1 円かかる 倉庫をどこに建設すればよい?  解き方: (1) 倉庫の位置変数をランダムに初期化 (2) C ij = 倉庫 i と工場 j の距離として計算 (3) Sinkhorn に (a, b, C) を入力して解いて総輸送コストを評価 (4) 倉庫の位置変数の勾配を求めて倉庫の位置を改善 → (2) へ ? ? ? どこに置く?
  165. 219 / 270 KYOTO UNIVERSITY 最新の手法でも配置問題に Sinkhorn が使われる  この問題はクラスタリングと見ることもできる

    [Cuturi+ ICML 2014]  工場がサンプル点で、倉庫がクラスタ中心 クラスタが含むサンプル数(倉庫に入る小麦粉)は決まっている  この発想から、クラスタリングの際に Sinkhorn 法を使う手法が 最新の機械学習でも使われている [Caron+ NeurIPS 2020] Marco Cuturi, Arnaud Doucet. Fast Computation of Wasserstein Barycenters. ICML 2014. Mathilde Caron, Ishan Misra, Julien Mairal, Priya Goyal, Piotr Bojanowski, Armand Joulin. Unsupervised Learning of Visual Features by Contrasting Cluster Assignments. NeurIPS 2020. クラスタ中心
  166. 220 / 270 KYOTO UNIVERSITY 最新の手法でも公平予測問題に Sinkhorn が使われる  他にも、二つの集合の距離を近づけるための微分可ロスとして使われる

     例えば、公平性の担保のため、男性についての予測と女性についての 予測分布が同じようにしたい  予測誤差 + 赤と青の最適輸送距離を最小化 [Oneto+ NeurIPS 2020] Luca Oneto, Michele Donini, Giulia Luise, Carlo Ciliberto, Andreas Maurer, Massimiliano Pontil. Exploiting MMD and Sinkhorn Divergences for Fair and Transferable Representation Learning. NeurIPS 2020. ニューラルネットワーク 入力 d 次元ベクトルが 出てくる Rd 予測器 出力 各丸は各サンプルの 埋め込み 赤: 女性サンプル 青: 男性サンプル
  167. 221 / 270 KYOTO UNIVERSITY この章のゴールの確認  エントロピー正則化最適輸送の定義と嬉しさを知る 「ぼやけた」輸送プランを好むように負のエントロピーが ペナルティとして加わる

    これにより最適解が一意になり微分可能になる  Sinkhorn アルゴリズムの動作と嬉しさを知る 双対問題の f と g を交互に最適化 f と g を変数変換することでより簡単な式で書けるようにする 実装が簡単かつ GPU で高速に計算できるようになる
  168. 223 / 270 KYOTO UNIVERSITY ニューラルネットワークによる推定は GAN の基本技術  その4:

    ニューラルネットワークによる推定は双対問題の特殊形を ニューラルネットワークの最適化に帰着する  Wasserstein GAN という名前で GAN(生成モデル)の 生成分布とデータ分布の距離をはかるのに使われている  長所 : コストについて微分可能 分布からのサンプルをミニバッチで繰り返し解くのに 向いている  短所 : 単一の入力には向いていない コストが距離で定義されていないと 使えない Tero Karras, Timo Aila, Samuli Laine, Jaakko Lehtinen. Progressive Growing of GANs for Improved Quality, Stability, and Variation. ICLR 2018.  GAN による顔写真の自動生成
  169. 225 / 270 KYOTO UNIVERSITY 三番目の定式化を考える  この章で考えるのは 3 番目の双対の最適化

    (f, g’) を変数にしたもの この章では主問題は考えない  f: 回収価格(高くしたい) g’: 引取コスト (低くしたい) 3 番目の定式化
  170. 226 / 270 KYOTO UNIVERSITY 倉庫と工場が座標に埋め込まれている場合を考える  倉庫・工場が座標を持っていて、輸送コストがその間の距離で 定義される場合を考える 

    倉庫 の座標を , 工場 の座標を と表し 倉庫 から工場 に 1 グラムの小麦粉を輸送するコストを とする  問題の入力: C は x, y から 陰に定まるので 受け取らない
  171. 227 / 270 KYOTO UNIVERSITY 座標が決まれば双対値は決まる  命題 と の関数性:

    最適解では ならば , ならば  証明: が実行可能で だが だと仮定する 任意の k に対して なので、 更新しても実行可能で、目的関数は良くなる よって最適解では ならば が成り立つ についても同様
  172. 228 / 270 KYOTO UNIVERSITY 関数の最適化問題に置き換える  倉庫や工場の座標について双対解の値は位置に定めてしまってよい  各倉庫や工場について値を持つのではなく、各座標に値を持つように

    問題を変える  流通企業が各倉庫個別に値段を決めるのではなく、全ての土地で 連続的に値段を定めるイメージ 位置 x での小麦粉の仕入れ価格・販売価格を定める ベクトルの最適化ではなく 関数の最適化 実際に取引が発生するのは各倉庫・工場の位置 g’ は微分ではなく 関数の名前なことに注意
  173. 229 / 270 KYOTO UNIVERSITY 最適解はリプシッツ連続  補題 リプシッツ連続性: 最適解では

    が成り立つ  証明: が実行可能で とする と更新しても実行可能、なぜなら これにより目的関数は増加するので、最適値では  についても同様に
  174. 230 / 270 KYOTO UNIVERSITY 最適解はでは f と g’ が一致

     命題 と の一致性: が成り立つ最適解が存在する  証明: は 以外の場所、 は 以外の場所は値によらず目的値は変わらないので では と仮定できる において とする 制約より よって と更新しても実行可能、なぜなら のリプシッツ性より を減少させると目的関数は増加するので最適解は とできる
  175. 231 / 270 KYOTO UNIVERSITY 最適化変数を f のみにできる  つまり、最適化変数から

    g’ を取り除くことができ  直感的には最適解では小麦の仕入れ価格と販売価格は各地で一致  また、制約条件は f がリプシッツであることにほかならない  C-transform の説明の時に、実質 f だけを自由変数とできると説明 したが、本質的にはそれと同じことを言っている C が距離公理を満たす場合は fC が簡単に表せる y 側でも f を使う
  176. 233 / 270 KYOTO UNIVERSITY 一般の距離空間でも最適化を f のみにできる  発展的注意:

    今回は簡単のためコストをユークリッド距離としたが、一般に 倉庫と工場がある空間 X に埋め込まれ、コストが X 上の距離関数 で定義するという拡張が考えられる  各点で同じ値を取る(関数性)・リプシッツ性・f と g’ の一致性の 証明には距離公理の性質しか使っていないので、一般の距離空間の 場合も f: X → R を X 上でのリプシッツ連続関数制約で最適化 する問題に帰着できる  X をグラフとし、グラフ上の最短経路距離コストを使うケースは頻出
  177. 234 / 270 KYOTO UNIVERSITY f は滑らかな範囲で x 上で大きく y

    上で小さくする  あらためて f だけの最適化問題を眺めてみる  f は X = {x 1 , ..., x n } 上で大きければ大きいほど Y = {y 1 , ..., y m } 上で小さければ小さいほどよい  ただし、f は「滑らか」(リプシッツ連続)でなければならない 正 正 小麦を購入する場所では安く 販売する場所では高く設定する
  178. 235 / 270 KYOTO UNIVERSITY 関数最適化のイメージ  例: 赤点は倉庫、青点は工場に対応 x

    i y i 背景が赤いほど f の値が大きい 背景が青いほど f の値が小さい 関数の変化は滑らか
  179. 236 / 270 KYOTO UNIVERSITY 双対問題は 2 クラス分類問題と見ることができる  このように見ると、X

    = {x 1 , ..., x n } を正例、Y = {y 1 , ..., y m } を 負例としたときの 2 クラス分類問題と見ることができる  a i , b i はサンプルの重み  激しく変化する分類器は制約違反 「正則化」としてリプシッツ関数の条件
  180. 237 / 270 KYOTO UNIVERSITY 分類問題をニューラルネットワークに任せる  「分類器」 f をニューラルネットワークでモデリングすることを考える

     ロス関数はシンプルに  f が 1-リプシッツ連続であることを課すのは難しい  ここでは素朴な解決策を 2 つ紹介
  181. 238 / 270 KYOTO UNIVERSITY パラメータの値を無理やり制限してリプシッツ性を課す  解決策 1: weight

    clipping [Arjovsky+ 2017]  訓練の際ニューラルネットワークの各パラメータの絶対値が 定数 γ > 0 を越えるたびに絶対値 γ にクリップする(γ:ハイパラ)  こうすると、ニューラルネットワーク f の表現できる関数は制限され 急激な変化をする関数は表せなくなる  γ の設定次第で、何らかの k > 0 について k-リプシッツであることが 保証できる  k-リプシッツな関数全てを表現できる訳ではないが、ニューラルネットワー クの柔軟性よりそれなりに豊富な関数が表現できることが期待できる Martin Arjovsky, Soumith Chintala, Léon Bottou. Wasserstein GAN. ICML 2017
  182. 239 / 270 KYOTO UNIVERSITY 各点での勾配が 1 から離れるとペナルティ  解決策

    2: gradient penalty [Arjovsky+ 2017]  f が微分可のとき各点での勾配のノルムが 1 以下 ↔ 1-リプシッツ  色んな点で勾配を評価して 1 から離れていたらペナルティを課す をロスとして最適化 Ishaan Gulrajani, Faruk Ahmed, Martín Arjovsky, Vincent Dumoulin, Aaron C. Courville.. Improved Training of Wasserstein GANs. NeurIPS 2017 ハイパラ正則化係数 色んな訓練サンプル 例えば X ∪ Y 勾配を 1 に近づける
  183. 240 / 270 KYOTO UNIVERSITY ニューラルネットワークの定式化は GAN に用いられる  f

    をニューラルネットワークで表すと f(x) は x について微分可能 → 最適輸送コストを下げるには x をどう動かせばよいかが分かる  ニューラルネットワークによる定式化は GAN(生成モデル)の 訓練によく用いられる: Wasserstein GAN  自然サンプルの集合を X ⊂ Rd とする。例えば画像の集合  GAN により、X に似たサンプル集合を生成したい  GAN により生成したサンプルの集合を Y ⊂ Rd とし、X と Y の距離を 最適輸送コストで測る  これを最小化するように GAN を訓練する
  184. 241 / 270 KYOTO UNIVERSITY ニューラルネットワークの推定は大規模データと相性よし  典型的な X は非常に大きい(数百万枚の画像など)

     生成モデルの生成結果 Y の候補は無限にある  ハンガリアン法や Sinkhorn 法では全ての対を一度に解かなければ ならず、非常に大きいデータには適していない  「分類器」 f を分類するという考え方に立てば、 X と Y から一部のデータをミニバッチとして取り出し f を徐々に 訓練していけばよい → ニューラルネットワークによる推定は大規模(or 無限)のデータ について解くのと相性がよい
  185. 242 / 270 KYOTO UNIVERSITY Wasserstein GAN: 生成モデル と f

    の交互訓練  1. 生成モデル G: Rr → Rd をランダムに初期化 これはノイズ z ∈ Rr を受け取り画像 G(z) ∈ Rd を返す 最適輸送問題を解くニューラルネットワーク f をランダムに初期化  2. G を使って画像 Y = {G(z 1 ), G(z 2 ), ..., G(z m )} を生成  3. f が X と Y を分類できるようにミニバッチで訓練して X と Y の(ミニバッチの)最適輸送コストを計算  4. f を使った最適輸送コストが小さくなるように G をアップデート  5. 2 へ戻る。このとき Y は変化するが f を一から訓練するのではなく、 少し G が変わったくらいでは最適な f もほぼ同じなので 以前のパラメータから訓練をスタートする(ホットスタート)
  186. 243 / 270 KYOTO UNIVERSITY 深層畳み込みネットワークを使うことでリアルな画像生成  Wasserstein GAN の論文や後続の論文では

    G や f として 深層畳み込みニューラルネットワークを使う  上は X を顔写真集合として訓練したてできた生成画像 非常にリアルな顔写真が生成できている 上の画像は段階的学習など別のテクニックも駆使して訓練されている Tero Karras, Timo Aila, Samuli Laine, Jaakko Lehtinen. Progressive Growing of GANs for Improved Quality, Stability, and Variation. ICLR 2018.
  187. 244 / 270 KYOTO UNIVERSITY この章のゴール  双対問題を関数の最適化に置き換える視点を知る 倉庫と工場に座標が付いていてコストが座標から決まる場合は 個々の変数を決めるのではなく空間全体の値を関数として

    求める見方もできる  関数の最適化問題は分類問題と見られることを理解する 関数を決める問題だと思うと、倉庫点を正例、工場点を負例とする 二クラス分類問題と見ることができる このとき関数はリプシッツ性という滑らかさが正則化として要求される
  188. 246 / 270 KYOTO UNIVERSITY スライス法は簡単かつ高速に計算ができる  その5: スライス法は一次元上の最適輸送が簡単に解けることを 利用する

     非常に簡単であることから、とりあえず最適輸送っぽいものを得たい ときには有用  長所 : 実装が簡単 非常に高速  短所 : 元の最適輸送問題の良い近似に なっているとは限らない https://www.programmersought.com/article/67174999352/
  189. 248 / 270 KYOTO UNIVERSITY 倉庫と工場が一次元に埋め込まれている場合を考える  この章では、倉庫と工場が一次元上に埋め込まれている場合を考える  主問題を考える

     入力: a 1 , ..., a n ∈ R, b 1 , ..., b m ∈ R, x 1 , ..., x n ∈ R, y 1 , ..., y m ∈ R 輸送コストが座標の差の絶対値 R x 1 x 2 x 3 y 1 y 2 y 3
  190. 249 / 270 KYOTO UNIVERSITY 一次元の場合は簡単に解ける  このとき最適輸送問題は非常に簡単に解ける  やり方:

    (1) まだ小麦粉が余っている倉庫のうち一番左にあるものを選ぶ (2) まだ小麦粉を要求している工場のうち一番左にあるものに  目一杯小麦粉を輸送する (3) 以上を全ての輸送が完了するまで行う  つまり貪欲法
  191. 254 / 270 KYOTO UNIVERSITY 実装も非常に簡単かつ高速  実装も非常にシンプルかつ高速  (1)

    x 1 , ..., x n  と y 1 , ..., y m をそれぞれソート (2) i = 1, j = 1 に初期化 (3) 倉庫 i から工場 j に min(a i , b j ) 輸送 (4) a i = 0 なら i ← i + 1 として (3) へ (5) a i > 0 なら j ← j + 1 として (3) へ  ソートがボトルネックなので O(n log n) 時間  n = m, a i = b j = 1(割当問題)のとき全て一対一でマッチするので x 1 , ..., x n , y 1 , ..., y m をソートして で一発で求まる
  192. 255 / 270 KYOTO UNIVERSITY 高次元でも一次元に帰着して簡単にする  一次元の場合非常に簡単なことを利用して、 高次元の場合も一次元に帰着する方法が提案されている 

    Sliced Wasserstein distance [Rabin+ 2011]: x 1 , ..., x n ∈ Rd, y 1 , ..., y m ∈ Rd のとき、ランダムな方向 Θ ∈ Rd に x 1 , ..., x n , y 1 , ..., y m を射影し、一次元にして輸送距離を計算 最終的な答えは何度か方向を変えた時の平均値  k 個方向を使うと O(k n log n) 時間 実装は一次元の場合と同等 Julien Rabin, Gabriel Peyré, Julie Delon, Marc Bernot. Wasserstein Barycenter and its Application to Texture Mixing. 2011. https://www.programmersought.com/article/67174999352/ ただしこの期待値が 元の空間でのコストの 近似になっている訳ではない
  193. 256 / 270 KYOTO UNIVERSITY 直線に限らず、どんな一次元に射影してもよい  Sliced Wasserstein distance

    [Rabin+ 2011] では直線に 沿って射影したが、一次元に射影できたらどういう射影でもよい  例えばデータが U 字型だったら U 字に沿った一次元、 一般にはデータのサポートに沿った一次元が良さそう  Generalized sliced Wasserstein [Kolouri+ NeurIPS 2019] 適当な射影関数 Θ: Rd → R を使って射影して一次元で解く ランダムな射影の平均だけでなく、一番青と赤が見分けられるような 射影を学習により見つける方法も提案(判別分析のようなイメージ) Julien Rabin, Gabriel Peyré, Julie Delon, Marc Bernot. Wasserstein Barycenter and its Application to Texture Mixing. 2011. Soheil Kolouri, Kimia Nadjahi, Umut Simsekli, Roland Badeau, Gustavo K. Rohde. Generalized Sliced Wasserstein Distances. NeurIPS 2019.
  194. 257 / 270 KYOTO UNIVERSITY 一次元から木への一般化  貪欲法で解けるのは一次元の場合だけでない 一般に、空間が木の場合でも貪欲法で厳密に解ける 

    辺に長さの情報が付いている木があり、ノードには倉庫と工場がある 倉庫 i から工場 j への輸送コストはその倉庫と工場があるノード を結ぶパスの長さの総和とする  これは一次元の場合の一般化 一次元はパスグラフ(一直線の木) に対応している 50 60 10 40 70 20 30 50 のコストは 20 + 30 + 50 = 100 パスグラフ
  195. 258 / 270 KYOTO UNIVERSITY 木の場合も同様に貪欲法で解ける  木の場合も貪欲法で解ける  やり方:

    (1) DFS 帰りがけ順でノードを走査する (2) 現在のノードの子孫でまだ小麦粉が残っている倉庫と  小麦粉を要求している工場があれば輸送する そのような倉庫か工場のどちらかがなくなるまで繰り返し なくなれば次のノードへ  線形時間で計算できる  一次元の場合は、ソートするのが木を構築するところに相当 ソートが終わればあとは線形時間なので、木の場合も一次元の場合も 同じ計算量
  196. 259 / 270 KYOTO UNIVERSITY 高次元でも木に帰着して簡単にする  一次元の場合と同様に、木に帰着する方法が提案されている  Tree-Sliced

    Wasserstein [Le+ NeurIPS 2019]: x 1 , ..., x n ∈ Rd, y 1 , ..., y m ∈ Rd のとき、ランダムな木(例えば クラスタリング木)のノードに x 1 , ..., x n , y 1 , ..., y m を射影し、 木上で輸送距離を計算 最終的な答えは何度か木を変えた時の平均値  k 個木を使うと O(k n log n) 時間  良い木を使うと元の輸送距離の近似にもなる Tam Le, Makoto Yamada, Kenji Fukumizu, Marco Cuturi. Tree-Sliced Variants of Wasserstein Distances. NeurIPS 2019. http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/117-hcpc-hierarchical-clustering-on-principal-components-essentials/
  197. 260 / 270 KYOTO UNIVERSITY 一次元の場合そのものが重要な場合もある  高次元を一次元や木に帰着する場合を主に扱ったが、一次元の 最適輸送そのものが重要なこともある 

    就活生の特徴量を受け取り自社の適正スコアを出す問題を考える 男性の適正スコアと女性の適正スコアは偏ってはいけないので 最適輸送コストを正則化に加える  データの方を男性と女性の見分けが付かなくなるように最適輸送コスト で補正するやり方が [Feldman+ KDD 2015] で提案されている 適正スコア x 1 x 2 x 3 y 1 y 2 y 3 Michael Feldman, Sorelle A. Friedler, John Moeller, Carlos Scheidegger, Suresh Venkatasubramanian. Certifying and Removing Disparate Impact. KDD 2015.
  198. 261 / 270 KYOTO UNIVERSITY この章のゴールの確認  一次元と木の最適輸送問題が簡単に解けることを知る コストが一次元や木で定義される場合は貪欲法で厳密に解ける 

    高次元の問題も一次元や木に帰着できることを知る 空間が一次元や木と限らなくても、空間を「スライス」することで 一次元の問題に帰着できる場合がある これにより簡単・高速に問題が解ける
  199. 263 / 270 KYOTO UNIVERSITY 今日は最適輸送の色んな解き方をみた  今日は最適輸送のアルゴリズムについてみた  その1:

    線形計画の定式化とシンプレックス法  その2: ハンガリアン法  その3: Sinkhorn アルゴリズム  その4: ニューラルネットワークによる推定  その5: スライシング  それぞれ、長所・短所があった
  200. 264 / 270 KYOTO UNIVERSITY おすすめ本紹介 1  さらに勉強したい人におすすめの本: 

    Gabriel Peyré and Macro Cuturi. Computational Optimal Transport with Applications to Data Science.  arXiv で電子版が公開されています:https://arxiv.org/abs/1803.00567  この講演のネットワークシンプレックス法と Sinkhorn アルゴリズムにつ いてはこの本を参考にした  Sinkhorn アルゴリズムについては特に 収束の話や別の見方など詳しい  動的な定式化・Gromov Wasserstein gradient flow など発展的な話題も多く載っている
  201. 265 / 270 KYOTO UNIVERSITY おすすめ本紹介 2  組合せアルゴリズム的な視点を深堀したい人は以下もおすすめ 

    B.コルテ, J. フィーゲン (著) / 浅野 孝夫, 浅野 泰仁, 小野 孝男, 平田 富夫 (訳) 「組合せ最適化」 → 第 9 章が最適輸送(最小費用フロー)について  藤重 悟 「グラフ・ネットワーク・組合せ論」 → 第 2.3 章が最適輸送(最小費用フロー)について
  202. 266 / 270 KYOTO UNIVERSITY 全体のゴールの確認  解きたい問題に合わせて適切なアルゴリズムを選べるようになる 基本は線形計画とネットワークシンプレックス法 割当問題はハンガリアン法でシンプルに

    コストについて微分したい・GPU を使いたいときは Sinkhorn 法 空間に埋め込まれた分布を最適化するときはニューラルネットワーク シンプルに答えを得たいときはスライス法  最適輸送の主・双対問題を自由自在に使いこなせるようになる 相補的な解を見つける視点 双対問題の物理的イメージと C-transform 双対問題が分類問題のようになっている視点
  203. 268 / 270 KYOTO UNIVERSITY 参考文献  Jason Altschuler, Jonathan

    Weed, Philippe Rigollet. Near-linear time approximation algorithms for optimal transport via Sinkhorn iteration. NeurIPS 2017.  Martin Arjovsky, Soumith Chintala, Léon Bottou. Wasserstein GAN. ICML 2017  Mathilde Caron, Ishan Misra, Julien Mairal, Priya Goyal, Piotr Bojanowski, Armand Joulin. Unsupervised Learning of Visual Features by Contrasting Cluster Assignments. NeurIPS 2020.  Marco Cuturi. Sinkhorn Distances: Lightspeed Computation of Optimal Transport. NeurIPS 2013.  Marco Cuturi, Arnaud Doucet. Fast Computation of Wasserstein Barycenters. ICML 2014.  Michael Feldman, Sorelle A. Friedler, John Moeller, Carlos Scheidegger, Suresh Venkatasubramanian. Certifying and Removing Disparate Impact. KDD 2015.  Ishaan Gulrajani, Faruk Ahmed, Martín Arjovsky, Vincent Dumoulin, Aaron C. Courville. Improved Training of Wasserstein GANs. NeurIPS 2017.
  204. 269 / 270 KYOTO UNIVERSITY 参考文献  Zhenyu Huang, Peng

    Hu, Joey Tianyi Zhou, Jiancheng Lv, Xi Peng. Partially View-aligned Clustering. NeurIPS 2020.  Tero Karras, Timo Aila, Samuli Laine, Jaakko Lehtinen. Progressive Growing of GANs for Improved Quality, Stability, and Variation. ICLR 2018.  Soheil Kolouri, Kimia Nadjahi, Umut Simsekli, Roland Badeau, Gustavo K. Rohde. Generalized Sliced Wasserstein Distances. NeurIPS 2019.  Tam Le, Makoto Yamada, Kenji Fukumizu, Marco Cuturi. Tree-Sliced Variants of Wasserstein Distances. NeurIPS 2019.  Francesco Locatello, Dirk Weissenborn, Thomas Unterthiner, Aravindh Mahendran, Georg Heigold, Jakob Uszkoreit, Alexey Dosovitskiy, Thomas Kipf. Object-Centric Learning with Slot Attention. NeurIPS 2020.  Luca Oneto, Michele Donini, Giulia Luise, Carlo Ciliberto, Andreas Maurer, Massimiliano Pontil. Exploiting MMD and Sinkhorn Divergences for Fair and Transferable Representation Learning. NeurIPS 2020.  Gabriel Peyré, Marco Cuturi. Computational Optimal Transport. 2019.
  205. 270 / 270 KYOTO UNIVERSITY 参考文献  Julien Rabin, Gabriel

    Peyré, Julie Delon, Marc Bernot. Wasserstein Barycenter and its Application to Texture Mixing. 2011.  Ivan Stelmakh, Nihar B. Shah, Aarti Singh. On Testing for Biases in Peer Review. NeurIPS 2019.  B.コルテ, J. フィーゲン (著) / 浅野 孝夫, 浅野 泰仁, 小野 孝男, 平田 富夫 (訳) 「組合せ最適化」  藤重 悟 「グラフ・ネットワーク・組合せ論」