Objective and Motivation
構造力学の数値解析として，モード解析またはモーダル解析と呼ばれる手法があります．これは，半離散化された運動方程式
$$ M \frac{\partial^2 u_I}{\partial t^2} + K u_I =0$$
を位相と振幅に分けて共通項を消去すると，
$$ K U_I = \omega^2 M U_I $$
の形の一般化固有値問題を解くことで，自由度Nに対して最大N次の固有振動数と固有振動モードを求めることができます．
このとき，固有値と固有ベクトルを求める必要があります．固有値をデジタルコンピュータで求めるには，
- QR法
- Arnoldi法
- Locally Optimal Block Preconditioned Conjugate Gradient(LOBPCG)法
等の方法が知られています．有名なLAPACKはQR法を用いており，実行時に自由度$N$のとき$O(N^2)$のメモリを確保せねばならず，条件分岐を伴う複雑なアルゴリズムを実行する必要があり，そのため大規模問題(例えば$N>100,000$)への適用には工夫が必要です．
一方で，計算コストを下げられる手法として，Krylov部分空間法を用いたアルゴリズムが知られます．上記のArnoldi法(Lanczos法)やLOBPCG法はKrylov部分空間法ベースの手法であり，特にLOBPCG法は求める固有値の数$k$(最大固有値または最小固有値から数えて$k$番目まで)の数が全固有値の20-30%を超えない範囲であれば適用できるマトリックスフリーな方法であることが知られています(https://github.com/lobpcg/blopex)．
モード解析においては，最小固有値から高々5個ないし10個程度の固有振動数と固有振動モードが計算できれば十分であることも多く，自由度Nに対してkは数%程度となることが通常ですので，LOBCPG法の弱点であるkの値の上限は全く問題になりません．一方で，マトリクス$M$や$K$は一般に大規模な疎行列であることが多いので，LOBPCG法の強みである，主にSpMVだけで計算を進められるメリットが生きます．
そこで，本記事では，一般化固有値問題を対象としたLOBPCGの実装方法をレビューするとともに，ソルバーを運用していて体験した本手法の強み/弱みを付記します．
その前に，そもそもLOBPCGをユーザである私が実装する必要があるか?
あなたにとっては，その必要はありません．各種言語で著者らの手による実装が公開されており，それを直接的ないし間接的に呼び出して使えばよいです．
まずはドキュメントを読み，それでわからなければコードを読めば何をやっているかわかるので，実装する必要はありません．もしそれでもわからなければ以下の文献をご一読ください．
A. V. Knyazev, I. Lashuk, M. E. Argentati, and E. Ovchinnikov, Block Locally Optimal Preconditioned Eigenvalue Xolvers (BLOPEX) in hypre and PETSc, SIAM Journal on Scientific Computing 25(5): 2224-2239 (2007) DOI
A. V. Knyazev, Toward the Optimal Preconditioned Eigensolver: Locally Optimal Block Preconditioned Conjugate Gradient Method, SIAM Journal on Scientific Computing 23(2): 517-541 (2001) DOI
ここまで読んで頂きありがとうございました．記事は以上で終わります．
LOBPCG法の実装
Single-Vector LOBPCG法の実装
では，部外者が帰ったようですので，実装に移ります．
上記Githubリンクでは2報の論文が紹介されていますが，そのうち
Knyazev, A. V., Argentati, M. E., Lashuk, I. & Ovtchinnikov, E. E. Block Locally Optimal Preconditioned Eigenvalue Xolvers (BLOPEX) in hypre and PETSc. SIAM J. Sci. Comput. 29, 2224–2239 (2007).
の論文がわかりやすくてよいです．
論文にしたがって，まずは最小固有値を求める場合のみ考えてみます．
LOBPCG法は，ざっくりいえばReyleigh商を最小化ないし最大化する固有値・固有ベクトルを求める最適化問題を，適当な初期推定値を改善することで求める方法です．
本文p. 3のSection 2.1.2に，最小固有値のみを求めるアルゴリズムが僅か2行で紹介されています．
実装としては，
(0) 初期推定値として$x^{(0)}$を与えます．これは一様乱数列を与えてもよいですが，すでにおおよその予測値があればそれを用いることが推奨されます．モード解析であれば，たとえばよりシンプルな串団子のモデルで事前解析を回しておいて，その結果をもとに最小固有値に対応する固有モードを初期推定値として与えるなどすると，局所解に寄り道しにくく，収束が非常に速くなります．
(1) 3項間漸化式(式2.1の1行目)により固有ベクトル$x$を更新します．
(2) 前処理行列$T$，なければ$T=I$を使って$w$ベクトルを更新します．以降，収束するまで(1)-(2)を繰り返します．
実装例はこのようになります．
subroutine LOBPCG_single_CRS(A,B,lambda,X,alpha,tol,debug)
type(CRS_),intent(in) :: A, B
real(real64),allocatable,intent(inout) :: X(:)
real(real64),allocatable,intent(inout) :: lambda
real(real64),intent(in) :: alpha
logical,optional,intent(in) :: debug
logical :: debug_mode
type(Random_) :: random
real(real64),allocatable :: r(:),rho
integer(int32) :: n,i
integer(int32) :: MAX_ITR = 1000000
real(real64),intent(in) :: tol
debug_mode = input(default=.false.,option=debug)
! BLOPEX IN HYPRE AND PETSC, SIAM Eq. (2.2)
! number of eigen value :: m
! initialize X and lambda
n = A%size()
X = 2.0d0*random%randn(n)-1.0d0
lambda = 0.0d0
!Single-vector version
do i=1,MAX_ITR
rho = dot_product(x, A%matmul(x) )/dot_product(x,B%matmul(x))
r = A%matmul(x) - rho*B%matmul(x)
x = x - alpha*r
if(debug_mode)then
print *, i, norm(r)
endif
if(norm(r) < tol )then
if(debug_mode)then
print *, "[OK] converged."
endif
exit
else
cycle
endif
enddo
if(i==MAX_ITR) print *, "[ERROR] LOBPCG NOT CONVERGED."
lambda = rho
end subroutine
なお，type :: CRS_ は適当なCRS形式のデータ型です．
前処理行列が単位行列でも，悪条件でない行列であれば素直に収束します(個人の感想です)．
複数の固有ベクトル/固有値の組を求める場合
この場合，Section 2.2のアルゴリズムによることになります．
反復すべきことはそう難しくはなく，
(1) m本のn次元固有ベクトル候補からなる行列X, n×m残差行列R，それらと直交するベクトルからなるn×m行列Pの3つを作成．
(2) V = [X, R, P] によりVを求める．ただし[,]は列方向にスタックするオペレータとする． また，Vの全ての縦ベクトルを直交化しておく．
(3) $\bar{A} = V^T A V$,$\bar{B} = V^T B V$により$\bar{A}$と$\bar{B}$を求め，一般化固有値問題$\bar{B}x = \lambda \bar{A}x$をLAPACKなどにより解く
(3) 求めた3m個の固有値$\lambda$を用いて，残差行列ベクトルRを更新する．
(4) $P=V \bar{B}$などにより行列Pを更新する．
(5) 残差が十分小さくなるまで(1)-(4)を繰り返す．
となります．
ただし，ここでも初期Xおよび初期Pの決定の仕方が問題となります．そこは永井さんにより公開されている以下の記事が非常に丁寧です．
要は，一旦はV = [X,R]で(1)-(4)を行い，RとPを求めた上で上記のV=[X,R,P]の反復を実行することになります．
以下がその実装例です．
subroutine LOBPCG_CRS(A,B,lambda,X,m,MAX_ITR,TOL,debug)
type(CRS_),intent(in) :: A, B
real(real64),allocatable,intent(out) :: X(:,:)
real(real64),allocatable,intent(out) :: lambda(:)
real(real64),intent(in) :: TOL
integer(int32),intent(in) :: m,MAX_ITR
logical,optional,intent(in) :: debug
logical :: debug_mode
type(Random_) :: random
real(real64),allocatable :: r(:,:), p(:,:),V(:,:),A_small(:,:),AV(:,:),&
b_small(:,:),lambda_small(:), XAX(:,:),AX(:,:),OR(:,:),BX(:,:),XBX(:,:),&
Bm_small(:,:),BV(:,:)
real(real64) :: initial_R
integer(int32) :: n,i,j,k
!integer(int32) :: MAX_ITR = 1000000
! References:
! Knyazev, A. V., Argentati, M. E., Lashuk, I. & Ovtchinnikov, E. E. Block Locally Optimal Preconditioned Eigenvalue Xolvers (BLOPEX) in hypre and PETSc. SIAM J. Sci. Comput. 29, 2224–2239 (2007).
! https://en.wikipedia.org/wiki/LOBPCG
! https://webpark1378.sakura.ne.jp/nagai/note_141009_LOBPCG.pdf
debug_mode = input(default=.false.,option=debug)
! BLOPEX IN HYPRE AND PETSC, SIAM Eq. (2.2)
! number of eigen value :: m
! debug
! initialize X and lambda
n = A%size()
!X = random%randn(n,m)
X = random%randn(n,m)! + random%randn(n,m) !+ random%randn(n,m)
do i=1,m
X(:,i) = X(:,i) / norm(X(:,i) )
enddo
call GramSchmidt(X,size(X,1),size(X,2),X )
lambda = zeros(m)
R = zeros(n,m)
P = zeros(n,m)
AV = zeros(n,3*m)
AX = zeros(n,m)
A_small = zeros(3*m,3*m)
V = zeros(n,3*m)
! m x m 固有値問題
do i=1,m
AX(:,i) = A%matmul(X(:,i) )
enddo
BX = zeros(size(AX,1),size(AX,2) )
do i=1,m
BX(:,i) = B%matmul(X(:,i) )
enddo
XAX = matmul(transpose(X(:,1:m) ),AX)
XBX = matmul(transpose(X(:,1:m) ),BX)
call LAPACK_EIG(XAX, XBX, b_small,lambda_small)
X = matmul(X,b_small)
do i=1,m
AX(:,i) = A%matmul(X(:,i) )
enddo
BX = zeros(size(AX,1),size(AX,2) )
do i=1,m
BX(:,i) = B%matmul(X(:,i) )
enddo
do i=1,m
R(:,i) = AX(:,i) - BX(:,i)*lambda_small(i)
enddo
V = zeros(n,2*m)
V(:,1 : m ) = X(:,:)
V(:,m+1 : 2*m ) = R(:,:)
call GramSchmidt(V,size(V,1),size(V,2),V )
do k=1,size(X,2)
V(:,k) = V(:,k) / norm(V(:,k) )
enddo
X(:,:) = V(:,1 : m )
R(:,:) = V(:,m+1 : 2*m )
AV = zeros(n,2*m)
do j=1,2*m ! n x n, n x 2m
AV(:,j)= A%matmul(V(:,j) )
enddo
BV = zeros(n,2*m)
do j=1,2*m ! n x n, n x 2m
BV(:,j)= B%matmul(V(:,j) )
enddo
A_small = matmul(transpose(V),AV)
Bm_small = matmul(transpose(V),BV)
call LAPACK_EIG(A_small,Bm_small, b_small,lambda_small)
do j=1,m
X(:,j) = matmul(V,b_small(:,j) )
enddo
AX = zeros(n,m)
do i=1,m
AX(:,i) = A%matmul(X(:,i) )
enddo
BX = zeros(size(AX,1),size(AX,2) )
do j=1,m
BX(:,j) = B%matmul(X(:,j) )
enddo
R = zeros(n,m)
do j=1,m
R(:,j) = AX(:,j) - BX(:,j)*lambda_small(j)
enddo
OR = zeros(n,2*m)
OR(:,1:m) = 0.0d0
OR(:,m+1:2*m) = R(:,:)
P = zeros(n,m)
do i=1,m
P(:,i) = matmul(OR,b_small(:,i) )
enddo
V = zeros(n,3*m)
AV = zeros(n,3*m)
BV = zeros(n,3*m)
AX = zeros(n,m)
BX = zeros(n,m)
do i=1,MAX_ITR
V(:,1 : m ) = X(:,:)
V(:,m+1 : 2*m ) = R(:,:)
V(:,2*m+1 : 3*m ) = P(:,:)
! Gram-Scmidtを計算する．
call GramSchmidt(V,size(V,1),size(V,2),V )
do k=1,size(X,2)
V(:,k) = V(:,k) / norm(V(:,k) )
enddo
do k=size(X,2)+size(R,2)+1,size(V,2)
V(:,k) = V(:,k) / norm(V(:,k) )
enddo
X(:,:) = V(:,1 : m )
R(:,:) = V(:,m+1 : 2*m )
P(:,:) = V(:,2*m+1 : 3*m )
!$OMP parallel do
do j=1,3*m
AV(:,j)=A%matmul(V(:,j) )
enddo
!$OMP end parallel do
!$OMP parallel do
do j=1,3*m
BV(:,j)=B%matmul(V(:,j) )
enddo
!$OMP end parallel do
A_small = matmul(transpose(V),AV)
Bm_small = matmul(transpose(V),BV)
! LAPACKの固有値求めるやつを呼ぶ
call LAPACK_EIG(A_small,Bm_small,b_small,lambda_small)
X = matmul(V,b_small(:,1:m) )
!$OMP parallel do
do j=1,m
AX(:,j) = A%matmul(X(:,j) )
enddo
!$OMP end parallel do
!$OMP parallel do
do j=1,m
BX(:,j) = B%matmul(X(:,j) )
enddo
!$OMP end parallel do
!$OMP parallel do
do j=1,m
R(:,j) = AX(:,j) - BX(:,j)*lambda_small(j)
enddo
!$OMP end parallel do
V(:,1 : m ) = 0.0d0
! matmul:: (n,3*m) x (n,3*m)
P = matmul(V,b_small(:,1:m) )
! Detect convergence
if(i==1)then
initial_R = maxval(abs(R) )
endif
if(debug_mode)then
print *, i, maxval(abs(R) )/initial_R
endif
if(maxval(abs(R) )/initial_R < tol )then
if(debug_mode)then
print *, "[OK] converged."
endif
lambda = lambda_small(1:m)
exit
else
cycle
endif
enddo
if(i==MAX_ITR) print *, "[ERROR] LOBPCG NOT CONVERGED."
end subroutine
! #################################################
上記コードを運用してみての感想
実際に計算を行うと，自由度$N>1,000,000$程度の悪条件でない疎行列であれば，ごかてい（なぜか変換できない)のパソコンであっという間に計算できます．また，ここでは詳細には記載しませんでしたが，SpMVを並列化すればアクセラレータやクラスタを用いて高速化することも容易で，役に立っています．
加えて，許容誤差を緩めた計算を回してからよりシビアな許容誤差を与えることで，固有モードが分離されていく様子が確認でき，解釈性にも優れると感じています．
ただし，悪条件行列であれば，当然のこと著しく収束性が悪化します．また，このコードのように前処理行列を用いない(これを用いないLOBPCGというのもなんか矛盾しますが)場合には，先述の本家実装と比べて収束性が悪い感じもします（再び，個人の感想です）．
また，解の近傍になるほど収束速度が著しく遅くなるので，その点も注意が必要です．自由度$N<10,000$程度の小規模な問題の場合，LAPACKで計算したほうが速い場合も稀によくあります（誤記ではありません）．
ということで，みなさまよき固有値ライフを！