近接特異値を持つ行列に対応したI-SVD法の並列化とその評価

全文

(1)Vol.2009-MPS-74 No.10 2009/7/13. 情報処理学会研究報告 IPSJ SIG Technical Report. ters. It is shown that a more efficiency of the parallelization results for matrices having many clusters of medium or small size.. 近接特異値を持つ行列に対応した I-SVD 法の並列化とその評価豊川博己†1 高田雅美†2. 1. はじめに上二重対角行列の特異値分解を行うアルゴリズムの 1 つとして，I-SVD 法3),4),7),9) があ. 木村欣司†1 中村佳正†1. る．これは行列の特異値を mdLVs 法によって計算し，求められた特異値に対する特異ベクトルをツイスト分解など別の方法で計算するアルゴリズムである．このアルゴリズムは，特異値分解のアルゴリズムとして従来から使われてきた QR 法，分割統治法，二分法と逆反. 行列の特異値分解を行うアルゴリズムの 1 つとして I-SVD 法がある．これは，特異値計算に mdLVs 法，特異ベクトル計算に dLV 型ツイスト分解を用いたアルゴリズムであり，計算速度や相対精度などにおいて QR 法に基づく標準解法より優れている．最近の研究において，I-SVD 法は内部で用いられている mdLVs 法と dLV 型ツイスト分解の処理を適切に複数の計算コアに割り当てることによって，マルチコアプロセッサ上で効率的な並列化を行うことができた．しかし，近接特異値を持つ行列に対して dLV 型ツイスト分解を行うと，特異ベクトルの直交性が悪くなることが知られている．本研究では，あまり大きくない大きさのクラスターを複数持つ行列について特異ベクトルの直交性を改善することによって，並列化された I-SVD 法の改良を試みる．また数値実験を行い，改良された並列版 I-SVD 法を評価する．. 復法の組み合わせなどよりも，計算速度や計算結果の相対誤差などにおいて優れた点がある．しかし，mdLVs 法が逐次型アルゴリズムであるために，I-SVD 法の並列化は困難と考えられてきた．そこで，I-SVD 法を特異値計算と個々の特異ベクトル計算の処理の集合とみなし，それらの処理を複数の計算コアに割り当てることで並列化を行う方法が考案された8) ．この方法をマルチコアプロセッサ上で実行した場合，効率的な並列化を行うことができる．しかし近接特異値を持つ行列に対しては，ツイスト分解によって計算される特異ベクトルの直交性が悪化する場合があるという問題がある．本研究では，並列化された I-SVD 法の特異ベクトル計算方法を一部変更することで，近. On Parallelization of the I-SVD Algorithm and its Evaluation for Clustered Singular Values. 接特異値を持つ行列の特異ベクトルの直交性の改善を試みる．また，マルチコアプロセッサ上で数値実験を行い，その性能を評価する．. Hiroki Toyokawa,†1 Kinji Kimura,†1 Masami Takata†2 and Yoshimasa Nakamura†1. 2. 特異値分解と I-SVD 法 2.1 特異値分解とは正方行列 M ∈ Rm×m に対する特異値分解 (SVD) とは，M を次のように行列 U, Σ, V の. A parallelization of the I-SVD algorithm is performed. Here, the I-SVD algorithm is a new singular value decomposition (SVD) algorithm consisting of the mdLVs scheme and the dLV twisted factorization. For some random matrices, this algorithm is better than today’s standard algorithms for SVD, with some respects. However, if some singular values are not separated well and form clusters, the orthogonality of the singular vectors computed by the twisted factorization may be worse. In this paper, the parallelization of the I-SVD is modified to improve the orthogonality of such singular vectors in terms of a reorthogozalization process. In order to evaluate this modified algorithm, numerical experiments are conducted with a multi-core processor. As a result, the I-SVD algorithm is parallelized for matrices whose singular values form clus-. 積に分解することである．. M = U ΣV T ただし，U, V ∈ Rm×m は直交行列，Σ ∈ Rm×m は非負の要素を持つ対角行列である．Σ †1 京都大学大学院情報学研究科 Graduate School of Informatics, Kyoto University †2 奈良女子大学大学院人間文化研究科 Graduate School of Humanities and Sciences, Nara Women’s University. 1. c 2009 Information Processing Society of Japan ⃝.

(2) Vol.2009-MPS-74 No.10 2009/7/13. 情報処理学会研究報告 IPSJ SIG Technical Report. の k 番目の要素は，M の k 番目の特異値 σk であり，U, V の k 番目の列ベクトルは，それ. I-SVD. ぞれ k 番目の左特異ベクトル，右特異ベクトルである．V T は V の転置を表す．任意の正方行列に対して特異値分解が可能である1) ．また，任意の m × n 長方行列に対. mdLVs. しても特異値分解は可能であるが，本研究では正方行列についてのみ扱うことにする．任意. . .

(3) . の m × m 正方行列 M はハウスホルダー変換1) などによって上二重対角行列 B に変換可能であるから，本研究で扱うアルゴリズムは任意の行列に対して適用することができる．. dLV. . 2.2 mdLVs 法と dLV 型ツイスト分解 mdLVs 法3) は上二重対角行列 B = (bi,j ) ∈ Rm×m の特異値を効率よく計算するアルゴリズムの 1 つである．mdLVs 法の反復が進むと，B の各成分は次のように収束する．. ( lim n→∞. (n) bk,k. =. σk2. −. ∞ ∑. ) 12 θ. (l) 2. ,. lim n→∞. (n) bk,k+1. 図 1 I-SVD 法の処理の流れ Fig. 1 Process flow of the I-SVD. =0. (1) 3.1 行列の分割による mdLVs 法の並列化. l=1. (n). ここで，θ(n) はシフトと呼ばれる正の値である．数値安定性を保証するためには，シフ. mdLVs 法で反復計算を続けていくと，ある副対角成分 bm1 ,m1 が 0 に近いある正の数. トは B の最小特異値 σm よりも小さくなければならない．シフトによって mdLVs 法は高. εc ≈ 0 より小さくなることがある．このとき，その副対角成分は 0 に収束したと見なすこ. 速に収束するが，シフト計算の逐次性は並列化が困難な理由にもなっている．mdLVs 法は. とができ，そこを境に行列 B (n) ∈ Rm×m は 2 つの小さな上二重対角行列 B1. QR 法や二分法のようなアルゴリズムと比べて収束が速く，QR 法や dqds 法などに比べて. と B2. 得られた特異値の相対精度が高いという特徴がある．. 呼ぶ．. (n). (n). ∈ Rm2 ×m2 ,. . dLV 型ツイスト分解は，上二重対角行列の与えられた特異値に対する特異ベクトルを計算するためのアルゴリズムの 1 つである．逆反復法のような他のアルゴリズムに比べて，反復計算なしで高速に特異ベクトルを計算できる．. B. 2.3 I-SVD 法. (n).   = . ∗. ∗ ∗. m2 = m − m1 に分けることができる．これを行列 B (n) の分割と. . .     →  ∗ ∗ . ∗. ε. ∗. I-SVD 法は mdLVs 法と dLV 型ツイスト分解によって構成されている7),9) ．まず B の特異. ∈ Rm1 ×m1. ∗ ∗.    = ∗ ∗  0. (. ). (n). B1. (n). B2. (2). ∗. を dLV 型ツイスト分解で計算する（図 1）．dLV 型ツイスト分解は並列に実行することが. (n) (n) (n) B1 と B2 の特異値は mdLVs 法によって独立に計算可能である．そのため，B1 と (n) B2 の特異値計算を異なる計算コアに割り当てることで，mdLVs 法を並列化することが可. 可能なため，特異ベクトル計算の部分を並列化することは難しくない．しかし mdLVs 法は. 能である．. 値 {σk } を mdLVs 法で計算する．次にそれらの特異値に対応する特異ベクトル対 {(uk , vk )}. 逐次型アルゴリズムであるため並列化が困難である．そのため，特異値分解全体の並列化性. 3.2 行列の減次が起きた際の dLV 型ツイスト分解の並列実行行列の分割の m1 = m − 1,. 5). 能は，実質的に mdLVs 法の逐次性によって制限される．. m2 = 1 の場合を特に減次と呼ぶ．この場合は最下行と最 √ ∑n (n) 2 2 σm − l=1 θ(l) に収束する．よって，減次が. 右列が削られ，B の (m, m) 成分 bm,m が. 3. 並列化の方法. 起こると特異値 σm が得られる．. この節では，I-SVD 法を並列化する方法について説明する．. I-SVD 法では，特異値 σk が計算されると，それに対応する特異ベクトル (uk , vk ) の計算を他の特異値や特異ベクトルとは独立に開始することができる．そのため減次が起こった際. 2. c 2009 Information Processing Society of Japan ⃝.

(4) Vol.2009-MPS-74 No.10 2009/7/13. 情報処理学会研究報告 IPSJ SIG Technical Report. 3.4 近接特異値に対する特異ベクトル計算.

(5)

(6) . 複数の特異値が十分に分離しておらずクラスタを形成している場合は，それらに対する特異ベクトルを dLV 型ツイスト分解で計算すると，直交性が悪化する場合がある．特異ベク.

(7) . トルの直交性を改良するために，近接特異値に対する特異ベクトルは dLV 型ツイスト分解とは異なる方法で計算される必要がある6),11) ．まず，すべての特異値を mdLVs 法によって相対精度を保ったまま計算し，それらの中から近接特異値のクラスタを探す．ここで，クラ. . /'012 3.4. スタの最大サイズ（クラスタに含まれる特異値の個数の最大値）を mc とする．次に，それ. /012 3.4. らのクラスタに含まれる近接特異値に対する特異ベクトルをすべて逆反復法で計算し，修正グラムシュミットの直交化法によってそれらを直交化する．この処理によって，近接特異値. !" .

(8) !"$# %'& ()*+,.-. に対する特異ベクトルの直交性は改善される．. . I-SVD 法の計算量は O(m2 ) である．しかしこの直交化処理によって，近接特異値に対す. 図 2 各コアに対するジョブの割り当て Fig. 2 Assigning jobs to each core. る特異ベクトルの計算に必要な計算量は O(m2c m) に増大する．どれくらい近接している特異値を「近接特異値」と見なすかの判定基準によって，全体の計算時間は大きく変わる可能性がある．. に，空いている計算コアを用いて対応する特異ベクトルの計算を開始することができる．. 3.5 マルチコアプロセッサ. 3.3 計算コアへのジョブの割り当て. マルチコアプロセッサとは，1 つのパッケージの中に複数の計算コアを含むプロセッサの. mdLVs 法と dLV 型ツイスト分解の両方を並列化するために，特異値分解を以下のよう. ことである．従来のマルチプロセッサシステムと比べて，マルチコアプロセッサは共有キャッ. なジョブの集合と見なすことにする．. (A). 行列の分割が起こるまで mdLVs 法によって特異値計算を行う. シュメモリによって効率良くメモリにアクセスできる．さらに，消費電力が小さいという利. (B). 求められた特異値に対する特異ベクトルを dLV 型ツイスト分解によって計算する. 点がある．. 3.6 プログラムの実装. ジョブが生成されると，未処理ジョブのリストに即座に追加される．この時点でアイドル. テストプログラムは以下のように実行される（図 3）．. 状態になっている計算コアがあれば，その計算コアがリストからそのジョブを取り出す．未処理ジョブのリストは待ち行列に似た働きをするが，ジョブ (A) はジョブ (B) よりも優先. (1). メインプロセスがワーキングスレッドを生成し，それらを各計算コアに割り当てる．そして，未処理ジョブのリストなどを準備する（初期化）．. 的に取り出されるようになっている．計算コアがジョブ (A) を取り出した場合は，分割が起こるまで mdLVs 法が実行される．. (2). ワーキングスレッドがジョブを実行する（mdLVs 法による特異値の計算・dLV 型ツ. 行列の分割が起こると，分割によって 2 つに分かれた行列それぞれに対するジョブ (A) が. イスト分解による特異ベクトルの計算）．この間メインプロセスは何もせず，すべて. 生成され，未処理ジョブのリストに追加される．行列の減次が起こると，求められた特異値. のジョブが終了するのを待機する．. に対するジョブ (B) が生成され，未処理ジョブのリストに追加される．. (3). すべての特異値が計算されたら，ワーキングスレッドの 1 つが近接特異値のクラスタを探す．そしてそれらに対する特異ベクトルを，逆反復法と修正グラムシュミット. 計算コアがジョブ (B) を取り出した場合は，dLV 型ツイスト分解が実行され特異ベクト. の直交化法を用いて計算するジョブを生成する（3.4 節の処理）．. ルが計算される．これを繰り返すことで，各ジョブは複数の計算コアで並列に実行され，I-SVD 法は並列. (4). すべての特異ベクトルが計算されたら，ワーキングスレッドは自ら終了し，メインプロセスは特異値分解の計算結果を返す（終了処理）．. 化される（図 2）．. 3. c 2009 Information Processing Society of Japan ⃝.

(9) Vol.2009-MPS-74 No.10 2009/7/13. 情報処理学会研究報告 IPSJ SIG Technical Report 表 1 テストコンピュータの性能 Table 1 Specification of the test computer. efg hFi. Fujitsu HX600 × 1 ノード Quad Core AMD Opteron 2.3GHz× 4 (4 コア × 4 プロセッサ) DDR2-667 32GBytes RedHat Enterprise Linux AS V4 g++ 3.4.6 and gfortran 4.1.2 (-O3 オプション) GotoBLAS 1.26. コンピュータ. CPU メモリ OS コンパイラ. BLAS. B. G8HI.

(10) LM; !#NPO0=R4. メインプロセスは初期化と終了処理のみを行う．初期化では，プログラムが実行されるコ. 6879 %&:';8<- =8>. 8;. %&('?"@BADC. ,+.-0/k]l. %&('"EF. ンピュータが持つ計算コアと同数のワーキングスレッドを生成する．ワーキングスレッドは. ,+.-0/.12354. 未処理ジョブのリストからジョブを取り出し，自律的にジョブの内容を実行する．. 4. 数値実験.

(11) "!#$ %&('. 表 1 にあるコンピュータ上で数値実験を行った．この数値実験では，表 2 にある 3 種類. )* . の行列を用いた．なお，B3 は以下のように表される．. .    B3 =    . ˆ W. . δ. ˆ W. δ. ... .. ... .. ˆ W. δ.    ,   . LM;0? "!#NPO=Q4. δ = 10−4. (3) 6879 %&:';8<- =8> %&('?"@BADC. ˆ W ˆ ) = {9, 8, . . . , 2, 1, 2, . . . , 8, 9}, diag(W. ˆ ) = {1, 1, . . . , 1} subdiag(W. %&('"EF. (4). 2. mdLVs 法でのステップサイズは δ (n) ≡ 1 とし，シフト θ(n) の計算には一般化ニュートン j HI. シフト（p = 2，減算なし）10) を用いた．このテストプログラムは C++と FORTRAN で記述した．プログラムの並列化には Pthreads (POSIX standard for threads) ライブラリ B. を使用しており，MPI や OpenMP は用いていない．. ;KJ0" : SUT8VWX. 表 3，表 4，表 5 は各行列サイズ m における I-SVD 法の計算時間を表す．以下では，1. USYZ[+]\^_,à, ;bcdX. コアで逐次実行した場合と比較した速度向上率のことを「相対速度」と表すものとする．. 図 3 テストプログラムの計算フロー Fig. 3 Process flow of the test program. 4.1 B1 （特異値が十分に分離している行列）の実験結果表 3，図 4 は B1 を特異値分解した際の計算時間と相対速度を表す．計算コアの個数が増. 4. c 2009 Information Processing Society of Japan ⃝.

(12) Vol.2009-MPS-74 No.10 2009/7/13. 情報処理学会研究報告 IPSJ SIG Technical Report 表 6 近接特異値の個数 Table 6 The number of the clustered singular values. 表 2 実験で用いた上二重対角行列 Table 2 The upper bidiagonal matrices used in the numerical experiments 行列. 対角成分 bk,k. 副対角成分 bk,k+1. B1 B2 B3. 2.001 ランダム 1, 2, . . . , 9. 2 ランダム 1, 10−4. B1 B2 num max num max 8500 0 134 100 17000 0 294 222 25500 0 475 329 num : 近接特異値の個数の合計 max : 近接特異値のクラスタの最大サイズ. サイズ m. σk の分布分離やや密集密集. 表 3 B1 の計算時間 Table 3 Computation time for B1 サイズ m. 1. 2. 8500 （相対速度）. 69.1 (1.00) 266.9 (1.00) 632.8 (1.00). 33.3 (2.07) 138.4 (1.93) 347.4 (1.82). 17000 （相対速度） 25500 （相対速度）. 計算コア数 4. 17.3 (3.99) 76.8 (3.48) 257.2 (2.46). 8. 16. 11.1 (6.24) 46.2 (5.78) 147.1 (4.30). 9.7 (7.14) 38.6 (6.91) 94.1 (6.73) 単位:sec. 加するに伴って，B1 の計算時間は減少していることが分かる．図 7 はサイズ 8500 × 8500 の B1 を特異値分解している最中に，実際に稼動している計算コアの個数の推移を表す（初期化と終了処理の計算時間は除かれている）．B1 の計算中には行列の分割は起こらない．そのためジョブ (B) は頻繁には現れず，特に特異値分解の始めの段階において，いくつかの計算コアがアイドル状態になっていることが分かる．. 4.2 B2 （近接特異値のクラスタを複数個含む行列）の実験結果表 4，図 5 は B2 を特異値分解した際の計算時間と相対速度を表す．図 8 はサイズ. 表 4 B2 の計算時間 Table 4 Computation time for B2 サイズ m. 1. 8500 （相対速度） 17000 （相対速度） 25500 （相対速度）. 71.8 (1.00) 281.3 (1.00) 697.9 (1.00). 計算コア数 2 4. 35.5 (2.03) 149.8 (1.88) 395.2 (1.77). 19.0 (3.78) 88.5 (3.18) 316.6 (2.20). 8500 × 8500 の B2 を特異値分解している最中に，実際に稼動している計算コアの個数 8. 16. 10.9 (6.58) 51.5 (5.46) 196.6 (3.55). 7.9 (9.13) 40.0 (7.04) 125.4 (5.56) 単位:sec. の推移を表す（初期化と終了処理の計算時間は除く）．表 6 から，B2 には小∼中サイズのクラスタが複数個存在することが分かる．クラスタが存在すると，0 に近い副対角成分 bk,k+1 が存在すると考えられる．そのため，行列の分割や減次が B1 に比べて多く起こる．そのため，B1 に比べてジョブ (B) が特異値分解の始めの段階で多く生成され，B2 の計算は B1 に比べて早く完了する．. 4.3 B3 （すべての特異値がいずれかのクラスタに含まれる行列）の実験結果表 5，図 6 は B3 を特異値分解した際の計算時間と相対速度を表す．図 9 はサイズ. 表 5 B3 の計算時間 Table 5 Computation time for B3 サイズ m. 1. 2. 8500 （相対速度） 17000 （相対速度） 25500 （相対速度）. 444 (1.00) 3335 (1.00) 14834 (1.00). 243 (1.83) 1728 (1.93) 11128 (1.33). 計算コア数 4. 139 (3.20) 1041 (3.20) 9011 (1.65). B3 num max 8500 500 17000 1000 25500 1500. 8500 × 8500 の B3 を特異値分解している最中に，実際に稼動している計算コアの個数の推移を表す（初期化と終了処理の計算時間は除く）．この行列では，すべての特異値がい 8. 16. 93 (4.79) 681 (4.90) 5035 (2.95). 86 (5.14) 641 (5.20) 2445 (6.07) 単位:sec. ずれかのクラスタに含まれている（表 6）．大きなクラスタに含まれる特異値に対応する特異ベクトルを計算するには，直交化のプロセスが必要なため，全体では非常に長い計算時間を要する（図 9）．そのため，B3 の計算時間は B1 , B2 に比べてかなり長くなる．. 5. 結. 論. 本研究では，近接特異値のクラスタを持つ行列にも対応できるように，既に考案されてい. 5. c 2009 Information Processing Society of Japan ⃝.

(13) Vol.2009-MPS-74 No.10 2009/7/13. 情報処理学会研究報告 IPSJ SIG Technical Report. る I-SVD 法の並列化法8) を改良した．また，その方法を数値実験によって検証した．その結果，特異値の分布によって分割や減次の起こる頻度が影響を受け，その結果，計算コアの稼動数の変動の様子が異なることが示された．また，近接特異値の個数やクラスタの最大サイズによって，並列版 I-SVD 法の計算時間が大きく変わることも分かった．. 参. 図 4 B1 を特異値分解する際の相対速度（対数目盛） Fig. 4 Relative speed for B1 (log scale). 考. 文. 献. 1) G. Golub and C. Van Loan : Matrix Computation, 3rd edition, John Hopkins Univ. Press, Baltimore (1996). 2) B. Groβer and B. Lang : An O(n2 ) algorithm for the bidiagonal SVD, Lin. Alg. Appl., Vol. 358, pp. 45-70 (2003). 3) M. Iwasaki and Y. Nakamura : Accurate computation of singular values in terms of shifted integrable schemes, Japan J. Indust. Appl. Math., Vol. 23, pp. 239-259 (2006). 4) 岩崎雅史, 阪野真也, 中村佳正 : 実対称 3 重対角行列の高精度ツイスト分解とその特異値分解への応用, 日本応用数理学会論文誌, Vol. 15, pp. 461-481 (2005). 5) T. Konda, M. Takata, M. Iwasaki, S. Tsujimoto and Y. Nakamura : A parallelization of singular value computation algorithm by the Lotka-Volterra System (in Japanese), IPSJ SIG Technical Rep., 2004-HPC-100, pp. 13-18 (2004). 6) T. Konda, M. Takata, M. Iwasaki and Y. Nakamura : A new singular value decomposition algorithm suited to parallelization and preliminary results, Proceedings of IASTED International Conference on Advances in Computer Science and Technology (ACST2006), pp. 79-85 (2006). 7) 中村佳正 : 可積分系の機能数理, 共立出版 (2006). 8) H. Toyokawa, K. Kimura, M. Takata and Y. Nakamura : On parallelism of the I-SVD algorithm with a multi-core processor, JSIAM Letters, to appear (2009). 9) M. Takata, K. Kimura, M. Iwasaki and Y. Nakamura : Performance of a new scheme for bidiagonal singular value decomposition of large scale, Proceedings of IASTED International Conference on Parallel and Distributed Computing and Networks (PDCN2006), pp. 304-309 (2006). 10) T. Yamashita, K. Kimura and Y. Nakamura : On subtraction free formula for the diagonal elements of inverse powers of symmetric positive definite tridiagonal matrices, in preparation (2009). 11) 片山幹基, 木村欣司, 高田雅美, 坪井洋明, 岩崎雅史, 中村佳正 : 特異値分解法 I-SVD における左特異ベクトル計算部の改善, 日本応用数理学会論文誌, Vol. 18, pp. 389-407 (2008).. 図 5 B1 を特異値分解する際の相対速度（対数目盛） Fig. 5 Relative speed for B2 (log scale). 図 7 B1 を特異値分解する際に稼動する計算コア個数の推移図 6 B3 を特異値分解する際の相対速度（対数目盛） Fig. 7 Changes of the number of working cores for B1 Fig. 6 Relative speed for B3 (log scale). 図 8 B2 を特異値分解する際に稼動する計算コア個数の図 9 B3 を特異値分解する際に稼動する計算コア個数の推移推移 Fig. 9 Changes of the number of working cores Fig. 8 Changes of the number of working cores for B3 for B2. 6. c 2009 Information Processing Society of Japan ⃝.

(14)