2013 年 2 月 18 日 ＠九州大学
最適化ワークショップ：広がっていく最適化

行列の数値的同時ブロック対角化アルゴリズム

前原 貴憲1,2

1. 国立情報学研究所
2. JST, ERATO, 河原林巨大グラフプロジェクト自己紹介• 名前：前原 貴憲（まえはら たかのり）• 略歴：–2012：東京大学 情報理工学系研究科 数理情報学専攻（指導教員：室田一雄教授）

2012–現在： 国立情報学研究所 ERATO 特任研究員• 興味： 最適化 ・ 数値計算

2/ 36

行列の同時ブロック対角化問題
Given：n × n 実正方行列 A1 , . . . , AN
Find：
直交行列 P s.t.
P ⊤ A1 P, . . . , P ⊤ AN P ：同じ形のブロック対角行列※ できるだけ細かい分解 を求めたい
A1

A2

AN

⇓ simultaneouslyP ⊤ A1 P P ⊤ A2 P P ⊤ AN P

3/ 36

本日の目標・発表のながれ• 同時ブロック対角化問題• モチベーション（応用先）• 分解手法• 課題・最近の進展

4/ 36

モチベーション（応用先）

5/ 36

同時ブロック対角化はどこであらわれるか• 量子物理 [Wigner 1931–]..
.
• 半正定値計画 [de Klerk-Dobre 2011] etc• 完全正計画 [Arima-Kim-Kojima 2012] etc• 非可換多項式計画 [Burgdorf-Klep-Povh 2011]• 構造物の分岐経路追跡 [Aiura-Kakimura-Murota 2011]• 独立成分分析 [Gutch-Krumsiek-Theis 2011]• 動的ネットワークの安定性解析 [Irving-Sorrentino 2012]..
.
古典的なテーマ，大規模問題 の前処理として最近注目

6/ 36

量子物理（特に結晶物理）[Wigner 1931–]

一般に：物理量の計算＝固有値計算

Hx = ϵx• H = α1 H 1 + · · · + αN H N
物理パラメタの線型結合• H1 , . . . , HN 同時ブロック対角化⇒ パラメタに寄らない普遍的構造• 結晶：高い対称性⇐⇒ 細かく分解可能（群表現論）
http://www.nasa.gov/centers/ames/multimedia/images/2006/crystal.html

7/ 36

動的ネットワークの安定性解析[Irving-Sorrentino 2012]

x(t) = A1 x(t − 1) + · · · + AN x(t − N )• x(t)：時刻 t で病気に感染している度合い• Ak ：k 日前の感染者からの影響
例： A1 ：第一近傍
A2 ：第二近傍..
.
• 低次元なら豊富な安定性解析が使える⇒ 低次元問題に分割＝同時ブロック対角化
http://www.dreamstime.com/stock-image-network-infection-image11173181

8/ 36

半正定値計画問題[Gatermann-Parrilo 2004], [Murota-Kanno-Kojima-Kojima 2010], ...

minimizesubject to⟨C, X⟩⟨Ai , X⟩ = bi (i = 1, . . . , N )X⪰O• 係数行列の同時ブロック対角化⇒ 問題のサイズ・非線型性低減

7 days → 7 mins
交差数最小化 [de Klerk-Dobre-Pasechnik 2009]9/ 36

独立成分分析[Jutten-Herault 1985], [Cardoso-Soulomiac 1993], ...
Y1

Y2

Yℓ

Given:
Find:

n 次元観測信号 X
互いに独立な信号 Y1 , . . . , Yℓ
従属U?

W

∗

P C• P =-Ij独立)

X = WY

http://www.kindredcocktails.com/

10/ 36

再掲 同時ブロック対角化はどこであらわれるか• 量子物理 [Wigner 1931–]..
.
• 半正定値計画 [de Klerk-Dobre 2011] etc• 完全正計画 [Arima-Kim-Kojima 2012] etc• 非可換多項式計画 [Burgdorf-Klep-Povh 2011]• 構造物の分岐経路追跡 [Aiura-Kakimura-Murota 2011]• 独立成分分析 [Gutch-Krumsiek-Theis 2011]• 動的ネットワークの安定性解析 [Irving-Sorrentino 2012]..
.
古典的なテーマ，大規模問題の前処理として最近注目

11/ 36

分解手法

12/ 36

分解手法：２つの研究の流れ• 数値線型計算：「SBD = 固有値計算の一般化」– 数値的に安定な計算– 最細性の保証はないが，実用的に十分– 代数に詳しくない人でも理解しやすい• 抽象代数：「SBD = 群表現の分解」– 古典的な手法– 最細分解の保証– 近年数値計算の手法と合流

13/ 36

手法の一覧
numerical linear algebraold

abstract algebra
one-by-one[folklore -1820]

The recipe[Schur 1905]

Jacobi-like[Bunse-Gerstner, Byers, Mehrmann 1990]

JADE[Cardoso-Souloumiacc 1993][Theis 2007]

MKKKM[Murota-Kanno-Kojima-Kojima
2010]
MM[Maehara-Murota 2012]new ?
14/ 36

数値線型計算の系列
numerical linear algebraoldabstract algebra
one-by-one[folklore -1820]

The recipe[Schur 1905]

Jacobi-like[Bunse-Gerstner, Byers, Mehrmann 1990]

JADE[Cardoso-Souloumiacc 1993][Theis 2007]

MKKKM[Murota-Kanno-Kojima-Kojima
2010]
MM[Maehara-Murota 2012]new ?
15/ 36

one-by-one method [２入力・可換]定理. ２つの可換な行列 A, X は同時対角化可能
証明. wlog. X = diag(x1 , . . . , xn ),AX − XA = [(xj − xi )aij ] = O∗ ∗x ∗ ∗
xy∗
XA「X を対角化して残った自由度で A を対角化する」

16/ 36

Jacobi-like method [２入力・可換](Bunse-Gerstner, Byers, Mehrmann 1990)

one-by-one は誤差に弱い ⇒ 最適化 の考えで同時対角化[]⊤⊤
minimize oﬀ (P AP ) + oﬀ (P BP )• Givens 回転（2 × 2 回転）を繰り返して目的関数減らす（行列が１つなら固有値計算法の Jacobi 法と一致）• 局所的には収束・大域的収束は不明（元々の研究動機：歪ハミルトン構造をもつシステム解析）

17/ 36

JADE [多入力・可換](Cardoso, Souloumiac 1993)[]⊤⊤
minimize oﬀ (P A1 P ) + · · · + oﬀ (P AN P )• Jacobi-like method の複数行列への拡張• 信号処理（独立成分分析）で必要に
1 次元信号源 = 同時対角化• 実用的に，非常にうまく いっている• 信号処理分野で標準的な手法の１つ（ソフトウェア公開 by Cardoso）

18/ 36

JADE’ [多入力・非可換 (ブロック対角化)](Theis 2007)[]⊤⊤
minimize oﬀ (P A1 P ) + · · · + oﬀ (P AN P )• 信号処理（独立空間分析）で必要に多次元信号源 = 同時ブロック対角化• 目的関数の性質：「最細分解」は局所最適解の１つ⇒ 非対角でも構わず対角化用ルーチンを使ううまくいかない例を発見 [Maehara–Gutch 2010]• 高重複度問題• 実型でない問題· · · （代数的手法で「苦労した」例）
19/ 36

抽象代数の系列
numerical linear algebra
oldabstract algebra
one-by-one[folklore -1820]

Schur lemma[Schur 1905]

Jacobi-like[Bunse-Gerstner, Byers, Mehrmann 1990]

JADE[Cardoso-Souloumiacc 1993][Theis 2007]

MKKKM[Murota-Kanno-Kojima-Kojima
2010]
MM[Maehara-Murota 2012]new ?
20/ 36対称性と同時ブロック対角化
行列 A1 , . . . , AN が 群 G の対称性をもつ：

Q⊤ Aj Q = Aj(Q ∈ G)G：直交行列からなる群定理．（Schur の補題）
- G を分解すると対応して Aj が分解される- 逆に，Aj が分解できるときは非自明な G が存在

G の分解：表現の既約分解• Schur lemma: 群論の標準的手法• MKKKM, MM: one-by-one の拡張（理論保証は代数的）
21/ 36

Schur lemma(Schur 1905)

Q⊤ Aj Q = Aj

(Q ∈ G)

1. A1 , . . . , AN の対称性 G を発見
2. G の既約分解を計算
3. A1 , . . . , AN を分解
1930 年付近で物理分野で流行りだす（Wigner が導入）
（cf. Heisenberg の行列力学 1925）• 典型的な結晶の「既約分解」の表を作成• 結晶に応じて表を組み合わせて既約分解（職人芸）
22/ 36

「対称性を人間が発見する」手法の限界
Q. これらの “対称性” は同じ？
3V × 3V 行列；(3i + j, 3k + l) 成分＝k 節点の第 l 次元が変化したときの i 節点の第 j 次元の変化

A. 違う！
• 幾何学的・組合せ的対称性は同じ• 人間の手による「対称性発見」ではカバーしきれない⇒ 数値抽象代数的手法の導入へ[Murota-Kanno-Kojima-Kojima 2010]

23/ 36MKKKM[Murota-Kanno-Kojima-Kojima 2010, Maehara-Murota 2011]

1. 行列を１つランダムサンプリングして対角化（A1 , . . . , AN の和・積・スカラー倍・転置）
2. 対角化する行列を全体に適用
3. 残りの行列を残った自由度でブロック対角化
cf: one-by-one
1. １つの行列を対角化
2. 対角化する行列を全体に適用
3. 残りの行列を残った自由度で対角化
24/ 36MKKKM 抽象代数的理解[Murota-Kanno-Kojima-Kojima 2010, Maehara-Murota 2011]

1. 行列を１つランダムサンプリングして対角化＝行列 ∗ 代数の代表元を取得
2. 対角化する行列を全体に適用＝単純成分分解

3. 残りの行列を残った自由度でブロック対角化＝規約成分分解

cf: one-by-one
1. １つの行列を対角化
2. 対角化する行列を全体に適用
3. 残りの行列を残った自由度で対角化
25/ 36MKKKM 正当性議論の流れ• A1 , . . . , AN のブロック対角化の構造「Artin-Wedderburn 構造定理」で完全に記述
定理．T := ⟨A1 , . . . , AN ⟩ ≃(Mn1 ⊗ Iµ1 ) ⊕ · · · ⊕ (Mnℓ ⊗ Iµℓ )• T からランダムサンプル：
ランダム行列の固有値は重複しない「固有値分解の一意性」より残る自由度： Mn ⊗ Iµ
※ one-by-one 同様，誤差に弱い

26/ 36

現行最良：MM 法（アイデア）[Maehara-Murota 2012]

one-by-one の別方向への拡張；X = diag(x1 , . . . , xn ),AX − XA = [(xj − xi )aij ] = O∗ ∗x ∗ ∗
xy∗
XA• A1 , . . . , AN と可換な対称 X を「見つけられたら」その固有値構造に応じて A1 , . . . , AN はブロック対角に• 固有値がバラバラ→細かい分解：ランダムサンプリング
27/ 36

MM[Maehara-Murota 2012]

1. Ai X − XAi = O (i = 1, . . . , N ) をみたす
対称行列をランダムサンプリング
2. X を対角化する直交行列を出力• ランダムサンプル：固有値が最も分離する行列（厳密な証明：Artin-Wedderburn）• 代数的には：交換子代数上の計算
T ′ := {X | Ai X − Xi A = O (i = 1, . . . , N )}⇒ 「対称性の群」を数値的に計算
ランダムサンプル：最も一般的な対称性
固有値分解：群表現の既約分解

28/ 36

誤差制御版 MM[Maehara-Murota 2012]

1. ∥Ai X − XAi ∥ ≤ ϵ (i = 1, . . . , N ) をみたす
対称行列をランダムサンプリング
2. X を対角化する直交行列を出力• 許容誤差に応じた分解（非対角ブロック成分に誤差）• 非対角を丸めると ϵ 近似ブロック対角行列• 適切な ϵ：方程式系の最小特異値（正規誤差を仮定 ⇒ 最小特異値 ∝ 誤差）

29/ 36

http://www.misojiro.t.u-tokyo.ac.jp/˜maehara/commdec/30/ 36

ここまでのまとめ行列の同時ブロック対角化問題
A1A2AN⇓ simultaneously
P ⊤ A1 P P ⊤ A2 P P ⊤ AN P

• 古典的問題，問題の大規模化に伴い近年活発に応用• 手法：数値線型計算型，抽象代数型
state of the art: 交換方程式の求解 [MM 2012]

31/ 36

課題・最近の進展

32/ 36

課題・最近の進展• 効率的なアルゴリズム– 前処理（preconditioning）– 交互射影型アルゴリズム• ランダム行列理論による評価– 定性的評価

33/ 36

ランダム行列理論による評価
Ai X − XAi = O (i = 1, . . . , N ) をみたす
対称行列をランダムサンプリング

X の固有値分布，特に |λi − λj | の分布を評価したい⇒ ランダム行列の理論の結果が使えそう
定理 [Dyson index]．
n × n 各成分標準正規のランダム対称行列の固有値密度∝ exp(−n∑λ2i /4)∏|λi − λj |i<j（i.e., λi ≃ λj になる確率は特に低い；反発効果）（cf. ランダムエルミートにすると |λi − λj |2 ）34/ 36

ランダム行列理論による評価
T ：行列 ∗ 代数，単純分解 T = T1 ⊕ · · · ⊕ Tℓ
X ∈ T の固有値分布（with 適切なサンプリング）：• 各ブロックのノルム ∝ n2k• 各ブロック内の分布 ∝ Dyson index• 全体：これらの同時分布⇒• ブロック数多は難しい（反発効果の無い固有値対増）• ブロックサイズ比大は難しい（固有値がゼロ付近に密集）• 複素数上で分解したほうが安定度が高い（反発効果の指数）（ToDo: もう少し定量的な結果？）

35/ 36

まとめ
行列の同時ブロック対角化問題
A1A2AN⇓ simultaneouslyP ⊤ A1 P P ⊤ A2 P P ⊤ AN P• 古典的問題，問題の大規模化に伴い近年活発に応用• 手法：数値線型計算型，抽象代数型
state of the art: 交換方程式の特異値分解 [MM 2012]• 課題・最近の進展(以上)
36/ 36