九州大学学術情報リポジトリ

全文

(1)

九州大学学術情報リポジトリ

Kyushu University Institutional Repository

超平面制約付き非線形方程式に対するニュートン反 復を利用した特異値分解法

矢谷, 健一

京都大学大学院情報学研究科

近藤, 弘一

同志社大学理工学部

岩崎, 雅史

京都府立大学生命環境学部

https://doi.org/10.15017/18719

出版情報:応用力学研究所研究集会報告. 21ME-S7 (29), 2010-03. Research Institute for Applied Mechanics, Kyushu University

バージョン:

権利関係:

(2)

応用力学研究所研究集会報告No.21ME-S7

「非線形波動研究の現状と将来 次の10年への展望」(研究代表者 矢嶋 徹)

共催 九州大学グローバルCOEプログラム

「マス・フォア・インダストリ教育研究拠点」

Reports of RIAM Symposium No.21ME-S7

Current and Future Research on Nonlinear Waves — Perspectives for the Next Decade

Proceedings of a symposium held at Chikushi Campus, Kyushu Universiy, Kasuga, Fukuoka, Japan, November 19 - 21, 2009

Co-organized by

Kyushu University Global COE Program

Education and Research Hub for Mathematics - for - Industry

Research Institute for Applied Mechanics Kyushu University

March, 2010 Article No. 29 (pp. 191-196)

超平面制約付き非線形方程式に対する ニュートン反復を利用した特異値分

解法

矢谷 健一 (YADANI Kenichi), 近藤 弘一 (KONDO Koichi), 岩崎 雅史 (IWASAKI Masashi)

(Received January 25, 2010)

(3)

超平面制約付き非線形方程式に対するニュートン反復を利用した特異値分解法

京都大学大学院情報学研究科 矢谷健一 (YADANI Kenichi) 同志社大学理工学部 近藤弘一 (KONDO Koichi) 京都府立大学生命環境学部 岩崎雅史 (IWASAKI Masashi)

概 要 行列の特異値問題を非線形方程式の求解問題に置き換え,これをニュートン法に似た反復 計算で数値的に解くことによって特異値分解を得る方法を紹介する.さらに,この特異値分解法のい くつかの理論的性質を明らかにし,実用的なハイブリッド版を示す.

1 はじめに

長方行列の特異値分解を数値的に計算する方法として,超平面制約付き非線形方程式に対する 求解手法を利用する方法が提案されている[1].超平面制約法と名付けられたこの特異値分解法は,

特異値と特異ベクトルの組を計算するためにニュートン型反復を採用し,すべての特異値と特異 ベクトルの組を求めるための仕組みも備えている.2節では,ニュートン型反復とその理論的性 質,3節では,提案法である超平面制約法と,その収束定理について説明する.4節では超平面制 約法の応用として,既存の特異値分解法と組み合わせたハイブリッド版,および部分的に多倍長 計算を利用した方法を紹介し,5節で結論を述べる.

2 特異対を求めるためのニュートン型反復

長方行列A∈RRRm×n(m≤n)の特異値分解は

A=UΣVT, Σ:= (diag(σ1,σ2, . . . ,σm)On−m,m)∈RRRm×n

と定義される.ここで,非負なσkAの特異値と呼ばれる.また,U= (uuu1uuu2···uuum)∈RRRm×m,V = (vvv1vvv2···vvvm)∈RRRn×nは直交行列であり,列ベクトルuuuk,vvvkはそれぞれ,Aの左特異ベクトル,右特 異ベクトルと呼ばれる.本稿では,特異値σkと対応する左特異ベクトルuuuk,右特異ベクトルvvvkの 組(σk,uuuk,vvvk)を特異対と呼ぶことにする.

右特異ベクトルvvv1,vvv2, . . . ,vvvmが既知ならば,残りの右特異ベクトルvvvm+1,vvvm+2, . . . ,vvvnは,既知の右 特異ベクトルが張る線形空間span(vvv1,vvv2, . . . ,vvvm)の直交補空間の正規直交基底として,例えばグラム- シュミット法などによって容易に得られる.そのため,以降m組の特異対(σ1,uuu1,vvv1),(σ2,uuu2,vvv2), . . . , (σm,uuum,vvvm)の計算方法を中心に議論する.

まず,1組の特異対を求めるための反復計算法を説明する[1].特異対(σk,uuuk,vvvk)は,連立非線 形方程式

Avvvuuu, ATuuuvvv, ∥uuu∥2=1, ∥vvv∥2=1 (2.1) の解(σ,uuu,vvv)になることに注意する.(2.1)の第3, 4式は,ベクトルuuu,vvvの正規化条件であるため,

別の制約条件を課して得られたuuu,vvvを最後に正規化してもよい.ここで,正規化条件の代わりに,uuu にのみ超平面ΓU={uuu|(zzzU,uuu) =CU}上に存在するという制約を課す.ただし,zzzU∈RRRm,CU∈RRR\{0} は超平面を定める任意パラメータであり,zzzUは超平面の法線ベクトルとなる.すなわち,(2.1)の 代わりに(m+n+1)変数の連立非線形方程式

Avvvuuu, ATuuuvvv, (zzzU,uuu) =CU (2.2) 1

(4)

を考えることにする.このとき,(σ,uuu,vvv) = (σk,αkuuuk,αkvvvk),αk=CU/(zzz,uuuk) は(2.2)の解となる.

(2.2)からσ を消去すると,

HH

HU(xxx) =000, HHHU(xxx):=

(

Avvv−σU(vvv)uuu ATuuu−σU(vvv)vvv

)

, xxx:=

( uuu vvv

)

, σU(vvv):=(vvv,wwwU)

CU , wwwU:=ATzzzU (2.3) を得る.非線形方程式(2.3)を数値的に解くためのニュートン反復は,を反復回数とすると

xxx(ℓ+1)=ΦΦΦU(xxx(ℓ)), xxx(ℓ):=

( uuu(ℓ) vvv(ℓ)

)

, =0,1,2, . . . , (2.4)

ΦΦΦU(xxx):=xxx−(JU(xxx))1HHHU(xxx), JU(xxx):=∂HHHU

xxx =

(−σU(vvv)Im A−uuuwwwUTCU1 AT −σU(vvv)In−vvvwwwUTCU1

)

(2.5) となる.なお,JU(xxx)はHHHU(xxx)のヤコビ行列であり,JU(xxx)の正則性および初期ベクトルxxx(0)の与 え方については後述する.

続いて,uuuではなくvvvにのみ制約条件を課した場合に得られる非線形方程式と,その解を求める ためのニュートン反復を導出する.任意パラメータzzzV ∈RRRn,CV ∈RRR\{0}を用いて,zzzV を法線ベク トルとする超平面ΓV ={vvv|(zzzV,vvv) =CV}を定め,(2.1)の第3, 4式の代わりに,vvvが超平面ΓV上に 存在するという制約を課す.このとき,(2.1)とは異なる連立非線形方程式

Avvvuuu, ATuuuvvv, (zzzV,vvv) =CV (2.6)

が得られ,(σ,uuu,vvv) = (σk,βkuuuk,βkvvvk),βk=CV/(zzzV,vvvk)は(2.6)の解となる.(2.6)からσを消去すると HHHV(xxx) =000, HHHV(xxx):=

(

Avvv−σV(uuu)uuu ATuuu−σV(uuu)vvv

)

, xxx:=

( uuu vvv

)

, σV(uuu):=(uuu,wwwV)

CU , wwwV :=AzzzV (2.7)

であり,(2.7)を解くためのニュートン反復は,

xxx(ℓ+1)=ΦΦΦV(xxx(ℓ)), xxx(ℓ):=

( uu u(ℓ) vvv(ℓ)

)

, =0,1,2, . . . , (2.8)

Φ Φ

ΦU(xxx):=xxx−(JV(xxx))1HHHV(xxx), JV(xxx):=∂HHHV

xxx =

(−σV(uuu)Im−uuuwwwVTCV1 A AT−vvvwwwTVCV1 −σV(uuu)In

)

となる.(2.4)同様,JV(xxx)はHHHV(xxx)のヤコビ行列であり,JV(xxx)の正則性および初期ベクトルxxx(0) の与え方については後述する.

以上より,特異対を計算するための2種類のニュートン反復(2.4), (2.8)が得られた.どちらの ニュートン反復でも,ℓ→においてuuu(ℓ)は左特異ベクトル,vvv(ℓ)は右特異ベクトルに収束するこ とが期待される.

ニュートン反復によって生成されるベクトル列{uuu(ℓ)},{vvv(ℓ)}の存在範囲に関する定理を次に示す.

定理1 (i) (2.4)において,uuu(0)ΓUならばuuu(ℓ)ΓU, ℓ=1,2, . . .. (ii) (2.8)において,vvv(0)ΓV ならばvvv(ℓ)ΓV, ℓ=1,2, . . ..

計算機上では,数値誤差のため,定理1がほとんど成立しないことに注意したい.ニュートン反復 (2.4), (2.8)に代えて,計算機上でもuuu(ℓ)ΓU,vvv(ℓ)ΓV が常に成立する反復式を導入する.ここで,

ΨΨΨU(xxx):= CU

(zzzU,uuu)xxx, ΨΨΨV(xxx):= CV

(zzzV,vvv)xxx

(5)

と定義された超平面ΓU,ΓV への射影作用素ΨΨΨU,ΨΨΨV を(2.4), (2.8)に追加すると,それぞれ,U型 ニュートン反復,V型ニュートン反復と呼ばれる

xxx(ℓ+1)=ΨΨΨU(ΦΦΦU(xxx(ℓ))), =0,1,2, . . . , (2.9) xxx(ℓ+1)=ΨΨΨV(ΦΦΦV(xxx(ℓ))), =0,1,2, . . . (2.10) が得られる.(2.9), (2.10)で生成されるベクトル列{uuu(ℓ)},{vvv(ℓ)}については,明らかに次の定理が 成立する.

定理2 (i) (2.9)において,uuu(ℓ)ΓU, ℓ=1,2, . . ..

(ii) (2.10)において,vvv(ℓ)ΓV, ℓ=1,2, . . ..

ニュートン型反復(2.9), (2.10)を利用した,特異対を計算するためのアルゴリズムは,以下のと おりである[1].

アルゴリズム1 (U型ニュートン型反復法)

1. function(σ,uuu,vvv) =hppair u(A,zzzU,uuu(0),vvv(0)) 2. for=0,1,2, . . .till convergence

3. xxx(ℓ+1):=ΨΨΨU(ΦΦΦU(xxx(ℓ))) 4. end for

5. if(σU(vvv(ℓ+1))<0)thenvvv(ℓ+1):=−vvv(ℓ+1)

6. uuu:=uuu(ℓ+1)/∥uuu(ℓ+1)2,vvv=vvv(ℓ+1)/∥vvv(ℓ+1)2,σ:=uuuTAvvv アルゴリズム2 (V型ニュートン型反復法)

1. function(σ,uuu,vvv) =hppair v(A,zzzV,uuu(0),vvv(0)) 2. for=0,1,2, . . .till convergence

3. xxx(ℓ+1):=ΨΨΨV(ΦΦΦV(xxx(ℓ))) 4. end for

5. if(σV(uuu(ℓ+1))<0)thenuuu(ℓ+1):=−uuu(ℓ+1)

6. uuu:=uuu(ℓ+1)/∥uuu(ℓ+1)2,vvv=vvv(ℓ+1)/∥vvv(ℓ+1)2,σ:=uuuTAvvv

アルゴリズム1, 2に関する4つの定理を続いて述べるが,証明については[2, 4]を参照されたい.

ニュートン型反復の途中で,HHHU(xxx),HHHV(xxx)のヤコビ行列JU(xxx),JV(xxx)の逆行列が現れるため,

JU(xxx),JV(xxx)の正則性を示しておく必要がある.ヤコビ行列JU(xxx),JV(xxx)の正則性および逆行列の 有界性について,次の定理が成り立つ.

定理3 (ヤコビ行列の正則性および逆行列の有界性)

(i) JU(xxx)が正則であるための必要十分条件は,

σU(vvv)nm { m

i=1

i2σU(vvv)2) 1 CU

m i=1

[(zzzU,uuuii((uuui,uuu)σi+ (vvvi,vvv)σU(vvv))

j̸=i

2j σU(vvv)2)]}

̸

=0 である.また,このとき(JU(xxx))−1は有界.

3

(6)

(ii) JV(xxx)が正則であるための必要十分条件は,

σV(vvv)nm { m

i=1

i2σV(uuu)2) 1 CV

m i=1

[(zzzV,vvvii((vvvi,vvv)σi+ (uuui,uuu))σV(uuu)

=i

2j σV(vvv)2)]}

̸

=0 である.また,このとき(JV(xxx))1は有界.

定理3から直ちに,解の点xxxkxxxk,βkxxxkにおけるヤコビ行列JUkxxxk),JVkxxxk)の正則性およ び逆行列(JUkxxxk))1,(JVkxxxk))1の有界性に関する定理が得られる.

定理4 (解の点におけるヤコビ行列の正則性および逆行列の有界性)

(i) Aが零および重複特異値をもたず,zzzUが(zzzU,uuuk)̸=0をみたすとき,JUkxxxk)は正則であり,

(JUkxxxk))1は有界である.

(ii) Aが零および重複特異値をもたず,zzzV が(zzzV,vvvk)̸=0をみたすとき,JVkxxxk)は正則であり,

(JVkxxxk))1は有界である.

よく知られているように,一般的なニュートン法は2次収束性をもつ.ニュートン型反復(2.9),

(2.10)に関しても,2次収束性が証明される.

定理5 (i) U型ニュートン型反復(2.9)において,(σU(vvv(ℓ)),uuu(ℓ),vvv(ℓ))は,いずれかの特異対に局 所2次収束する.

(ii) V型ニュートン型反復(2.10)において,(σV(uuu(ℓ)),uuu(ℓ),vvv(ℓ))は,いずれかの特異対に局所2次 収束する.

また,ニュートン型反復(2.9), (2.10)の任意パラメータzzzU,zzzV をそれぞれ適切に選ぶと,ある特 定の特異対を求めることができる.

定理6 (i) 超平面ΓUの法線ベクトルzzzUUk=span(uuu1,uuu2, . . . ,uuuk)の直交補空間Ukから選ぶ と,U型ニュートン型反復(2.9)において,xxx(ℓ)の極限xxx= ((uuu)T(vvv)T)Tuuuは,uuu∈/Uk

をみたす.

(ii) 超平面ΓV の法線ベクトルzzzVVk=span(vvv1,vvv2, . . . ,vvvk)の直交補空間Vkから選ぶと,V型 ニュートン型反復(2.10)において,xxx(ℓ)の極限xxx= ((uuu)T(vvv)T)Tvvvは,vvv∈/Vkをみたす.

3 特異値分解を求めるための超平面制約法

非線形方程式のすべての解をニュートン法によって求めるには,解の数だけ適切な初期値を用 意しなければならないが,一般的には極めて困難である.ところが,定理6に示したように,任 意パラメータzzzU,zzzV を適切に選べば,ニュートン型反復の収束先を制限できることに注意したい.

すなわち,zzzU,zzzV を,それぞれ既得の左特異ベクトル,右特異ベクトルと直交するように選びなが らニュートン型反復を繰り返すと,すべての特異対を計算できる.よって,アルゴリズム1, 2を 利用した特異値分解アルゴリズムは以下のようになる.

アルゴリズム3 (U型超平面制約法) 1. function(Σ,U,V)=hpsvd u(A)

(7)

2. fork=1,2, . . . ,m

3. Randomly selectzzzU from the orthogonal complement ofUk1=span(uuu1,uuu2, . . . ,uuuk1) 4. Randomly selectuuu(0)∈RRRm,vvv(0)∈RRRn

5. (σk,uuuk,vvvk):=hppair u(A,zzzU,uuu(0),vvv(0)) 6. if(hppair udid not converge)then goto3 7. end for

8. Compute the orthonormal basis of{vvv1,vvv2, . . . ,vvvm}and set them as{vvvm+1,vvvm+2, . . . ,vvvn} 9. Σ:= (diag(σ1,σ2, . . . ,σm)Om,nm),U:= (uuu1uuu2 ··· uuum),V := (vvv1vvv2 ··· vvvn)

アルゴリズム4 (V型超平面制約法) 1. function(Σ,U,V)=hpsvd v(A) 2. fork=1,2, . . . ,m

3. Randomly selectzzzV from the orthogonal complement ofVk1=span(vvv1,vvv2, . . . ,vvvk1) 4. Randomly selectuuu(0)∈RRRm,vvv(0)∈RRRn

5. (σk,uuuk,vvvk):=hppair v(A,zzzV,uuu(0),vvv(0)) 6. if(hppair vdid not converge)then goto3 7. end for

8. Compute the orthonormal basis of{vvv1,vvv2, . . . ,vvvm}and set them as{vvvm+1,vvvm+2, . . . ,vvvn} 9. Σ:= (diag(σ1,σ2, . . . ,σm)Om,nm),U:= (uuu1uuu2 ··· uuum),V := (vvv1vvv2 ··· vvvn)

定理6から,アルゴリズム3, 4について,以下の収束定理が成り立つ.

定理7 (超平面制約法の収束定理) Aは非零かつ重複でない特異値をもつとする.

(i) アルゴリズム3において,hppair uが常に収束するならば,Aの特異値分解が得られる.

(ii) アルゴリズム4において,hppair vが常に収束するならば,Aの特異値分解が得られる.

超平面制約法の数値特性については,[3]を参照されたい.

4 超平面制約法の高速・高精度化

ニュートン型反復は2次収束性をもつ(定理5)が,反復終盤の2次収束の段階に達するまでに 多くの反復が必要となることが数値的に確認されている.反復序盤の無駄を省くため,他の高速 特異値分解法とのハイブリッド法が提案されている[3].すなわち,既存の特異値分解アルゴリズ ムfast svdで得られた特異対によってxxx(0)を定め,ニュートン型反復を実行する.U型超平面 制約法とのハイブリッド法は以下のようになる.

アルゴリズム5 (ハイブリッド法)

1. function(Σ,U,V)=hybridsvd u(A)

2. (σ1,σ2, . . . ,σm,uuu1,uuu2, . . . ,uuum,vvv1,vvv2, . . . ,vvvm):=<fast svd>(A) 3. fork=1,2, . . . ,m

4. if((σk,uuuk,vvvk)is not accurate)then

5. Randomly selectzzzUfrom the orthogonal complement ofUk1=span(uuu1,uuu2, . . . ,uuuk1) 6. (σk,uuuk,vvvk):=hppair u(A,zzzU,uuuk,vvvk)

7. if(hppair udid not converge)then goto5 8. end if

9. end for

5

(8)

10. Compute the orthonormal basis of{vvv1,vvv2, . . . ,vvvm}and set them as{vvvm+1,vvvm+2, . . . ,vvvn} 11. Σ:= (diag(σ1,σ2, . . . ,σm)Om,nm),U:= (uuu1uuu2 ··· uuum),V := (vvv1vvv2 ··· vvvn)

V型超平面制約法とのハイブリッド法もアルゴリズム5とほとんど同じなので,ここでは省略す る.数値実験の結果から,超平面制約法と比べてハイブリッド法で得られる特異値分解の精度に 大きな変化は見られないが,ニュートン型反復の反復回数削減に成功していることが分かる[3]. 同時に,ハイブリッド法は既存法に比べて高精度な特異値分解法であることが確認されている.

この他,[5]では,ニュートン型反復の一部分にのみ多倍長計算を利用する超平面制約法および ハイブリッド法が提案されている.ニュートン型反復(2.9)において,uuu,vvvを多倍長数として保持 し,最も計算量の多い(2.5)中の(JU(xxx))1HHHUの計算は倍精度計算で,それ以外は多倍長で計算す る.ニュートン型反復以外の部分は倍精度演算のみを用いる.仮数部96ビットの多倍長数を部分 的に利用した場合,求められた特異値が真の特異値と,53ビットすべてが一致する,つまり,倍 精度の範囲で正確な特異値が得られることが数値的に示されている[5].

5 結論

本稿では,[1]で提案された,非線形方程式の数値解法に基づく行列の特異値分解法について,

その収束定理と高速・高精度化する方法を示した.2節では,特異対を計算するためのニュートン 型反復の理論的性質を述べ,特に,この反復が2次収束することを示した.3節では,ニュートン 型反復の任意パラメータを適切に選べば収束先が制限されることに着目し,超平面制約法の収束 定理を示した.4節では,ニュートン型反復の反復回数を削減するためのハイブリッド法,および,

部分的に多倍長精度計算を利用する超平面制約法およびハイブリッド法を紹介した.

参考文献

[1] 近藤弘一,杉本昌平,岩崎雅史: “非線形方程式の解法による行列の特異値分解アルゴリズム”, 日本応用数理学会論文誌19(2009), 81–103.

[2] K. Yadani, K. Kondo, M. Iwasaki: “A singular value decomposition algorithm based on solving hyperplane constrained nonlinear systems”, to appear in Appl. Math. Comp.

[3] K. Yadani, K. Kondo, M. Iwasaki: “Numerical performance of hyperplane constrained method and its hybrid method for singular value decomposition”, submitted.

[4] K. Yadani, K. Kondo, M. Iwasaki: “On the convergence of the V-type hyperplane constrained method for singular value decomposition”, to appear in JSIAM Letters.

[5] K. Yadani, K. Kondo, M. Iwasaki: “Mixed double-multiple precision version of hyperplane con- strained method for singular value decomposition”, to appear in JSIAM Letters.

Updating...

参照

Updating...

関連した話題 :