Conv. check
Simple: 100
TR : 90
1.0E-16 1.0E-14 1.0E-12 1.0E-10 1.0E-08 1.0E-06 1.0E-04 1.0E-02 1.0E+00
0.0 0.2 0.4 0.6 0.8 1.0 1.2
演算量[Gflop]
残差
性能評価 ( Ⅴ ) 〜 Lanczos 法 収束性 (2)
θ
28θ
30θ
29TR
simple
• 密集固有値が ある 場合： bcsstk21
• TR-Lanczos は密集固有値に有効
Conv. check
Simple: 20
TR : 70
５ . Jacobi-Davidson 法
• Jacobi-Davidson 法を考える動機
• Jacobi-Davidson 法の原理
• アルゴリズム
• リスタート
• デフレーション
• 性能評価
Jacobi-Davidson 法を考える動機
• Lanczos 法／ Arnoldi 法で解きにくい問題
–
スペクトルの端にある固有値で，かつ他の固有値からの分離が 良い場合は解きやすい。–
それ以外の場合，効率的に計算を行うには，スペクトル変換など の技法が必要–
しかし，(A –
σM)
–1 による乗算は，コストが掛かる場合が多い。スペクトル変換を用いずに，内部固有値／分離の悪い固有値を 効率的に求められる解法が必要
Jacobi-Davidson 法の原理 (1)
• 修正ベクトルと部分空間の拡張
–
部分空間S
m と正規直交基底{v
1, v
2, … , v
m}
とが与えられてい るとする。V
m= [v
1, v
2, … , v
m]
とし，S
mに関するA
のRitz
値，Ritz
ベクトルの組を(
θj, u
j= V
ms
j)
( j = 1, … , m)
とする。–
いま，求めたい固有値λに対する近似固有ベクトルとしてu
j を 選び，u
j を正しい固有ベクトルにするための修正ベクトルをt
と 書くと，– Jacobi-Davidson
法 では，ベクトルu
j の直交補空間の中でt
を求 め，これによってS
m を拡張する。
A ( u
j+ t ) =
λ( u
j+ t )
すなわち，(A –
λI ) t = – (A –
λI ) u
jJacobi-Davidson 法の原理 (2)
• 修正方程式
– (A –
λI )
の作用をu
j⊥ に制限した演算子は( I – u
ju
jt)(A –
λI )( I – u
ju
jt)
だから，
t
を求めるための方程式はとなる。さらに，λは未知だから θj で置き換えると，
これを修正方程式と呼ぶ。
– JD
法では，t
を求めた後，これを{v
1, v
2, … , v
m}
に対して正規 直交化してv
m+1 を求め，部分空間を拡大していく。( I – u
ju
jt)(A –
λI )( I – u
ju
jt) t = – (A –
λI ) u
j( I – u
ju
jt)(A –
θjI )( I – u
ju
jt) t = – (A –
θjI ) u
jJacobi-Davidson 法の原理 (3)
• 修正方程式の解法
–
修正方程式の係数行列はu
j⊥ への制限のため特異であるが，右辺も
u
j⊥ に属する。したがって方程式は無矛盾–
初期値t
0= 0
から出発し，GMRES
法，CGS
法などのKrylov
部 分空間法を用いて解けば，自動的にu
j⊥ に属する解が得られる。–
経験上，厳密解を求める必要はなく，GMRES
法の５ステップ分 程度の精度で十分–
前処理を行う場合は前処理行列K
も部分空間u
j⊥ に制限する 必要がある。アルゴリズム
58初期ベクトル
t = v
0を設定DO m = 1, 2, …
DO i = 1, 2, …, m–1
t := t – (t
tv
i) v
iModified Gram-Schmidt
法END DO
v
m= t /
‖t
‖2,
v
mA= Av
mDO i = 1, 2, …, m
M
i,m= v
itv
mAEND DO
m
×m
行列M
の最大固有値θ，固有ベクトルs
を求める。
u = Vs
（V = [ v
1, v
2, … , v
m]
）u
A= V
As
（V
A= A[v
1, v
2, … , v
m]
）r = u
A–
θu
‖
r
‖2 ≦εならλ=
θ，x = u
として停止
( I – u u
t)(A –
θI )( I – u u
t) t = –r
を解いてt
を求める。END DO
リスタート
• 必要性
–
１ステップ当たりのメモリ所要量・演算量を一定値に抑えるため，JD
法でもリスタートを行う。• 実現方法
– JD
法では，部分空間S
mがKrylov
部分空間の構造を持つ必要 がないため，リスタートは容易– Ritz
ベクトルのうち，求めたい固有値に近いRitz
値を持つベクトルを複数本選び，部分空間の次元を削減して計算を続行すれば よい。
デフレーション
60• 必要性
–
複数の固有値・固有ベクトルを求める場合，収束した固有ベクト ルを基底から取り除いてリスタートすると，計算が効率化できる。–
この操作をデフレーションと呼ぶ。• 実現方法
–
収束した固有ベクトルx
1, … , x
lを基底から取り除く。–
アルゴリズムにおいて，A –
θI
をで置き換える。
–
これにより，計算はすべてX
l⊥の中で行える。( I – X
lX
lt)(A –
θI )( I – X
lX
lt)
ただしX
l= [x
1, … , x
l]
性能評価 〜 JD 法 対称・収束特性
cf) Sleijpen, Jacobi-Davidson algorithms for various eigenproblems(1999),p.15 Fig.4-1 Iteration number
Residual norm
Iteration number
Residual norm
• テスト行列 N=1000, A
ii=i, A
i-1,i=0.5, A
1000,1=0.5
• 収束判定 ‖ Ax – x θ‖ < 10
-8Arnoldi JD with GMRES(5)
疎行列固有値計算全般
• Z. Bai, J. Demmel, J. Dongarra, A.Ruhe and H. Van der Vorst (eds.):
“Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide”, SIAM, Philadelphia, 2000.
• F. Chatlin: “Valeurs Propres de Matrices”, Masson, Paris, 1988. （W.
Ledermann “Eigenvalues of Matrices”, John-Wiley and Sons, Chichester, 1993.）
• J. W. Demmel: “Numerical Linear Algebra”, SIAM, Philadelphia, 1996.
• G. H. Golub and C. F. van Loan: “Matrix Computations”, 3rd edition, The Johns Hopkins University Press, 1989.
• B. N. Parlett: “The Symmetric Eigenvalue Problem”, Prentice-Hall, Englewood Cliffs, 1980.
• Y. Saad: “Numerical Methods for Large Eigenvalue Problems”, Halsted Press, New York, 1992.
• J. H. Wilkinson: “The Algebraic Eigenvalue Problem”, Claredon Press, Oxford, 1965.
参考文献（１ /5 ）
参考文献 (2/5)
Lanczos
法(1)
• J. Cullum and W. E. Donath: “A Block Lanczos Algorithm for Computing the q
Algebraically Largest Eigenvalues and a Corresponding Eigenspace of Large Sparse Real Symmetric Matrices”, Proc. of the 1974 IEEE Conference on Decision and
Control, Phoenix, Arizona, pp. 505-509 (1974).
• J. Cullum and R. A. Willoughby: “Lanczos Algorithms for Large Symmetric Eigenvalue Computations”, Volume 1, Theory, Birkhauser, Boston, 1985.
• R. G. Grimes, J. G. Lewis and H. D. Simon: “A Shifted Block Lanczos Algorithm for Solving Sparse Symmetric Generalized Eigenproblems”, SIAM Journal on Matrix Analysis and Applications, Vol. 15, No. 1, pp. 228-272 (1994).
• C. Lanczos: “An Iteration Method for the Solution of the Eigenvalue Problem of Linear Differential and Integral Operators”, J. Research of the National Bureau of Standards, Vol. 45, No. 4, pp. 255-282 (1950).
• 森正武, 杉原正顕, 室田一雄: “線形計算”, 岩波書店, 1994.
参考文献 (3/5)
Lanczos 法 (2)
• B. N. Parlett and D. Scott: “The Lanczos Algorithm with Selective
Orthogonalization”, Mathematics of Computation, Vol. 33, pp. 217-238 (1979).
• H. D. Simon: “The Lanczos Algorithm with Partial Reorthogonalization”, Mathematics of Computation, Vol. 42, pp. 115-136 (1984).
• R. Underwood: “An Iterative Block Lanczos Method for the Solution of Large Sparse Symmetric Eigenproblems”, Report STAN-CS-75-495, Department of Computer Science, Stanford University, Stanford, California (1975).
• K. Wu and H. Simon: “Thick-Restart Lanczos Method for Large Symmetric
Eigenvalue Problems”, SIAM Journal on Matrix Analysis and Applications, Vol. 22, No. 2, pp. 602-616 (2000).
参考文献 (4/5)
Arnoldi
法• W. E. Arnoldi: “The Principle of Minimized Iterations in the Solution of the Matrix Eigenvalue Problem”, Quart. J. Applied Mathematics, Vol. 9, pp. 17-29 (1951).
• R. B. Lehoucq and D. C. Sorensen: “Deflation Techniques for an Implicitly
Restarted Arnoldi Iteration”, SIAM Journal on Matrix Analysis and Applications, Vol. 17, No. 4, pp. 789-821 (1996).
• R. B. Lehoucq, K.J. Maschhoff.: “implementation of an implicitly restarted block Arnoldi method”, Preprint MCS-P649-0297, Argonne National Laboratory,IL,1997.
• R. B. Lehoucq, D. C. Sorensen and C. Yang: “ARPACK User’s Guide”, SIAM, Philadelphia, 1998.
• D. C. Sorensen: “Implicit Application of Polynomial Filters in a k-Step Arnoldi Method”, SIAM Journal on Matrix Analysis and Applications, Vol. 13, No. 1, pp.
357-385 (1992).