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で計算したほうが速い場合も稀によくあります（誤記ではありません）．

ということで，みなさまよき固有値ライフを！