ウォームスタートの学習による L 凸関数最小化の高速化

JCCA-DMIA-2024 ミニシンポジウム（組合せ最適化）の発表スライド．
学会URL：https://sites.google.com/view/jcca-dmia-2024/

Transcript

1. ウォームスタートの学習による L 凸関数最小化の高速化 坂上 晋作（東京大学） 大城 泰平 氏（北海道大学）との共同研究 JCCA 組合せ最適化ミニシンポジウム @ 山形大学 S & Oki. Discrete-convex-analysis-based framework for warm-starting algorithms with predictions, NeurIPS 2022. S & Oki. Rethinking warm-starts with predictions: Learning predictions close to sets of optimal solutions for faster L-/L♮ -convex function minimization, ICML 2023. 論文： https://speakerdeck.com/ssakaue スライド：
2. 概要 2 予測を使ってどれぐらい計算量を改善できるか？ どのくらい問題サンプルがあれば良い予測を学習できるか？ 計算量は初期解と最適解集合の距離に依存 最適解集合 初期解 似た問題を繰り返し解くなら 最適解集合に近い初期解を予測できると良さそう 予測 L 凸関数最小化を最急降下法で解く
3. 3 入力：二部グラフ 𝐺 = 𝐿 ∪ 𝑅, 𝐸 ，重み 𝑤 ∈ ℤ𝐸（ 𝐿 = 𝑅 = 𝑛/2，実行可能） 出力：最大重み完全二部マッチング 𝑀 ⊆ 𝐸 ユーザー 𝐿 アイテム 𝑅 1 2 3 4 5 6 7 8 4 2 3 3 2 2 1 1 1 1 双対 LP（変数 𝑝 = 𝑠, 𝑡 ∈ ℝ𝐿∪𝑅 ） min σ 𝑖∈𝐿 𝑠𝑖 − σ 𝑗∈𝑅 𝑡𝑗 s. t. 𝑠𝑖 − 𝑡𝑗 ≥ 𝑤𝑖𝑗 (𝑖, 𝑗) ∈ 𝐸 主 LP（変数 𝑥 ∈ ℝ𝐸 ） max σ𝑒∈𝐸 𝑤𝑒 𝑥𝑒 s. t. σ 𝑒∈𝑁(𝑣) 𝑥𝑒 = 1 𝑣 ∈ 𝐿 ∪ 𝑅 , 𝑥 ≥ 0 最大重み二部マッチング
4. 4 双対 LP（変数 𝑝 = 𝑠, 𝑡 ∈ ℝ𝐿∪𝑅 ） min σ 𝑖∈𝐿 𝑠𝑖 − σ 𝑗∈𝑅 𝑡𝑗 s. t. 𝑠𝑖 − 𝑡𝑗 ≥ 𝑤𝑖𝑗 (𝑖, 𝑗) ∈ 𝐸 Schrijver. Combinatorial Optimization: Polyhedra and Efficiency. Springer, 2003. cf. Schrijver ’03 (Section 18.5b) ハンガリー法（双対版） 初期実行可能解 𝑝∘ = 𝑠, 𝑡 ∈ ℤ𝐿∪𝑅 を設定 以下を反復： タイト部分グラフ (𝑉, 𝐸∗) を見つける If 𝐸∗ に完マあり（ 𝑆 + 𝑇 = 𝑛/2） 𝐸∗の完全マッチングを一つ出力して終了 (𝑉, 𝐸∗) の最小頂点被覆 (𝑆, 𝑇) を見つける 𝑝 ← (𝑠 + 𝜆𝟏𝑆 , 𝑡 + 𝜆𝟏𝑅∖𝑇 ) ユーザー 𝐿 アイテム 𝑅 1 2 3 4 5 6 7 8 4 2 3 3 2 2 1 1 1 1
5. 5 双対 LP（変数 𝑝 = 𝑠, 𝑡 ∈ ℝ𝐿∪𝑅 ） min σ 𝑖∈𝐿 𝑠𝑖 − σ 𝑗∈𝑅 𝑡𝑗 s. t. 𝑠𝑖 − 𝑡𝑗 ≥ 𝑤𝑖𝑗 (𝑖, 𝑗) ∈ 𝐸 Schrijver. Combinatorial Optimization: Polyhedra and Efficiency. Springer, 2003. cf. Schrijver ’03 (Section 18.5b) ハンガリー法（双対版） 初期実行可能解 𝑝∘ = 𝑠, 𝑡 ∈ ℤ𝐿∪𝑅 を設定 以下を反復： タイト部分グラフ (𝑉, 𝐸∗) を見つける If 𝐸∗ に完マあり（ 𝑆 + 𝑇 = 𝑛/2） 𝐸∗の完全マッチングを一つ出力して終了 (𝑉, 𝐸∗) の最小頂点被覆 (𝑆, 𝑇) を見つける 𝑝 ← (𝑠 + 𝜆𝟏𝑆 , 𝑡 + 𝜆𝟏𝑅∖𝑇 ) 𝑂(𝑚 𝑛) 時間 Hopcroft−Karp の重み無し 最大二部マッチング・最小頂点被覆 総計算量：𝑂(𝑚 𝑛 × 反復回数) 初期解 𝑝∘を上手く選べば高速化可能？
6. 計算量上界 6 実行可能な双対初期解 𝑝∘ ∈ ℤ𝑉 が与えられた場合のハンガリー法の計算量は 𝑂(𝑚 𝑛 × 𝑝∗ − 𝑝∘ 1 ) 𝑝∗ = (ある)双対最適解 実行可能とは限らない予測 Ƹ 𝑝 ∈ ℝ𝑉 が与えられた場合も 𝑂(𝑚 𝑛 × 𝑝∗ − Ƹ 𝑝 1 ) （∵ 任意の Ƹ 𝑝 から 𝑝∗ − 𝑝∘ 1 ≤ 3 𝑝∗ − Ƹ 𝑝 1 を満たす 𝑝∘ を 𝑂(𝑚) 時間の貪欲法で計算可） 事前に 𝑝∗ − Ƹ 𝑝 1 が小さい Ƹ 𝑝 ∈ ℝ𝑉を学習できれば最速で 𝑂(𝑚 𝑛) 時間 補足 • 通常のハンガリー法の最悪時計算量は 𝑂 𝑚𝑛 • 理論最速は Brand et al. ’23 の 𝑚1+𝑜(1)（あまり実用的ではない） Brand et al. A deterministic almost-linear time algorithm for minimum-cost flow. FOCS, 2023 Dinitz et al. Faster matchings via learned duals. NeurIPS, 2021 最大重み二部マッチングの高速化（Dinitz et al. ’21）
7. 7 仮定 • 台集合 𝑉 = 𝐿 ∪ 𝑅 固定（ 𝐿 = 𝑅 = 𝑛/2）かつ実行可能 • ∃𝐶 > 0， −𝐶, +𝐶 𝑛 の中に双対最適解が存在（ 𝑤 ∞ ≤ 𝐶/𝑛 なら成立） • 各インスタンス 𝜋 に対し最適解 𝑝∗(𝜋) ∈ −𝐶, +𝐶 𝑛 を一意に決める（タイブレイク） Dinitz et al. Faster matchings via learned duals. NeurIPS, 2021 最大重み二部マッチングの高速化（Dinitz et al. ’21） 𝑁 十分大で学習した Ƹ 𝑝 を使うと，期待計算量 ≈ 𝑂(𝑚 𝑛 × 𝓓上最良の ℓ𝟏 予測誤差) 学習可能性 𝓓上最良の ℓ𝟏 予測誤差 𝔼𝜋∼𝒟 𝑝∗(𝜋) − Ƹ 𝑝 1 ≤ min ො 𝑝∗∈ −𝐶,+𝐶 𝑛 𝔼𝜋∼𝒟 𝑝∗ 𝜋 − Ƹ 𝑝∗ 1 + 𝜀 仮定を満たすインスタンスが分布 𝒟 から i.i.d. で 𝑁 = ෩ Ω 𝑛3𝐶2/𝜀2 個 与えられれば，高確率で以下を満たす Ƹ 𝑝 を poly(𝑛, 𝑁) 時間で計算可能：
8. 8 Dinitz et al. ’21 を L 凸関数最小化に拡張し改良 計算量上界 𝑂(𝑇init + 𝑇loc × min 𝑝∗:optimal 𝑝∗ − Ƹ 𝑝 ∞ ) 二部マッチングの場合 𝑇init ≤ 𝑇loc = 𝑂(𝑚 𝑛) ①．問題クラスの拡張 ③．サンプル複雑度改善（𝑛3 → 𝑛） ④．𝑝∗ は Ƹ 𝑝 の最近最適解（タイブレイク不要） ②．ℓ1 から ℓ∞ に 学習可能性 L 凸関数 𝑔 が分布 𝒟 から i.i.d. で 𝑁 = ෩ Ω 𝑛𝐶2/𝜀2 個与えられれば， 高確率で以下を満たす Ƹ 𝑝 を poly(𝑛, 𝑁) 時間で計算可能： ４つの改善 𝔼𝑔∼𝒟 min 𝑝∗:optimal for 𝑔 𝑝∗ − Ƹ 𝑝 ∞ ≾ min ො 𝑝∗∈ −𝐶,+𝐶 𝑛 𝔼𝑔∼𝒟 min 𝑝∗:optimal for 𝑔 𝑝∗ − Ƹ 𝑝∗ ∞ + 𝜀 本研究：L 凸関数最小化への拡張と改良 予測 Ƹ 𝑝 から初期解を得る時間 一反復の時間（局所的に重み無し問題を解く時間） ① ② ③ ④
9. 𝑝 𝑞 𝑝 ∨ 𝑞 𝑝 ∧ 𝑞 9 関数 𝑔: ℤ𝑉 → ℝ ∪ {+∞} が L 凸 ⇔ 以下を満たす： 劣モジュラ性．𝑔 𝑝 + 𝑔 𝑞 ≥ 𝑔 𝑝 ∨ 𝑞 + 𝑔(𝑝 ∧ 𝑞) 𝟏 方向線形性．∃𝑟 ∈ ℝ，𝑔 𝑝 + 𝟏 = 𝑔 𝑝 + 𝑟 Murota. Discrete Convex Analysis. Discrete Mathematics and Applications. SIAM, 2003. Murota ’03 L 凸関数・L 凸集合
10. 𝑝 𝑞 𝑝 ∨ 𝑞 𝑝 ∧ 𝑞 10 重み付きマトロイド交叉，最小凸費用テンション問題なども扱える 関数 𝑔: ℤ𝑉 → ℝ ∪ {+∞} が L 凸 ⇔ 以下を満たす： 劣モジュラ性．𝑔 𝑝 + 𝑔 𝑞 ≥ 𝑔 𝑝 ∨ 𝑞 + 𝑔(𝑝 ∧ 𝑞) 𝟏 方向線形性．∃𝑟 ∈ ℝ，𝑔 𝑝 + 𝟏 = 𝑔 𝑝 + 𝑟 𝑔 𝑝 = σ 𝑖∈𝐿 𝑠𝑖 − σ 𝑗∈𝑅 𝑡𝑗 + σ 𝑖𝑗∈𝐸 𝛿𝑤𝑖𝑗 (𝑠𝑖 − 𝑡𝑗 ) 線形, 並進不変 𝛿𝑤 𝑥 = 0 if 𝑥 ≥ 𝑤 and +∞ otherwise. 一般に，凸関数 𝜓𝑖𝑗 : ℝ → ℝ ∪ {+∞} に対し 𝜓𝑖𝑗 (𝑝𝑖 − 𝑝𝑗 ) は L 凸 Murota. Discrete Convex Analysis. Discrete Mathematics and Applications. SIAM, 2003. 二部マッチングの双対 LP min σ 𝑖∈𝐿 𝑠𝑖 − σ 𝑗∈𝑅 𝑡𝑗 s. t. 𝑠𝑖 − 𝑡𝑗 ≥ 𝑤𝑖𝑗 (𝑖, 𝑗) ∈ 𝐸 Murota ’03 L 凸関数・L 凸集合 ※ L 凸 + L 凸 = L 凸
11. 1: 𝑝 ← 𝑝∘ 2: while not converged 3: 𝑑 ← argmin 𝑔 𝑝 + 𝑑′ 𝑑′ ∈ 0, +1 𝑉 ∪ 0, −1 𝑉 4: if 𝑔 𝑝 + 𝑑 − 𝑔 𝑝 = 0 5: return 𝑝 6: 𝜆 ← 1 (or use the long-step rule) 7: 𝑝 ← 𝑝 + 𝜆𝑑 11 重み無し問題 定理 (Murota and Shioura ’14) 最急降下法は 𝜇(𝑝∘) + 1 反復で停止し最適解を出力． 𝜇 𝑞 ≔ min 𝑝∗ − 𝑞 ∞ ± 𝑝∗ ∈ argmin 𝑔 (𝑞 ∈ ℤ𝑉, 𝑥 ∞ ± ≔ +𝑥 ∨ 0 ∞ + −𝑥 ∨ 0 ∞) Murota and Shioura. Exact bounds for steepest descent algorithms of L-convex function minimization. Oper. Res. Lett., 2014. 𝑝∘ 二部マッチングの場合はハンガリー法の双対版に対応 最急降下法
12. 予測ありの計算量上界 12 計算量上界 任意の予測 Ƹ 𝑝 ∈ ℝ𝑉 に対し ℓ∞ ± -射影 ො 𝑞 を以下で定義： ො 𝑞 ∈ argmin 𝑞 − Ƹ 𝑝 ∞ ± 𝑞 ∈ conv(dom 𝑔) 𝜇 の連続版 ҧ 𝜇: ℝ𝑉 → ℝ を定義： ℓ∞ ± -射影の時間を 𝑇init ，最急降下方向の計算時間 (Step 3) を 𝑇loc とすると ① 𝑝∘ = ො 𝑞 ∈ dom 𝑔 が成立（射影の四捨五入は実行可能解） ② 𝑝∘ を初期点とする最急降下法の反復回数は高々 2 ҧ 𝜇 Ƹ 𝑝 + 2 ҧ 𝜇 𝑞 ≔ min 𝑝∗ − 𝑞 ∞ ± 𝑝∗ ∈ conv(argmin 𝑔) (≤ 𝜇 𝑞 ) 総計算量 = 𝑂 𝑇init + 𝑇loc × ҧ 𝜇 Ƹ 𝑝 ≾ 𝑂(𝑇init + 𝑇loc × min 𝑝∗:optimal 𝑝∗ − Ƹ 𝑝 ∞ )
13. 13 以下全てのケースで 𝑇init ≤ 𝑇loc ，総計算量は 𝑂(𝑇loc × ҧ 𝜇 Ƹ 𝑝 ) 重み付き二部マッチング • Hopcroft−Karp の最大マッチング（最小被覆）アルゴリズムで 𝑇loc = 𝑂(𝑚 𝑛) 重み付きマトロイド交叉 • Cunningham の重み無しマトロイド交叉アルゴリズム で 𝑇loc = 𝑂(𝑛𝑟1.5 ⋅ 𝑇ind ) • Chakrabarty et al. ’19 のアルゴリズムで 𝑇loc = 𝑂(𝑛𝑟 log 𝑟 ⋅ 𝑇ind ) 最小凸費用テンション問題 • Dinic の最大流（最小カット）アルゴリズムで 𝑇loc = 𝑂 𝑚𝑛2 • Brand et al. ’23 のアルゴリズムで 𝑇loc = 𝑚1+𝑜(1) 𝑇loc = 重み無しの問題を解く計算時間 様々な場合の計算量 Brand et al. A deterministic almost-linear time algorithm for minimum-cost flow. FOCS, 2023 Chakrabarty et al. Faster matroid intersection. FOCS, 2019.
14. 14 ҧ 𝜇 Ƹ 𝑝; 𝑔 ≔ min 𝑝∗ − Ƹ 𝑝 ∞ ± 𝑝∗ ∈ conv(argmin 𝑔) 学習可能性 𝑉 上の L 凸関数（インスタンス）が分布 𝒟 から i.i.d. で 𝑁 = ෩ Ω 𝑛 𝐶/𝜀 2 個 与えられれば，高確率で以下を満たす Ƹ 𝑝 を poly(𝑛, 𝑁) 時間で計算可能： 𝑁 十分大で Ƹ 𝑝 を学習した場合 期待計算量 ≈ 𝑂(𝑇loc × 𝓓上最良の ҧ 𝜇 予測誤差)． 𝓓上最良の ҧ 𝜇 予測誤差（タイブレイク不要） 𝔼𝑔∼𝒟 ҧ 𝜇 Ƹ 𝑝; 𝑔 ≤ min ො 𝑝∗∈ −𝐶,+𝐶 𝑛 𝔼𝑔∼𝒟 ҧ 𝜇( Ƹ 𝑝∗; 𝑔) + 𝜀 予測の学習
15. 15 Cutkosky. Anytime online-to-batch, optimism and acceleration. ICML, 2019 Cesa-Bianchi et al. On the generalization ability of on-line learning algorithms. IEEE Trans. Inf. Theory, 2004. Zinkevich ’03 Zinkevich. Online convex programming and generalized infinitesimal gradient ascent. ICML, 2003. 学習方法：オンライン劣勾配降下法 𝑓𝑡 Ƹ 𝑝 ≔ ҧ 𝜇 Ƹ 𝑝; 𝑔𝑡 = min 𝑝∗ − Ƹ 𝑝 ∞ ± 𝑝∗ ∈ conv(argmin 𝑔𝑡 ) は凸関数． 1: Fix Ƹ 𝑝1 ∈ 𝐾 2: For 𝑡 = 1, … , 𝑁 3: Observe 𝑓𝑡 4: Compute 𝑧𝑡 ∈ 𝜕𝑓𝑡 Ƹ 𝑝𝑡 5: Ƹ 𝑝𝑡+1 ← Proj −𝐶,+𝐶 𝑛 Ƹ 𝑝𝑡 − 𝜂𝑧𝑡 6: Return 1 𝑁 σ𝑡=1 𝑁 Ƹ 𝑝𝑡 劣勾配計算： 相補性＋Danskin の定理 リグレット上界 ＋ online-to-batch により学習可能性の保証が得られる．
16. 16 10 頂点の二部マッチングのインスタンスをランダムに 𝑡 = 1,2, … , 1000 まで生成 訓練データ数 予測 Ƹ 𝑝 の計算方法 4 種を比較： • OSD + 損失関数 3 種： ҧ 𝜇( Ƹ 𝑝; 𝑔𝑡 ) (提案法), ℓ1 -誤差 (Dinitz et al. ’21), ℓ∞ -誤差 • Cold-start 訓練データの増加に対するハンガリー法の反復回数の変化 実験 訓練データ数 枝重みノイズ弱 枝重みノイズ強
17. 17 L 凸関数最小化に対する最急降下法の計算量は 予測 Ƹ 𝑝 から最適解集合への ℓ∞ ± 距離 ҧ 𝜇 Ƹ 𝑝; 𝑔 に比例 ҧ 𝜇 Ƹ 𝑝; 𝑔 は Ƹ 𝑝 ∈ ℝ𝑉 について凸関数であり，オンライン劣勾配法によって 観測された 𝑔1 , 𝑔2 , … の最適解集合に近くなるように学習可能 M 凸関数最小化についても類似の結果が得られる (Oki and S ’23) （タイブレイクの扱いについては未解決） Oki and S. Faster discrete convex function minimization with predictions: The M-convex case, NeurIPS, 2023. まとめ 今後の展開 • L 凸，M 凸以外の問題への適用（特に一般マッチング） • 平均計算時間以外の指標（CVaR, 一定時間内に解ける問題数, など） • 予測モデルの学習（新たなインスタンスに対し適切な予測を生成）
18. Supplementary
19. 19 1 2 3 4 5 6 7 8 4 2 3 3 2 2 4 4 4 4 0 0 0 0 1 1 1 1 𝐿 𝑅 初期実行可能解 𝑝 = 𝑠, 𝑡 ∈ ℤ𝐿∪𝑅 を設定 双対 LP（変数 𝑝 = 𝑠, 𝑡 ∈ ℝ𝐿∪𝑅 ） min σ 𝑖∈𝐿 𝑠𝑖 − σ 𝑗∈𝑅 𝑡𝑗 s. t. 𝑠𝑖 − 𝑡𝑗 ≥ 𝑤𝑖𝑗 (𝑖, 𝑗) ∈ 𝐸 Schrijver. Combinatorial Optimization: Polyhedra and Efficiency. Springer, 2003. cf. Schrijver ’03 (Section 18.5b) ハンガリー法（双対版）
20. 20 1 2 3 4 5 6 7 8 4 2 3 3 2 2 4 4 4 4 0 0 0 0 1 1 1 1 𝐿 𝑅 双対 LP（変数 𝑝 = 𝑠, 𝑡 ∈ ℝ𝐿∪𝑅 ） min σ 𝑖∈𝐿 𝑠𝑖 − σ 𝑗∈𝑅 𝑡𝑗 s. t. 𝑠𝑖 − 𝑡𝑗 ≥ 𝑤𝑖𝑗 (𝑖, 𝑗) ∈ 𝐸 Schrijver. Combinatorial Optimization: Polyhedra and Efficiency. Springer, 2003. cf. Schrijver ’03 (Section 18.5b) ハンガリー法（双対版） 初期実行可能解 𝑝 = 𝑠, 𝑡 ∈ ℤ𝐿∪𝑅 を設定 タイト部分グラフ (𝑉, 𝐸∗) を見つける
21. 21 1 2 3 4 5 6 7 8 4 2 3 3 2 2 4 4 4 4 0 0 0 0 𝑆 1 1 1 1 (𝑉, 𝐸∗) の最小頂点被覆 (𝑆, 𝑇) を見つける 𝐿 𝑅 双対 LP（変数 𝑝 = 𝑠, 𝑡 ∈ ℝ𝐿∪𝑅 ） min σ 𝑖∈𝐿 𝑠𝑖 − σ 𝑗∈𝑅 𝑡𝑗 s. t. 𝑠𝑖 − 𝑡𝑗 ≥ 𝑤𝑖𝑗 (𝑖, 𝑗) ∈ 𝐸 Schrijver. Combinatorial Optimization: Polyhedra and Efficiency. Springer, 2003. cf. Schrijver ’03 (Section 18.5b) ハンガリー法（双対版） 4 初期実行可能解 𝑝 = 𝑠, 𝑡 ∈ ℤ𝐿∪𝑅 を設定 タイト部分グラフ (𝑉, 𝐸∗) を見つける
22. 22 1 2 3 4 5 6 7 8 4 2 3 3 2 2 4 4 5 4 1 1 1 1 1 1 1 1 𝑝 ← (𝑠 + 𝜆𝟏𝑆 , 𝑡 + 𝜆𝟏𝑅∖𝑇 ) 𝐿 𝑅 双対 LP（変数 𝑝 = 𝑠, 𝑡 ∈ ℝ𝐿∪𝑅 ） min σ 𝑖∈𝐿 𝑠𝑖 − σ 𝑗∈𝑅 𝑡𝑗 s. t. 𝑠𝑖 − 𝑡𝑗 ≥ 𝑤𝑖𝑗 (𝑖, 𝑗) ∈ 𝐸 Schrijver. Combinatorial Optimization: Polyhedra and Efficiency. Springer, 2003. cf. Schrijver ’03 (Section 18.5b) ハンガリー法（双対版） (𝑉, 𝐸∗) の最小頂点被覆 (𝑆, 𝑇) を見つける 初期実行可能解 𝑝 = 𝑠, 𝑡 ∈ ℤ𝐿∪𝑅 を設定 タイト部分グラフ (𝑉, 𝐸∗) を見つける
23. 23 1 2 3 4 5 6 7 8 4 2 3 3 2 2 4 4 5 4 1 1 1 1 1 1 1 1 𝐿 𝑅 双対 LP（変数 𝑝 = 𝑠, 𝑡 ∈ ℝ𝐿∪𝑅 ） min σ 𝑖∈𝐿 𝑠𝑖 − σ 𝑗∈𝑅 𝑡𝑗 s. t. 𝑠𝑖 − 𝑡𝑗 ≥ 𝑤𝑖𝑗 (𝑖, 𝑗) ∈ 𝐸 Schrijver. Combinatorial Optimization: Polyhedra and Efficiency. Springer, 2003. cf. Schrijver ’03 (Section 18.5b) ハンガリー法（双対版） 𝑝 ← (𝑠 + 𝜆𝟏𝑆 , 𝑡 + 𝜆𝟏𝑅∖𝑇 ) (𝑉, 𝐸∗) の最小頂点被覆 (𝑆, 𝑇) を見つける 初期実行可能解 𝑝 = 𝑠, 𝑡 ∈ ℤ𝐿∪𝑅 を設定 タイト部分グラフ (𝑉, 𝐸∗) を見つける 上記を繰り返し |𝐸∗| を増やす ⋮
24. 24 1 2 3 4 5 6 7 8 4 2 3 3 2 2 4 4 5 5 2 1 2 2 𝑆 𝑇 1 1 1 1 𝐿 𝑅 双対 LP（変数 𝑝 = 𝑠, 𝑡 ∈ ℝ𝐿∪𝑅 ） min σ 𝑖∈𝐿 𝑠𝑖 − σ 𝑗∈𝑅 𝑡𝑗 s. t. 𝑠𝑖 − 𝑡𝑗 ≥ 𝑤𝑖𝑗 (𝑖, 𝑗) ∈ 𝐸 Schrijver. Combinatorial Optimization: Polyhedra and Efficiency. Springer, 2003. cf. Schrijver ’03 (Section 18.5b) ハンガリー法（双対版） 上記を繰り返し |𝐸∗| を増やす ⋮ 𝑝 ← (𝑠 + 𝜆𝟏𝑆 , 𝑡 + 𝜆𝟏𝑅∖𝑇 ) (𝑉, 𝐸∗) の最小頂点被覆 (𝑆, 𝑇) を見つける 初期実行可能解 𝑝 = 𝑠, 𝑡 ∈ ℤ𝐿∪𝑅 を設定 タイト部分グラフ (𝑉, 𝐸∗) を見つける
25. 25 1 2 3 4 5 6 7 8 4 2 3 3 2 2 4 4 5 5 2 1 2 2 1 1 1 1 𝐿 𝑅 If 𝐸∗ に完全マッチングあり（ 𝑆 + 𝑇 = 𝑛/2） 双対 LP（変数 𝑝 = 𝑠, 𝑡 ∈ ℝ𝐿∪𝑅 ） min σ 𝑖∈𝐿 𝑠𝑖 − σ 𝑗∈𝑅 𝑡𝑗 s. t. 𝑠𝑖 − 𝑡𝑗 ≥ 𝑤𝑖𝑗 (𝑖, 𝑗) ∈ 𝐸 Schrijver. Combinatorial Optimization: Polyhedra and Efficiency. Springer, 2003. cf. Schrijver ’03 (Section 18.5b) ハンガリー法（双対版） 上記を繰り返し |𝐸∗| を増やす ⋮ 𝑝 ← (𝑠 + 𝜆𝟏𝑆 , 𝑡 + 𝜆𝟏𝑅∖𝑇 ) (𝑉, 𝐸∗) の最小頂点被覆 (𝑆, 𝑇) を見つける 初期実行可能解 𝑝 = 𝑠, 𝑡 ∈ ℤ𝐿∪𝑅 を設定 タイト部分グラフ (𝑉, 𝐸∗) を見つける
26. 26 1 2 3 4 5 6 7 8 4 2 3 3 2 2 4 4 5 5 2 1 2 2 1 1 1 1 𝐿 𝑅 双対 LP（変数 𝑝 = 𝑠, 𝑡 ∈ ℝ𝐿∪𝑅 ） min σ 𝑖∈𝐿 𝑠𝑖 − σ 𝑗∈𝑅 𝑡𝑗 s. t. 𝑠𝑖 − 𝑡𝑗 ≥ 𝑤𝑖𝑗 (𝑖, 𝑗) ∈ 𝐸 Schrijver. Combinatorial Optimization: Polyhedra and Efficiency. Springer, 2003. cf. Schrijver ’03 (Section 18.5b) ハンガリー法（双対版） If 𝐸∗ に完全マッチングあり（ 𝑆 + 𝑇 = 𝑛/2） 𝐸∗の完全マッチングを一つ出力 上記を繰り返し |𝐸∗| を増やす ⋮ 𝑝 ← (𝑠 + 𝜆𝟏𝑆 , 𝑡 + 𝜆𝟏𝑅∖𝑇 ) (𝑉, 𝐸∗) の最小頂点被覆 (𝑆, 𝑇) を見つける 初期実行可能解 𝑝 = 𝑠, 𝑡 ∈ ℤ𝐿∪𝑅 を設定 タイト部分グラフ (𝑉, 𝐸∗) を見つける
27. 27 𝑂(𝑚 𝑛) 時間 （Hopcroft−Karp の重み無しアルゴリズム） 総計算量（初期解 𝑝 は所与として） 𝑂(𝑚 𝑛 × 反復回数) 双対 LP（変数 𝑝 = 𝑠, 𝑡 ∈ ℝ𝐿∪𝑅 ） min σ 𝑖∈𝐿 𝑠𝑖 − σ 𝑗∈𝑅 𝑡𝑗 s. t. 𝑠𝑖 − 𝑡𝑗 ≥ 𝑤𝑖𝑗 (𝑖, 𝑗) ∈ 𝐸 Schrijver. Combinatorial Optimization: Polyhedra and Efficiency. Springer, 2003. cf. Schrijver ’03 (Section 18.5b) ハンガリー法（双対版） If 𝐸∗ に完全マッチングあり（ 𝑆 + 𝑇 = 𝑛/2） 𝐸∗の完全マッチングを一つ出力 上記を繰り返し |𝐸∗| を増やす ⋮ 𝑝 ← (𝑠 + 𝜆𝟏𝑆 , 𝑡 + 𝜆𝟏𝑅∖𝑇 ) (𝑉, 𝐸∗) の最小頂点被覆 (𝑆, 𝑇) を見つける 初期実行可能解 𝑝 = 𝑠, 𝑡 ∈ ℤ𝐿∪𝑅 を設定 タイト部分グラフ (𝑉, 𝐸∗) を見つける
28. 28 𝑂(𝑚 𝑛) 時間 （Hopcroft−Karp の重み無しアルゴリズム） 総計算量（初期解 𝑝 は所与として） 𝑂(𝑚 𝑛 × 反復回数) 双対 LP（変数 𝑝 = 𝑠, 𝑡 ∈ ℝ𝐿∪𝑅 ） min σ 𝑖∈𝐿 𝑠𝑖 − σ 𝑗∈𝑅 𝑡𝑗 s. t. 𝑠𝑖 − 𝑡𝑗 ≥ 𝑤𝑖𝑗 (𝑖, 𝑗) ∈ 𝐸 Schrijver. Combinatorial Optimization: Polyhedra and Efficiency. Springer, 2003. cf. Schrijver ’03 (Section 18.5b) ハンガリー法（双対版） If 𝐸∗ に完全マッチングあり（ 𝑆 + 𝑇 = 𝑛/2） 𝐸∗の完全マッチングを一つ出力 上記を繰り返し |𝐸∗| を増やす ⋮ 𝑝 ← (𝑠 + 𝜆𝟏𝑆 , 𝑡 + 𝜆𝟏𝑅∖𝑇 ) (𝑉, 𝐸∗) の最小頂点被覆 (𝑆, 𝑇) を見つける 初期実行可能解 𝑝 = 𝑠, 𝑡 ∈ ℤ𝐿∪𝑅 を設定 タイト部分グラフ (𝑉, 𝐸∗) を見つける 初期解 𝑝 を上手く選べば高速化可能？
29. 𝑝 𝑞 𝑝 ∨ 𝑞 𝑝 ∧ 𝑞 29 関数 𝑔: ℤ𝑉 → ℝ ∪ {+∞} が L 凸 ⇔ 以下を満たす： 劣モジュラ性．𝑔 𝑝 + 𝑔 𝑞 ≥ 𝑔 𝑝 ∨ 𝑞 + 𝑔(𝑝 ∧ 𝑞) 𝟏 方向線形性．∃𝑟 ∈ ℝ，𝑔 𝑝 + 𝟏 = 𝑔 𝑝 + 𝑟 重要な性質 • L 凸 + L 凸 = L 凸 • 𝑔 が L 凸関数 ⇒ dom 𝑔 ≔ 𝑝 ∈ ℤ𝑉 𝑔 𝑝 < +∞}，及び argmin 𝑔 は L 凸集合 • L 凸集合の不等式系表現：𝑆 = 𝑝 ∈ ℤ𝑉 | 𝑝𝑗 − 𝑝𝑖 ≤ 𝛾𝑖𝑗 ∀𝑖, 𝑗 ∈ 𝑉; 𝑖 ≠ 𝑗 for some 𝛾𝑖𝑗 ∈ ℤ ∪ {+∞} • ℤ𝑉 を ℝ𝑉 に置き換えると 𝑆 の凸包 集合 𝑆 ⊆ ℤ𝑉 が L 凸 ⇔ 以下を満たす： 𝑝, 𝑞 ∈ 𝑆 ⇒ 𝑝 ∨ 𝑞, 𝑝 ∧ 𝑞 ∈ 𝑆 and 𝑝 ∈ 𝑆 ⇒ 𝑝 ± 𝟏 ∈ 𝑆 Murota. Discrete Convex Analysis. Discrete Mathematics and Applications. SIAM, 2003. Murota ’03 L 凸関数・L 凸集合
30. 予測ありの計算量上界 30 計算量上界 任意の予測 Ƹ 𝑝 ∈ ℝ𝑉 に対し ℓ∞ ± -射影 ො 𝑞 を以下で定義： ො 𝑞 ∈ argmin 𝑞 − Ƹ 𝑝 ∞ ± 𝑞 ∈ conv(dom 𝑔) 𝜇 の連続版 ҧ 𝜇: ℝ𝑉 → ℝ を定義： ℓ∞ ± -射影の時間を 𝑇init , 最急降下方向の計算時間 (Step 3) を 𝑇loc とすると 総計算量 𝑂(𝑇init + 𝑇loc × ҧ 𝜇 Ƹ 𝑝 ) （任意の 𝑝∗ ∈ argmin 𝑔 に対し ҧ 𝜇 Ƹ 𝑝 ≤ 𝜇 Ƹ 𝑝 ≤ 𝑝∗ − Ƹ 𝑝 ∞ ± ≤ 2 𝑝∗ − Ƹ 𝑝 ∞ ） ❶ 𝑝∘ = ො 𝑞 ∈ dom 𝑔 が成立（射影の四捨五入は実行可能解） ❷ 𝑝∘ を初期点とする最急降下法の反復回数は高々 2 ҧ 𝜇 Ƹ 𝑝 + 2 ҧ 𝜇 𝑞 ≔ min 𝑝∗ − 𝑞 ∞ ± 𝑝∗ ∈ conv(argmin 𝑔) (≤ 𝜇 𝑞 )
31. 31 𝑞 ∈ conv dom 𝑔 ⇒ 𝑞 ∈ dom 𝑔 を示せば十分． L 凸集合の不等式表現： • dom 𝑔 = 𝑝 ∈ ℤ𝑉 | 𝑝𝑗 − 𝑝𝑖 ≤ 𝛾𝑖𝑗 ∀𝑖, 𝑗 ∈ 𝑉; 𝑖 ≠ 𝑗 （𝛾𝑖𝑗 ∈ ℤ ∪ {+∞}） • conv(dom 𝑔) = 𝑝 ∈ ℝ𝑉 | 𝑝𝑗 − 𝑝𝑖 ≤ 𝛾𝑖𝑗 ∀𝑖, 𝑗 ∈ 𝑉; 𝑖 ≠ 𝑗 （𝛾𝑖𝑗 ∈ ℤ ∪ {+∞}） 𝑞 ∈ conv dom 𝑔 より 𝑞𝑗 − 𝑞𝑖 ≤ 𝛾𝑖𝑗 ． ⌊𝑞𝑗 ⌉ − 𝑞𝑖 ≤ ⌈𝑞𝑗 −𝑞𝑖 ⌉ ≤ 𝛾𝑖𝑗 より 𝑞 ∈ dom 𝑔． さらに 𝛾𝑖𝑗 ∈ ℤ ∪ {+∞} より ⌈𝑞𝑗 −𝑞𝑖 ⌉ ≤ 𝛾𝑖𝑗 ． ① の証明
32. 𝑝∗ − 𝑝∘ ∞ ± ≤ 𝑝∗ − ො 𝑞 ∞ ± + 1 ≤ 𝑝∗ − Ƹ 𝑝 ∞ ± + Ƹ 𝑝 − ො 𝑞 ∞ ± + 1 ≤ 2 ҧ 𝜇 Ƹ 𝑝 + 1 32 (Murota and Shioura ’14) より反復回数は 𝜇 𝑝∘ + 1． 整数点上では 𝜇 = ҧ 𝜇 だから 𝜇 𝑝∘ = ҧ 𝜇 𝑝∘ ≤ 𝑝∗ − 𝑝∘ ∞ ± ．右辺を上から抑える： 𝜇 𝑝∘ ≤ 2 ҧ 𝜇 Ƹ 𝑝 + 1 を示す（𝑝∘ = ො 𝑞 , ො 𝑞 ∈ argmin 𝑞 − Ƹ 𝑝 ∞ ± 𝑞 ∈ conv(dom 𝑔) ）． 𝑝∗ ∈ conv(argmin 𝑔) を ҧ 𝜇 Ƹ 𝑝 = 𝑝∗ − Ƹ 𝑝 ∞ ± を達成するように定義． 𝑝∗ − Ƹ 𝑝 ∞ ± = ҧ 𝜇 Ƹ 𝑝 conv dom 𝑔 ⊇ conv(argmin 𝑔) より Ƹ 𝑝 − ො 𝑞 ∞ ± ≤ Ƹ 𝑝 − 𝑝∗ ∞ ± = ҧ 𝜇 Ƹ 𝑝 𝑝∘ = ො 𝑞 三角不等式 ② の証明
33. 33 敵対者：凸関数 𝑓𝑡 : ℝ𝑛 → ℝ (𝑡 = 1, … , 𝑇) を任意に決定 プレイヤー：𝑓1 , … , 𝑓𝑡 を観測し 𝑥𝑡+1 を右手順で計算 （𝐾 ⊆ ℝ𝑛 は ℓ2-直径 𝐷 の凸集合） 1: Fix 𝑥1 ∈ 𝐾 2: For 𝑡 = 1, … , 𝑇 3: Observe 𝑓𝑡 4: Compute 𝑧𝑡 ∈ 𝜕𝑓𝑡 𝑥𝑡 5: 𝑥𝑡+1 ← Proj𝐾 (𝑥𝑡 − 𝜂𝑧𝑡 ) Cutkosky. Anytime online-to-batch, optimism and acceleration. ICML, 2019 Cesa-Bianchi et al. On the generalization ability of on-line learning algorithms. IEEE Trans. Inf. Theory, 2004. Zinkevich ’03 Zinkevich. Online convex programming and generalized infinitesimal gradient ascent. ICML, 2003. リグレット上界 max 𝑡 𝑧𝑡 2 ≤ 𝐿 ならば，𝜂 = 𝐷/(𝐿 𝑇) とすると σ𝑡=1 𝑇 𝑓𝑡 𝑥𝑡 ≤ min 𝑥∗∈𝐾 σ𝑡=1 𝑇 𝑓𝑡 𝑥∗ + 𝑂(𝐷𝐿 𝑇) Online-to-Batch (Cesa-Bianchi et al. ’04, Cutkosky ’19) 𝑓𝑡 が i.i.d., max 𝑡 𝑧𝑡 1 ≤ 𝐿1 , 𝐾 の ℓ∞ -直径が 𝐷∞ なら，確率 1 − 𝛿 で 𝔼 𝑓 1 𝑇 σ𝑡 𝑥𝑡 ≤ min 𝑥∗∈𝐾 𝔼 𝑓(𝑥∗) + 𝑂 𝐷𝐿 𝑛+𝐷∞𝐿1 log(1/𝛿) 𝑇 オンライン劣勾配降下法（OSD）
34. 学習可能性 34 時刻 𝑡 の損失関数を ҧ 𝜇 ⋅; 𝑔𝑡 ≔ min 𝑝∗ −⋅ ∞ ± 𝑝∗ ∈ conv(argmin 𝑔𝑡 ) とする • ҧ 𝜇 ⋅; 𝑔𝑡 は ℝ𝑉 上の凸関数 • 劣勾配 𝑧𝑡 ∈ 𝜕 ҧ 𝜇 ⋅; 𝑔𝑡 s.t. 𝑧𝑡 1 = 𝑂(1) を多項式時間で計算可（後述） • 予測のドメイン 𝐾 = −𝐶, +𝐶 𝑛 の ℓ2 -直径は 2𝐶 𝑛，ℓ∞ -直径は 2𝐶 ҧ 𝜇 ⋅; 𝑔1 , … , ҧ 𝜇 ⋅; 𝑔𝑇 に対する OSD の出力を Ƹ 𝑝𝑡 ， Ƹ 𝑝 = 1 𝑇 σ𝑡 Ƹ 𝑝𝑡 とすると 𝑇 = Ω 𝐶/𝜀 2(𝑛 + log(1/𝛿)) とすれば，誤差は高々 𝜀． 𝔼𝑔∼𝒟 ҧ 𝜇 Ƹ 𝑝; 𝑔 ≤ min ො 𝑝∗∈ −𝐶,+𝐶 𝑛 𝔼𝑔∼𝒟 ҧ 𝜇 Ƹ 𝑝∗; 𝑔 + 𝑂 𝐶(𝑛+ log(1/𝛿)) 𝑇 . 𝑔𝑡 ∼ 𝒟 は L 凸関数 ҧ 𝜇 についての理論保証
35. 35 Proj𝐾 で 𝑥∗ ∈ 𝐾 に近づくから 𝑥𝑡+1 − 𝑥∗ 2 = Proj𝐾 (𝑥𝑡 − 𝜂𝑧𝑡 ) − 𝑥∗ 2 ≤ 𝑥𝑡 − 𝜂𝑧𝑡 − 𝑥∗ 2 = 𝑥𝑡 − 𝑥∗ 2 + 𝜂2 𝑧𝑡 2 + 2𝜂 𝑧𝑡 , 𝑥𝑡 − 𝑥∗ . 凸性と上式より 𝑓𝑡 𝑥𝑡 − 𝑓𝑡 𝑥∗ ≤ − 𝑧𝑡 , 𝑥𝑡 − 𝑥∗ ≤ 𝑥𝑡−𝑥∗ 2− 𝑥𝑡+1−𝑥∗ 2 2𝜂 + 𝜂 2 𝑧𝑡 2. 𝑡 について和を取り，− 𝑥𝑇+1 − 𝑥∗ 2≤ 0 を無視して σ𝑡 𝑓𝑡 𝑥𝑡 − 𝑓𝑡 𝑥∗ ≤ 𝑥1−𝑥∗ 2 2𝜂 + 𝜂 2 σ𝑡 𝑧𝑡 2 ≤ 𝐷2 2𝜂 + 𝜂𝐿2 2 𝑇 = 1 2 𝐷𝐿 𝑇 𝜂 = 𝐷/(𝐿 𝑇) OSD のリグレット上界の証明 1: Fix 𝑥1 ∈ 𝐾 2: For 𝑡 = 1, … , 𝑇 3: Observe 𝑓𝑡 4: Compute 𝑧𝑡 ∈ 𝜕𝑓𝑡 𝑥𝑡 5: 𝑥𝑡+1 ← Proj𝐾 (𝑥𝑡 − 𝜂𝑧𝑡 )
36. 36 ҧ 𝜇 ⋅; 𝑔𝑡 = min 𝑝∗ −⋅ ∞ ± 𝑝∗ ∈ conv(argmin 𝑔𝑡 ) の劣勾配を以下の手順で計算 Step 2: Step 1 の不等式系から，以下を満たすグラフ 𝐺 = (𝑉, 𝐸, {෥ 𝑤𝑖𝑗 Ƹ 𝑝 }) を構成 ҧ 𝜇 Ƹ 𝑝; 𝑔𝑡 = “𝐺 の最短路長” Step 3: Danskin の定理 • 𝜙 が Ƹ 𝑝 について凸 • max がコンパクト集合上 ならば適用可 Step 1: conv(argmin 𝑔𝑡 ) の不等式系表現を構成 = − max 𝑃:simple path {𝜙 Ƹ 𝑝; 𝑃 ≔ − σ 𝑖𝑗∈𝑃 ෥ 𝑤𝑖𝑗 ( Ƹ 𝑝)} Step 1 Step 2 Step 3 𝜕max 𝜙 Ƹ 𝑝; 𝑃 𝑃: simple path = conv ∇𝜙 Ƹ 𝑝; 𝑃∗ 𝑃∗: shortest path for Ƹ 𝑝 (𝑉, 𝐸, {෥ 𝑤𝑖𝑗 Ƹ 𝑝 }) 上の最短路 𝑃∗ を見つけて ∇𝜙 Ƹ 𝑝; 𝑃∗ = σ𝑖𝑗∈𝑃∗ ∇෥ 𝑤𝑖𝑗 ( Ƹ 𝑝) を計算すれば OK 劣勾配の計算：全体の流れ
37. 37 Murota. Combinatorial relaxation algorithm for the maximum degree of subdeterminants: Computing Smith-Mcmillan form at infinity and structural indices in Kronecker form. Appl. Algebra Eng. Commun. Comput., 1995. 最大重み完全二部マッチングの場合で説明 双対 LP（変数 𝑝 = 𝑠, 𝑡 ∈ ℝ𝐿∪𝑅 ） min σ 𝑖∈𝐿 𝑠𝑖 − σ 𝑗∈𝑅 𝑡𝑗 s. t. 𝑠𝑖 − 𝑡𝑗 ≥ 𝑤𝑖𝑗 (𝑖, 𝑗) ∈ 𝐸 主 LP（変数 𝑥 ∈ ℝ𝐸 ） max σ𝑒∈𝐸 𝑤𝑒 𝑥𝑒 s. t. σ 𝑒∈𝑁(𝑣) 𝑥𝑒 = 1 𝑣 ∈ 𝐿 ∪ 𝑅 , 𝑥 ≥ 0 主双対実行可能解の組 (𝑥, 𝑝) が最適 ⟺ 以下の相補性条件を満たす： 𝑥𝑒 𝑠𝑖 − 𝑡𝑗 − 𝑤𝑖𝑗 = 0 𝑒 = 𝑖, 𝑗 ∈ 𝐸 任意の主最適解 𝑥 を用いて conv(argmin 𝑔𝑡 ) を以下のように表現可： cf. Murota ’95 𝑝 = 𝑠, 𝑡 ∈ ℝ𝐿∪𝑅 𝑠𝑖 − 𝑡𝑗 ≥ 𝑤𝑖𝑗 for 𝑖, 𝑗 ∈ 𝐸 𝑠𝑖 − 𝑡𝑗 ≤ 𝑤𝑖𝑗 for 𝑖, 𝑗 ∈ 𝐸 s. t. 𝑥𝑖𝑗 > 0 conv(argmin 𝑔𝑡 ) の不等式系表現の構成
38. 38 ҧ 𝜇 Ƹ 𝑝; 𝑔𝑡 = min 𝑝∗ − Ƹ 𝑝 ∞ ± 𝑝∗ ∈ conv(argmin 𝑔𝑡 ) = 𝐺 上の最短路長 𝑖 𝑗 • 𝑉：台集合 • 𝐸：全有向辺 𝑖, 𝑗 の集合 • ෥ 𝑤𝑖𝑗 ( Ƹ 𝑝)：𝛾𝑖𝑗 − Ƹ 𝑝𝑗 + Ƹ 𝑝𝑖 証明概略 左辺の LP 表現 ＝ (𝑉, 𝐸, ෥ 𝑤𝑖𝑗 ( Ƹ 𝑝)) 上の最短路問題の双対 LP で定義すると 一般に conv(argmin 𝑔𝑡 ) は L 凸集合の凸包であり， ある 𝛾𝑖𝑗 ∈ ℤ ∪ {+∞} を用いて差分不等式系表現可能： 𝑝 ∈ ℝ𝑉 𝑝𝑗 − 𝑝𝑖 ≤ 𝛾𝑖𝑗 (𝑖, 𝑗 ∈ 𝑉, 𝑖 ≠ 𝑗) グラフ 𝐺 = (𝑉, 𝐸, {෥ 𝑤𝑖𝑗 Ƹ 𝑝 }) を グラフの構成
39. 39 𝑃∗ の始点・終点以外の成分は打ち消し合い，結局劣勾配は 𝑧𝑡 = 0, … , −1, 0, … , +1, 0, … ∈ 𝜕 ҧ 𝜇 Ƹ 𝑝; 𝑔𝑡 (𝑉, 𝐸, ෥ 𝑤𝑖𝑗 ( Ƹ 𝑝)) 上の最短路 𝑃∗ を見つけて以下を計算： σ 𝑖𝑗∈𝑃∗ ∇෥ 𝑤𝑖𝑗 ( Ƹ 𝑝) = σ 𝑖𝑗∈𝑃∗ ∇(𝛾𝑖𝑗 − Ƹ 𝑝𝑗 + Ƹ 𝑝𝑖 ) 幾何的解釈 −𝑧𝑡 方向への移動は，conv(argmin 𝑔𝑡 ) の不等式表現 𝑝 ∈ ℝ𝑉 𝑝𝑗 − 𝑝𝑖 ≤ 𝛾𝑖𝑗 (𝑖, 𝑗 ∈ 𝑉, 𝑖 ≠ 𝑗) の中で違反度 Ƹ 𝑝𝑗 − Ƹ 𝑝𝑖 − 𝛾𝑖𝑗 > 0 が最も大きい 𝑖𝑗 を改善 𝑖th 𝑗th conv(argmin 𝑔𝑡 ) Ƹ 𝑝 𝑝𝑗 − 𝑝𝑖 ≤ 𝛾𝑖𝑗 𝑖 𝑗 具体的な計算方法