半正定値計画問題に対する 行列補完理論の高速実装

Presentation on theme: "半正定値計画問題に対する 行列補完理論の高速実装"— Presentation transcript:

1 半正定値計画問題に対する 行列補完理論の高速実装
東京工業大学　数理・計算科学専攻 　 山下 真 この研究は笹川研究助成の助成を受けています 2012/12/15 SOTA:「最適化の理論と応用」研究部会

2 行列補完の簡単な例 半正定値制約がある 半正定値計画問題(SemiDefinite Program, SDP)の例 2012/12/15
Notation 内積計算に必要なのは青の要素のみ これだけで計算できれば高速化できる 半正定値制約がある 2012/12/15 SOTA:「最適化の理論と応用」研究部会

3 SOTA:「最適化の理論と応用」研究部会
組合せ最適化のときにあらわれる疎構造 格子ネットワーク(p,q)上の最大クリーク問題 内積計算に必要な要素のみを青で表示 ほとんどがゼロ要素⇒こういった特徴を利用して高速化 半正定値条件を取り入れる⇒行列補完 2012/12/15 SOTA:「最適化の理論と応用」研究部会

4 SOTA:「最適化の理論と応用」研究部会
行列補完の簡単な例の続き 黒なしで計算 黒を適切に補完 2012/12/15 SOTA:「最適化の理論と応用」研究部会

5 SOTA:「最適化の理論と応用」研究部会
今日の話の方向性 大規模な半正定値計画問題に対して 行列補完を利用して高速に求解 行列補完自体は SDPA-C に実装済み 今回のメインは SDPA-C の改良にあたる 話の内容の多くは、行列補完について 2012/12/15 SOTA:「最適化の理論と応用」研究部会

6 SOTA:「最適化の理論と応用」研究部会
インデックス 半正定値計画問題の標準形と主双対内点法 行列補完理論と主双対内点法 実装の高速化 行列補完計算の改善 Schur complement マルチスレッド化 主探索方向マルチスレッド化 数値実験結果 2012/12/15 SOTA:「最適化の理論と応用」研究部会

7 SOTA:「最適化の理論と応用」研究部会
半正定値計画問題の主な応用 制御理論 組合せ最適化問題の緩和 ロバスト最適化 量子化学における電子構造計算 量子計算 (物理状態は半正定値行列で表現できるため) 2012/12/15 SOTA:「最適化の理論と応用」研究部会

8 半正定値計画問題の標準形 (SemiDefinite Programs, SDP)
主問題と双対問題で定義 行列補完理論で本質的なのは双対側 2012/12/15 SOTA:「最適化の理論と応用」研究部会

9 SOTA:「最適化の理論と応用」研究部会
主双対内点法の特徴 主問題と双対問題を同時に求解 多項式時間アルゴリズム 多くのソフトウェア SDPA, CSDP, SeDuMi, SDPT3, … 主問題の変数が密行列 Schur 補完行列がボトルネック 2012/12/15 SOTA:「最適化の理論と応用」研究部会

10 SOTA:「最適化の理論と応用」研究部会
主双対内点法の枠組み 実行可能領域 探索方向 ステップ長 中心パス 最適解 2012/12/15 SOTA:「最適化の理論と応用」研究部会

11 SOTA:「最適化の理論と応用」研究部会
中心パスの点と探索方向計算 最適性条件 中心パスの点 探索方向 の計算はニュートン法に基づく 2012/12/15 SOTA:「最適化の理論と応用」研究部会

12 SOTA:「最適化の理論と応用」研究部会
探索方向計算の続き 変数を左に集めて 変形していくと まとめると 2012/12/15 SOTA:「最適化の理論と応用」研究部会

13 SOTA:「最適化の理論と応用」研究部会
行列補完理論で 計算式が変わる部分 Schur 補完方程式 主問題側変数の計算 行列補完理論では密行列を避けて疎行列(下三角)で計算する 2012/12/15 SOTA:「最適化の理論と応用」研究部会

14 SOTA:「最適化の理論と応用」研究部会
インデックス 半正定値計画問題の標準形と主双対内点法 行列補完理論と主双対内点法 実装の高速化 行列補完計算の改善 Schur complement マルチスレッド化 主探索方向マルチスレッド化 数値実験結果 2012/12/15 SOTA:「最適化の理論と応用」研究部会

15 SOTA:「最適化の理論と応用」研究部会
行列補完理論と主双対内点法 変数行列の疎構造に基づいてコレスキー分解 Schur 補完行列を効率的に計算できる 詳しくは、以下の論文 Exploiting sparsity in semidefinite programming via matrix completion I: general framework Fukuda, Kojima, Murota, Nakata SIAM J. Optim, 2000 Exploiting sparsity in semidefinite programming via matrix completion II: implementation and numerical results Nakata, Fujisawa, Fukuda, Kojima, Murota Math. Program. B, 2003 2012/12/15 SOTA:「最適化の理論と応用」研究部会

16 SOTA:「最適化の理論と応用」研究部会
行列補完理論のポイント どれだけの要素が必要か⇒ Aggregate Sparsity Pattern どう行列を分解するか ⇒ Chordal Graph & 極大クリーク どう補完するか ⇒ 2012/12/15 SOTA:「最適化の理論と応用」研究部会

17 Aggregate Sparsity Pattern
双対側の変数行列の非ゼロパターン 1 5 3 2 4 6 7 グラフ表現 例 2012/12/15 SOTA:「最適化の理論と応用」研究部会

18 SOTA:「最適化の理論と応用」研究部会
Chordal Graph グラフが Chordal であるとは、4以上のサイクルに弦(chord)があること ここで、極大 Clique は (Clique は、すべての頂点間に枝がある部分のこと) Point :: 極大 Clique に沿って、行列が分解される 1 5 3 2 4 6 7 1 5 3 2 4 6 7 length4 length5 Chordal 化 2012/12/15 SOTA:「最適化の理論と応用」研究部会

19 行列の分解 1 5 3 2 4 6 7 2012/12/15 青はAggregate, 赤は Chordal で追加
Grone et al. 1984 2012/12/15 SOTA:「最適化の理論と応用」研究部会

20 Chordal Graph と Cholesky分解
1 5 3 2 4 6 7 1 5 3 2 4 6 7 2012/12/15 SOTA:「最適化の理論と応用」研究部会

21 行列の分解 1 5 3 2 4 6 7 2012/12/15 青はAggregate, 赤は Chordal で追加
Grone et al. 1984 黒の要素をどうする？ 2012/12/15 SOTA:「最適化の理論と応用」研究部会

22 SOTA:「最適化の理論と応用」研究部会
行列補完の簡単な例 1 2 3 正則転置 正定値 正則 2012/12/15 SOTA:「最適化の理論と応用」研究部会

23 SOTA:「最適化の理論と応用」研究部会
行列補完のステップ１ 1 5 3 2 4 6 7 2012/12/15 SOTA:「最適化の理論と応用」研究部会

24 SOTA:「最適化の理論と応用」研究部会
行列補完のステップ２ 下三角行列 これらを使うと、以下のように補完される 以下を満たすことを計算で確認できる 2012/12/15 SOTA:「最適化の理論と応用」研究部会

25 SOTA:「最適化の理論と応用」研究部会
行列補完の性質 補完行列は密だが逆行列が疎構造を持つ . 2012/12/15 SOTA:「最適化の理論と応用」研究部会

26 SOTA:「最適化の理論と応用」研究部会
Cholesky分解に似た分解 Point:: 2012/12/15 SOTA:「最適化の理論と応用」研究部会

27 SOTA:「最適化の理論と応用」研究部会
双対側は通常のCholesky分解 Sparse Cholesky 分解を適用 2012/12/15 SOTA:「最適化の理論と応用」研究部会

28 SOTA:「最適化の理論と応用」研究部会
行列補完と主双対内点法 主双対内点法のボトルネック Schur 補完行列 主問題側の探索方向 密行列を避けて疎行列で計算する 2012/12/15 SOTA:「最適化の理論と応用」研究部会

29 SOTA:「最適化の理論と応用」研究部会
Schur 補完行列の計算式の変更 2012/12/15 SOTA:「最適化の理論と応用」研究部会

30 SOTA:「最適化の理論と応用」研究部会
主問題側探索方向の計算 各列ごとに計算 2012/12/15 SOTA:「最適化の理論と応用」研究部会

31 SOTA:「最適化の理論と応用」研究部会
SDPA-C 行列補完理論を実装したSDPソルバー SDPA(SemiDefinie Programming Algorithm) – Completion version SDPA 従来の主双対内点法を実装 で公開 Debian/Ubuntu/Linux Mint パッケージ $ sudo apt-get install sdpa libsdpa-dev sdpam ソースコードは git でバージョン管理 ４年分の全部の記録を管理 （複数サーバでの数値実験に便利） SDPA と SDPA-Cで ソースを比較すると １万９千行のうち １万行ぐらい違う 2012/12/15 SOTA:「最適化の理論と応用」研究部会

32 SOTA:「最適化の理論と応用」研究部会
インデックス 半正定値計画問題の標準形と主双対内点法 行列補完理論と主双対内点法 実装の高速化 行列補完計算の改善 Schur complement マルチスレッド化 主探索方向マルチスレッド化 数値実験結果 2012/12/15 SOTA:「最適化の理論と応用」研究部会

33 SOTA:「最適化の理論と応用」研究部会
SDPA-CとSDPAの比較 SDPA-C 6.2.1 SDPA 7.3.8 SDPA-C 7.3.8 行列補完理論 ○ × Sparse Cholesky 分解 自作ルーチン MUMPS CHOLMOD &MUMPS Schur マルチスレッド dXマルチスレッド Callable Library Matlab 対応 SDPA-C 6.2.1 SDPA 7.3.8 SDPA-C 7.3.8 行列補完理論 ○ × Sparse Cholesky 分解 自作ルーチン MUMPS CHOLMOD &MUMPS Schur マルチスレッド dXマルチスレッド Callable Library Matlab 対応 2012/12/15 SOTA:「最適化の理論と応用」研究部会

34 SOTA:「最適化の理論と応用」研究部会
行列分解計算式の検討 . Point 2012/12/15 SOTA:「最適化の理論と応用」研究部会

35 SOTA:「最適化の理論と応用」研究部会
2012/12/15 SOTA:「最適化の理論と応用」研究部会

36 SOTA:「最適化の理論と応用」研究部会
2012/12/15 SOTA:「最適化の理論と応用」研究部会

37 これと同じことを極大クリークの情報で 帰納法を行えば一般的にできる クリークの情報にアクセスできる 数値計算ライブラリ
2012/12/15 SOTA:「最適化の理論と応用」研究部会

38 Sparse Cholesky 分解のライブラリ
CHOLMOD アルゴリズムは supernodal 分解 行列演算のマルチスレッドで性能が伸びにくい 行列補完理論に採用 MUMPS a MUltifrontal Massively Parallel sparse direct Solver アルゴリズムは multiple frontal method 行列演算のマルチスレッドで高速化される Schur complement 方程式の求解に採用 2012/12/15 SOTA:「最適化の理論と応用」研究部会

39 SOTA:「最適化の理論と応用」研究部会
Supernode 分解 . Sparse Cholesky 分解を効率的に行える supernode から クリークの情報が得られる 2012/12/15 SOTA:「最適化の理論と応用」研究部会

40 SOTA:「最適化の理論と応用」研究部会
内点法の前処理 2012/12/15 SOTA:「最適化の理論と応用」研究部会

41 SOTA:「最適化の理論と応用」研究部会
数式改善による高速化 格子上の最大クリーク問題 クリーク数 438, 平均29.89, 最大59, 合計 13090 SDPA-C6.2.0 SDPA-C7.3.8 Schur 要素計算 秒 秒 Schur Cholesky 185.36秒 161.03秒 dX の計算 316.00秒 242.51秒 その他 22.22秒 17.93秒 合計時間 秒 秒 1.88倍高速化 時間は短縮されたが、まだまだ遅い マルチスレッドによる並列計算で高速化 2012/12/15 SOTA:「最適化の理論と応用」研究部会

42 Schur 補完行列のマルチスレッド化 各列の計算は独立している スレッド割り当て⇒終わったスレッドが次の列へ 4 3 2 1 1 2 3
このままだとスレッドが衝突する 2012/12/15 SOTA:「最適化の理論と応用」研究部会

43 スレッドの衝突 行列積担当の BLAS もマルチスレッド そのままマルチスレッド化すると（４コアのとき） スレッド1 Schur 補完行列
スレッドが増えすぎて スレッドの切り替えに 時間を取られてしまう Schur 補完行列 スレッド2 スレッド3 スレッド4 列計算スレッド BLASスレッド 2012/12/15 SOTA:「最適化の理論と応用」研究部会

44 SOTA:「最適化の理論と応用」研究部会
スレッドの衝突回避 Schur 補完行列計算直前にBLASのスレッドを無効化 Schur 補完行列計算終了後にBLASのスレッドを戻す スレッド1 スレッド1 Schur 補完行列 スレッド2 スレッド2 スレッド3 スレッド3 スレッド4 スレッド4 列計算スレッド BLASスレッド 2012/12/15 SOTA:「最適化の理論と応用」研究部会

45 SOTA:「最適化の理論と応用」研究部会
dX計算のマルチスレッド化 各列ごとに計算 各列ごとにスレッドで並列計算 今回もBLASのスレッドを無効化 2012/12/15 SOTA:「最適化の理論と応用」研究部会

46 SOTA:「最適化の理論と応用」研究部会
マルチスレッドの効果 SDPA-C 6.2.0(1) SDPA-C 7.3.8(1) SDPA-C 7.3.8(2) SDPA-C 7.3.8(4) Schur 要素計算 秒 秒 秒 1.78倍 731.98秒 2.86倍 Schur Cholesky 185.36秒 161.03秒 85.24秒 1.88倍 47.77秒 3.37倍 dX の計算 316.00秒 242.51秒 131.03秒 1.85倍 83.54秒 2.90倍 その他 22.22秒 17.93秒 23.88秒 25.86秒 合計時間 秒 秒 秒 1.78倍 889.15秒 2.82倍 5.31倍高速化 （ ）の数字はスレッド数 最大クリーク問題 2012/12/15 SOTA:「最適化の理論と応用」研究部会

47 SOTA:「最適化の理論と応用」研究部会
インデックス 半正定値計画問題の標準形と主双対内点法 行列補完理論と主双対内点法 実装の高速化 行列補完計算の改善 Schur complement マルチスレッド化 主探索方向マルチスレッド化 数値実験結果 2012/12/15 SOTA:「最適化の理論と応用」研究部会

48 SOTA:「最適化の理論と応用」研究部会
数値実験環境とテスト問題 CPU Xeon X5365(3.0GHz), Memory 48GB, RedHat Linux テスト問題1 格子ネットワーク(p,q)で生成したSDP緩和 最大クリーク問題 最大カット問題 テスト問題2 量子化学におけるSpin Glassの基底状態エネルギー計算 DIMACS challenge の torusg-3-15のタイプ ↑ただし、前処理をする↑ 青が格子の情報で生成される, 赤が変数 2012/12/15 SOTA:「最適化の理論と応用」研究部会

49 テスト問題の Sparsity Pattern
最大カット問題 クリーク数 1607, 平均8.04, 最大26, 合計 12924 SpinGlass(3D) クリーク数 591, 平均29.97, 最大773, 合計 17714 2012/12/15 SOTA:「最適化の理論と応用」研究部会

50 最大クリーク問題(p=400,q=10) SDPA-C が非常に高速に最大クリーク問題を解く
クリーク数 438, 平均29.89, 最大59, 合計 13090 SDPA-C 6.2.0 SDPA 7.3.8(4) SDPA-C 7.3.8(4) SeDuMi 1.3 Schur 262.04 - Cholesky 403.28 248.55 dX 734.41 175.24 その他 31.37 47.65 合計 SeDuMi: SDPA-C が非常に高速に最大クリーク問題を解く 2012/12/15 SOTA:「最適化の理論と応用」研究部会

51 最大カット問題(p=500,q=10) 2012/12/15 SDPA-C 6.2.0 SDPA 7.3.8(4)
クリーク数 1607, 平均8.04, 最大26, 合計 12924 SDPA-C 6.2.0 SDPA 7.3.8(4) SDPA-C 7.3.8(4) SeDuMi 1.3 Schur 105.25 16.17 317.56 - Cholesky 502.43 214.43 226.06 dX 61.21 315.73 その他 14.38 17.46 合計 683.27 876.81 入力行列が単純なため、マルチスレッドのオーバーヘッドが大きい SDPA の Schur も非常に高速 2012/12/15 SOTA:「最適化の理論と応用」研究部会

52 SOTA:「最適化の理論と応用」研究部会
SpinGlass(p=15) クリーク数 591, 平均29.97, 最大773, 合計 17714 SDPA-C 6.2.0 SDPA 7.3.8(4) SDPA-C 7.3.8(4) SeDuMi 1.3 Schur 565.11 7.99 261.40 - Cholesky 39.47 6.64 11.19 dX 496.16 126.90 252.24 その他 209.26 194.87 35.570 合計 336.40 560.40 SDPA-C7.3.8は6.2.0よりも高速 ただしSDPAは、さらに高速 2012/12/15 SOTA:「最適化の理論と応用」研究部会

53 SpinGlass(p=25) このサイズではSDPAよりSDPA-Cのほうが高速 2012/12/15 SDPA-C 6.2.0
クリーク数 3173, 平均29.50, 最大1798, 合計 93609 SDPA-C 6.2.0 SDPA 7.3.8(4) SDPA-C 7.3.8(4) SeDuMi 1.3 Schur 220.19 - Cholesky 520.39 945.49 dX その他 516.45 合計 実験なし このサイズではSDPAよりSDPA-Cのほうが高速 2012/12/15 SOTA:「最適化の理論と応用」研究部会

54 SOTA:「最適化の理論と応用」研究部会
SpinGlass まとめ 計算時間の単位は秒 p n クリーク数 平均サイズ 最大サイズ SDPA7.3.8 SDPA-C7.3.8 10 1000 155 25.69 294 11.85 20.77 15 3375 191 29.97 773 336.40 560.40 18 5832 1118 28.13 913 20 8000 1737 25.75 1080 25 15625 3173 29.50 1798 SDPA-Cのほうが計算時間の上昇が緩やか 大きい問題ほどSDPA-Cが有利 ちなみに、CPU温度は Core i7 3770K で室温約15度の場合 通常 約30度, 最大 約80度, SDPA 約70度, SDPA-C 約60度 2012/12/15 SOTA:「最適化の理論と応用」研究部会

55 SOTA:「最適化の理論と応用」研究部会
まとめ 行列補完理論の数式を改善して高速化 大きい問題ほど効果が高くなる マルチスレッドでボトルネックを高速化 ただし、性能は疎構造に強く依存 2012/12/15 SOTA:「最適化の理論と応用」研究部会