5.3 Bečka の動的順序
5.3.1 アルゴリズム
Bečka の⼿法には 2 つの構成要素がある.1 つは𝐴 = 𝐴 𝐴 ⋯ 𝐴 に対して𝐶 = 𝐴 𝐴の⾮対⾓
ブロックの⼤きさ (フロベニウスノルム) を推定する⼿法,もう 1 つは,推定したノルムを⽤いて,貪 欲法によって,収束を早めてかつ並列に選択可能なペアを列挙する⼿法である.前者を求める理由は,
3 章の式 (3.36) の通り,ブロックヤコビ法の⼀次収束性が,消去した⾮対⾓ブロック分だけの⾮対⾓ブ ロックのフロベニウスノルムの⾃乗和が減少することによるためである.すなわち,最もフロベニウ スノルムの⼤きい⾮対⾓ブロックを選択すれば,局所的に最も収束を早めることになる.当然これは 最急降下法と同様に,あるステップにおいて局所的に最適な解を⽤いているだけであって,反復回数 を最⼩化するとは限らないが,計算可能かつ容易であり,また,実験的にもうまくいく⼿法である.
ここで問題となることが𝐶の計算であり,実際にすべての要素を計算するためには 1 並列ステップ の 1 ノード当たり約𝑛 /𝑝の演算が必要となるが,これは他の演算のオーダーである𝑂(𝑛 𝑘/𝑝 )より
も⼤きくなってしまう.そこで何らかの⽅法でこれを近似的に求める必要がある.
Bečka の⼿法では𝐶の⾮対⾓ブロックのフロベニウスノルム 𝐴 𝐴 を計算する代わりに,1 ベク トル𝟙を⽤いて𝑒 = 𝟙
‖𝟙‖ によって
𝑔, ≡ 𝐴 𝐴 𝑒 (5.6)
を計算する.この計算は,1 つの⾮対⾓ブロックに対して,𝑂(𝑛𝑙 )となる.⾮対⾓ブロックは ( ) 個あるが,𝑝並列で計算できるため,1 並列ステップの 1 ノード当たりの演算量は𝑂 となる.つ まり,⾮対⾓要素すべてを計算するよりも演算量を抑えることができる.また,次の定理により,𝑔 は 𝐴 𝐴 の良い近似になっていることが期待される.
定理 5.1(⼭本ら [120], Theorem 3, 4). 𝐶は𝑝 × 𝑞の実⾏列であり,𝑝 ≥ 𝑞かつ‖𝐶‖ = 𝑓を満たす.ま
た,ℳ , , は𝑝 × 𝑞のフロベニウスノルムが𝑓の実⾏列の集合であるとする.すなわち,正の対⾓⾏列
Σ = diag(𝜎 , 𝜎 , … , 𝜎 )と,𝑛次の直交⾏列の集合𝒪(𝑛)について,
ℳ , , = 𝑈 Σ
𝟘 , 𝑉 𝑈 ∈ 𝒪(𝑝), 𝑉 ∈ 𝒪(𝑞) (5.7)
に対して,
ℳ , , ≡ ∪ ( ) ℳ , , (5.8) と定義する.このとき,𝐶 ∈ ℝ × をℳ , , からランダムにサンプルされた⾏列だとすると,次を満 たす:
E ‖𝐶𝑒‖ 𝐶 ∈ ℳ , , = 𝑓 , (5.9)
V ‖𝐶𝑒‖ 𝐶 ∈ ℳ, , = 2
𝑞 + 2 𝜎 − 𝜎 . (5.10)
系 5.2(⼭本ら [120], Corollary 5). いま,特異値のサンプル平均 ̂𝑢とサンプル分散 ̂𝑠を次のように定義 する:
̂𝑢 = 1
𝑞 𝜎 , (5.11)
̂𝑠 = 1
𝑞 − 1 𝜎 − ̂𝑢 = 1
𝑞(𝑞 − 1) 𝜎 − 𝜎 . (5.12)
このとき,式 (5.10) より,
V ‖𝐶𝑒‖ 𝐶 ∈ ℳ , ,
E ‖𝐶𝑒‖ 𝐶 ∈ ℳ , , = 2 𝑞 + 2⋅
∑ 𝜎 − 𝜎
∑ 𝜎
= 2
𝑞 + 2⋅ 𝑞(𝑞 − 1) ̂𝑠 𝑞 ̂𝑢
= 2(𝑞 − 1) 𝑞(𝑞 + 2) ⋅ ̂𝑠
̂𝑢 ≈ 2 𝑞
̂𝑠
̂𝑢. (5.13)
以上の定理において,ランダムにサンプルされた⾏列𝐶は𝐴 𝐴の⾮対⾓ブロック𝐶, = 𝐴 𝐴 を想 定しており,定理における⾏列サイズ𝑝, 𝑞はそれぞれ𝐴 と𝐴 の列ブロック幅𝑛 と𝑛 である.そこ で,𝐶, がℳ , , からランダムにサンプルされたものであれば,𝑔, の期待値はほしかった値 𝐴 𝐴 と⼀致し,また,𝑔, の分散は,{𝜎 , 𝜎 , … , 𝜎 }の分散に⽐例してかつ√𝑞に反⽐例する.そのため,𝑛 が⼤きければ𝑔, はよい近似となっていると考えられる.ただし,𝐶, が本当にℳ , , からランダム にサンプルされたものと仮定できるのか,については議論の必要がある.実際には元の⾏列𝐶 = 𝐴 𝐴 がランダムにサンプルされたものであったとしても反復途中の⾏列の⾮対⾓ブロックがどのような構 造を持っているのか,解析は困難である.そのため⾮対⾓ブロックのフロベニウスノルムが過⼩評価 され,ヤコビ法の収束性が悪化したり,偽収束が発⽣したりする可能性がある.例えば,𝐶, の0に近 い固有値に対応する固有空間に𝟙ベクトルが含まれる場合,ノルムが過⼩評価される.しかしながら,
反復の途中でこのような過⼩評価される⾮対⾓ブロックがあったとして,𝐶, が消去されなかったと しても,他のブロックの更新によって𝐶, の値も書き換わり,固有ベクトルが変化する.そこで,過⼩
評価されたブロックが再び過⼩評価されるという状況が起きなければ,偽収束に⾄ることはない.ま た以下の実験においては偽収束は観測されなかった.
Bečka は他に数種類の重みの定義を作成している [121] が,収束性との関係が最もわかりやすいも のがこの定義である.他の定義においては,𝐴 と𝐴 それぞれの張る空間同⼠がどれだけ互いに傾斜 しあっているかの度合いを測るため,𝐴, 𝐴 の代わりに,それらの列を正規化した𝑈 と𝑈 を⽤い,
𝑈 𝑈 を近似的に計算する.この場合においても近似⼿法に上記と同じ⼿法を⽤いることができる.
このように重みを定義できたとき,並列に消去可能でかつ重みが⼤きいものを探索したい.いま重 み𝑔, をフロベニウスノルムの近似としたので,この⼆乗和が最⼤になるような並列ペアを探索でき れば,局所的に最も収束に近づくものを,近似的に求められると考えられる.これはグラフ理論の最
⼤重みマッチングを⽤いれば定式化できる.いまグラフ𝐺を𝑤個の頂点を持つ完全グラフと定義し,
頂点𝑖と𝑗の間の重み𝑔̄, = 𝑔, とする.このとき,𝐺上の最⼤重みマッチング𝑀 とは,𝐺上のす べてのマッチングの集合𝔐について
𝑀 = arg max
∈𝔐 ( , )∈
̄
𝑔, (5.14)
を求めることである.最⼤重みマッチングはマッチングの⼀種であり,マッチングは枝で結ばれる頂 点が 2 度以上現れないような,つまり頂点のペアを集めたグラフである.そこで頂点のペア(𝑖, 𝑗)を消 去する⾮対⾓ブロック(直交化する列ブロックのペア)だと考えると,最⼤重みマッチングによって 並列に消去可能でかつ重みの和が最⼤のペアを列挙できる.
最⼤重みマッチングは基本的なグラフ理論の問題の 1 つであり,Edmond による Blossom アルゴリ ズム [122] が有名である.Duan ら [123] に詳しいが,完全グラフ上の最⼤重みマッチングの計算量は
𝑂(𝑤 )となり,また,⾏列計算とは異なり,計算が複雑で並列化が難しく,また計算効率も低い.そ
こで,Bečka は最⼤重みマッチングの厳密解を求める代わりに,⾼速な近似解法によって近似解を求 める⽅法を選択している.これはそもそも重みが近似的に求められていること,最⼤重みマッチン グの厳密解であっても,ヤコビ法の反復計算全体の厳密解ではないことから,最⼤重みマッチング の厳密解を求める意味が薄いためである.最⼤重みマッチングの近似解法については Avis らによる Survey [124] や,Duan ら [123] に詳しい.Bečka は,貪欲法に基づく近似解法を⽤いている.これは 実装が単純でかつ,𝑂(𝑤 log 𝑤)の計算量であり,厳密解の重みの和に対して1/2以上の重みを持つ
1 subroutine Greedy-Matching([ ,, ,, ,, … , ,, … , , ]):
2 = []
3 = [ , , … , ]
4 = sort([ ,, ,, ,, … , ,, … , , ])
5 do
6 , ← the first element of .
7 remove , from .
8 if in and in :
9 remove and from .
10 append , to .
11 while is not empty
12 return
図5.17 貪欲法を⽤いた最⼤重みマッチングの近似解法の擬似プログラム
解を得られることが知られている [124, Section II.B].Duan らの⼿法のように,任意精度の解を得られ るものも存在するが,ここではアルゴリズムの単純さと実⾏速度を考えて,貪欲法を⽤いる.このア ルゴリズムを図 5.17に⽰す.