九州大学学術情報リポジトリ
Kyushu University Institutional Repository
固有値分解のための超平面制約法に対する理論解析
吉武, 奈緒美
京都府立大学生命環境学部
岩崎, 雅史
京都府立大学生命環境学部
近藤, 弘一
同志社大学理工学部
https://doi.org/10.15017/23417
出版情報：応用力学研究所研究集会報告. 22AO-S8 (33), pp.214-219, 2011-03. Research Institute for Applied Mechanics, Kyushu University
バージョン：
権利関係：
応用力学研究所研究集会報告No.22AO-S8
「非線形波動研究の新たな展開 — 現象とモデル化 —」（研究代表者 筧 三郎）
共催 九州大学グローバルCOEプログラム
「マス・フォア・インダストリ教育研究拠点」
Reports of RIAM Symposium No.22AO-S8
Development in Nonlinear Wave: Phenomena and Modeling Proceedings of a symposium held at Chikushi Campus, Kyushu Universiy,
Kasuga, Fukuoka, Japan, October 28 - 30, 2010 Co-organized by
Kyushu University Global COE Program
Education and Research Hub for Mathematics - for - Industry
Research Institute for Applied Mechanics Kyushu University
March, 2011 Article No. 33 (pp. 214 - 219)
固有値分解のための超平面制約法に対 する理論解析
吉武 奈緒美（ YOSHITAKE Naomi ），岩崎 雅史
（ IWASAKI Masashi ），近藤 弘一（ KONDO Koichi ）
（Received 16 January 2011）
固有値分解のための超平面制約法に対する理論解析
京都府立大学生命環境学部 吉武 奈緒美 (YOSHITAKE Naomi) 京都府立大学生命環境学部 岩崎 雅史 (IWASAKI Masashi) 同志社大学理工学部 近藤 弘一 (KONDO Koichi) 概 要 行列の固有値分解を求めるためのアルゴリズムとして超平面制約法が提案されている[2]． 特異値分解用の超平面制約法については数値安定性，収束性，計算精度などの解析が[4, 5, 6]におい て進められているが，固有値分解用に関してはアルゴリズムが定式化された段階で未解明な部分が少 なくない．本報告では，固有値分解用の超平面制約法に関するいくつかの性質を明らかにする．
1 はじめに
行列A∈Cn×nおよびスカラーλ∈Cに対して，
Ax=λx (1.1)
を満たす0でないベクトルx∈Cnが存在するならば，λがAの固有値，xがAの固有ベクトルとな る．固有値λと固有ベクトルxの組(λ,x)はAの固有対と呼ばれる．本報告では，n個の固有値を λ1,λ2, . . . ,λn，対応する固有ベクトルをx1,x2, . . . ,xnとして，固有対(λ1,x1),(λ2,x2), . . .(λn,xn) を求める固有値分解について考える．
固有値分解の応用分野は多岐にわたるため，高速かつ高精度に固有値分解を求めるためのアル ゴリズム開発は重要課題といえる．固有値分解のためのアルゴリズムとして，ハウスホルダーQR 法，分割統治法，ヤコビ法などがあるが，応用分野の要求を完全には満たせていないのが現状で
ある[1, 3]．著者らは，[2]において高精度に固有値分解を求めるためのアルゴリズムとして超平
面制約法を提案した．超平面制約法では，固有ベクトルxの存在範囲を超平面に制約してxに関 する非線形2次方程式を導き，それをニュートン型反復によって解く．[2]では超平面制約法が高 精度に固有値分解できることを数値的に示しているが，超平面制約法の性質に関する理論的な考 察はなされていない．本報告では，超平面制約法におけるニュートン型反復について再検討し，そ の収束性について明らかにする．2節では，まず固有値分解用の超平面制約法について概説する．
3節では，考え得るニュートン型反復をすべて挙げ，そのうち有効なニュートン型反復がどれにな るかを調べる．4節で有効なニュートン型反復の収束性について証明する．まとめと今後の課題は 5節で述べる．
2 固有値分解のための超平面制約法
例えば，(1.1)を満たすAの固有対(λ,x)を求める際に，制約条件kxk=1を付加した Ax=λx, kxk2=1
を解いてもよい．また，任意のα∈Rに対してαxもAの固有ベクトルとなるので，kxk2=1の 代わりに(z,x) =C(6=0)を課してもよい．つまり，
Ax=λx, (z,x) =C (2.1)
からAの固有対(λ,x)を求めてもよい．ただし，zは超平面P={x|(z,x) =Cの法線ベクトルで あり，設定方法については[2]を参照されたい．
以下に，(2.1)の解(λ,x)に関する定理を示す. 1
定理2.1 zが(z,xk)6=0, k=1,2, . . . ,nを満たすならば，(2.1)の解は (λ,x) = (λk,αkxk), αk= C
(z,xk), k=1,2, . . . ,n となる．
証明 (2.1)の第1式にλ =λk,x=αkxkを代入すると，Axk=λkxkなので(λk,αkxk)はAの固有 対である．(2.1)の第2式にx=αkxkを代入すると(z,αkxk) =Cであり，αk はスカラーなので αk(z,xk) =Cとなる．よって，αkについて整理するとαk=C/(z,xk)を得られる．
ここで，(2.1)の第1式の両辺に左側からzHをかける．ただし，上付き添え字のHは共役転置
を表す．λ はスカラーなので(AHz)Hx=λ(z,x)となる．(z,x) =Cを踏まえてλについて整理す ると，λ= (AHz,x)/Cが得られる．λ はxの関数とみなせるので，改めて
λ(x) =(w,x)
C , x=AHz
と書くことにする．このとき，Aの固有ベクトルxが解となる非線形方程式
F(x) =Ax−λ(x)x=0 (2.2)
が導かれる．超平面制約法では，(2.2)を解くためにニュートン型反復を使用する．具体的な漸化 式の1つは，
x(`+1)= xˆ(`+1)
(z,xˆ(`+1)), `=0,1, . . . , `max, ˆ
x(`+1)=x(`)−J(x(`))−1F(x(`)), J(x(`)) =A−λ(x(`))I−x(`)wH
C , λ(x(`)) =(w,x(`))
C .
(2.3)
である．ここで，J(x(`))はF(x(`))のヤコビ行列，Iはn×nの単位行列である．適当なx(0)に対 して，数列{x(`)}`=0,1,...が極限`→∞で収束するならば，lim`→∞x(`)は(2.2)の解となる．
3 最適なニュートン型反復の選択
(2.2)を解くために(2.3)ではF(x(`)),J(x(`)),λ(x(`))を
F(x(`)) =Ax(`)−λ(x(`))x(`), J(x(`)) =A−λ(x(`))I−x(`)wH
C , λ(x(`)) =(w,x(`))
C
(3.1)
としている．ここで，F(x(`)),J(x(`)),λ(x(`))の定め方は(3.1)以外に考え得ることに注意したい．
(3.1)と異なる4つの候補を以下に示す．
F(x(`)) =A Cx(`)
(z,x(`))−λ(x(`)) Cx(`) (z,x(`)), J(x(`)) = C
(z,x(`)) (
A−λ(x(`))I−(w,x(`)) (z,x(`))
) , λ(x(`)) =(w,x(`))
C .
(3.2)
F(x(`)) =Ax(`)−λ(x(`))x(`), J(x(`)) =A−λ(x(`))I+
(λ(x(`))z C −w
C )H
x(`), λ(x(`)) =(w,x(`))
(z,x(`)).
(3.3)
F(x(`)) =A x(`)
kx(`)k−λ(x(`)) x(`) kx(`)k, J(x(`)) = 1
kx(`)k(A−λ(x(`))I), λ(x(`)) =(w,x(`))
(z,x(`)).
(3.4)
F(x(`)) =A x(`)
kx(`)k−λ(x(`)) x(`) kx(`)k, J(x(`)) = A
kx(`)k− Ax(`)
kx(l)k3xH+λ(x(`))x(`)
kx(`)k3 xH− (Ax(`))H (x(`),x(`))
x(`) kx(`)k, λ(x(`)) = (x,Ax(`))
(x(`),x(`)).
(3.5)
以降，(3.1), (3.2), (3.3), (3.4), (3.5)のようにF(x(`)),J(x(`)),λ(x(`))を定めたニュートン型反復を それぞれタイプ1，タイプ2，タイプ3，タイプ4，タイプ5と呼ぶことにする．タイプ1では固 定された超平面上に，タイプ2では固定された任意の超平面上に固有ベクトルの存在範囲が制約 される．タイプ3では固有ベクトルの存在範囲が超平面上に制約されるが，超平面は反復ごとに 変化する．タイプ4では固有べクトルが球面に，タイプ5では固有対が球面に制約される．
各タイプにおいてヤコビ行列J(x(`))が正則と仮定すると，タイプ3，タイプ4，タイプ5は2回 目の反復でx(2)=0となる．よって，少なくともタイプ3，タイプ4，タイプ5ではすべての固有 対を求めることはできない．ゆえに，タイプ1，タイプ2が有効なニュートン型反復と結論づけら れる．タイプ1，タイプ2における{x(`)}`=0,1,...の存在範囲に関する定理は以下のとおりである．
定理3.1 x(`)が(z,x(`)) =Cを満たすとする．ヤコビ行列J(x(`)) が正則ならば，x(`+1)はまた
(z,x(`+1)) =Cを満たす．つまり，J(x(`))が超平面上に存在するならば，x(`+1)も超平面上に存在
する．
4 ニュートン型反復の収束性
ニュートン型反復に含まれるヤコビ行列J(x(`))を正則と仮定すると，以下の定理が得られる．
定理4.1 ヤコビ行列J(x(`))が正則ならば，kJ(x(`))−1k ≤M1となる正の定数M1が存在する．
証明 J(x(`))が正則なので
kJ(x(`))−1k ≤ 1
kdet(J(x))kkadj(J(x))k
であり，kdet(J(x))kとkadj(J(x))kは有界となる．よって，kJ(x(`))−1k ≤M1である．
一般的に，ニュートン法に含まれるヤコビ行列が正則ならば，ニュートン法は2次収束するこ とが知られている．ニュートン型反復が核となる超平面制約法は，ニュートン法と同様の収束性
3
をもつことが予想される．本節では，x(`)からx(`+1)へ変換を2つの写像Φ,Ψにわけて，タイプ 1のニュートン型反復の収束性を明らかにする．写像Φ,Ψは以下のように定義する．
Φ(x) =x−J(x)−1F(x), Ψ(x) = C
(z,x)x.
まず，F(x)とJ(x)に関する2つの補助定理を与える．
補助定理4.2 F(x), J(x)に対して，
F(x+∆x) =F(x) +J(x)∆x−λ(∆x)∆x (4.1) が成り立つ．
証明 F(x) =Ax−λ(x)xより，
F(x+∆x) = (Ax−λ(x)x) + (
A−λ(x)−x wH C
)
∆x−λ(∆x)∆x
なので，(4.1)が得られる．
補助定理4.3 k∆xk →0のとき，
Φ(x+∆x) =Φ(x)−J(x)−1J0(∆x)J(x)−1F(x)−J(x)−1λ(∆x)∆x
−(J(x)−1J0(∆x))2J(x)−1F(x) +η(kxk3) (4.2) が成り立つ．ただし，
J0(∆x) =λ(∆x) +∆xwH
C , (4.3)
η(kxk3) =
O(kxk3) O(kxk3)
... O(kxk3)
である．
証明 J(x) =A−λ(x)−1F(x)と(4.3)より，
J(x+∆x) =J(x) +J0(∆x) (4.4) となる．Φ(x) =x−J(x)−1F(x)なので，補助定理4.2と(4.4)より，
Φ(x+∆x) =x+∆x−(
J(x) +J0(∆x))−1
(F(x) +J(x)∆x−λ(∆x)∆x) (4.5) となる．(4.5)において(J(x) +J0(∆x))−1= (In−J(x)−1J0(∆x))−1J(x)−1であり，さらに，
(In−J(x)−1J0(∆x))−1J(x)−1=
∑
∞ i=0(J(x)−1J0(∆x))iJ(x)−1
なので，(4.2)が得られる．
次に，x(`+1)=Φ(x(`))によって生成される{x(`)}`=0,1,...の収束性に関する定理を導く．
定理4.4 (z,x)6=0とする．x(0)がαkxkに十分近いならば，x(`+1)=Φ(x(`))によって生成される
{x(`)}`=0,1,...は`→∞においてαkxkに2次収束する．
証明 x∗=αkxkとする．補助定理4.3においてx=x∗,∆x=x(`)−x∗とすると，F(x∗) =0,Φ(x∗) = x∗,x(`+1)=Φ(x(`+1))なので，
x(`+1)−x∗=−J(x)−1λ(x(`)−x∗)(x(`)−x∗) +η(kxk3) となる．ここで，定理4.1よりkJ(x)−1k ≤M1であり，
|λ(x(`)−x∗)|= wH
C (x(`)−x∗)(x(`)−x∗)
≤M2kx(`)−x∗k となる正の定数M2が存在するので，
kx(`+1)−x∗k ≤M1M2kx(`)−x∗k2 (4.6)
が得られる．(4.6)は，`→∞においてx(`)がx∗に2次収束することを表している．
最後に，x(`+1)=Ψ(Φ(x(`)))(漸化式(2.3))によって生成される{x(`)}`=0,1,...の収束性に関する 定理を与える．
定理4.5 (z,x)6=0とする．x(0)がαkxkに十分近いならば，x(`+1)=Ψ(Φ(x(`)))によって生成さ
れる{x(`)}`=0,1,...は`→∞においてαkxkに2次収束する．
証明 Ψ(x+∆x)のテイラー級数を考える．(z,x+∆x) = (z,x) + (z,∆x)なので，
Ψ(x+∆x) = C (z,x)
(
1+(z,∆x) (z,x)
)−1
(x+∆x) となる．また，
(
1+(z,∆x) (z,x)
)−1
=1−
((z,∆x) (z,x)
) +
((z,∆x) (z,x)
)2
− ···
を利用すると，
Ψ(x+∆x) = C
(z,x)x+ C
(z,x)∆x−C(z,∆x)
(z,x)2 x−C(z,∆x)
(z,x)2 ∆x+C(z,∆x)2
(z,x)3 x+η(kxk3) (4.7) が得られる．ここで，∆x(`):=x(`)−x∗とおく．補助定理4.3においてx=x∗,∆x=∆x(`)とする と，F(x∗) =0,Φ(x∗) =x∗,Φ(x∗) =x∗より，
Φ(x(`)) =x∗−λ(∆x(`))(J(x∗))−1∆x(`)+η(kxk3) (4.8) となる．(4.7)においてx=x∗,∆x=−λ(∆x(`))(J(x∗))−1∆x(`)+η(kxk3)とし，(z,x∗) =Cおよび (4.8)を使うと，
Ψ(Φ(x(`))) =x∗−λ(∆x(`)) (
In−x∗zH C
)
(J(x∗)−1∆x(`)+η(kxk3) になる．よって，
kx(`+1)−x∗k=
J(x)−1λ(x(`)−x∗) (
In−x∗zH C
)
(x(`)−x∗) +η(kx(`)−x∗k3)
が成り立つ．ゆえに，kIn−x∗zH/Ck ≤M3となる正の定数M3が存在するので，定理4.4を踏まえ ると，
kx(`+1)−x∗k ≤M1M2M3kx(`)−x∗k2 (4.9)
が導かれる．(4.9)は`→∞においてx(`)がx∗に2次収束することを意味する．
5
5 まとめと今後の課題
本報告では，固有値分解のための超平面制約法について，いくつかの数値特性を明らかにした．
まず，超平面制約法の核となるニュートン型反復として5つのタイプが考えられるが，そのうち2 つのタイプが有効であることを示した．続いて，ニュートン型反復に含まれるヤコビ行列は常に 正則と仮定して，1つの有効なタイプに絞ってその局所2次収束性を証明した．
もう1つの有効なニュートン型反復について，今後その収束性を明らかにしたい．また，ヤコ ビ行列の正則性に関する議論や誤差解析による計算精度の検証なども今後の研究課題である．
参考文献
[1] J. W. Demmel: “Applied Numerical Linear Algebra”, SIAM, Philadelphia, 1997.
[2] K. Kondo, S. Yasukouchi and M. Iwasaki: “Eigendecomposition algorithms solving sequentially quadratic systems by Newton method”, JSIAM Letters1(2009), 40–43.
[3] LAPACK, http://www.netlib.org/lapack/.
[4] K.Yadani, K.Kondo and M.Iwasaki: “A singular value decompositon algorithm based on solving hyperplane constrained nonlinear systems”, Appl. Math. Comp.216(2010), 779–790.
[5] K. Yadani, K. Kondo and M. Iwasaki: “On the convergence of the V-type hyperplane constrained method for singular value decomposition”, JSIAM Letters2(2010), 21–24.
[6] K. Yadani, K. Kondo and M. Iwasaki: “Mixed double-multiple precision version of hyperplane constrained method for singular value decomposition”, JSIAM Letters2(2010), 25–28.
