FMM(Fast Multipole Method)ライブラリの紹介 パート1

この記事の概略

• BIE (Boundary Integral Equation) Method の紹介
• fmm2d ($2$次元版Fast Multipole Methodライブラリ)の紹介
• $2$次元Laplace方程式の境界値問題をサンプル問題として，そのためのfmm2dに用意された関数を呼び出すラッパー関数をCとFortranで書いたので紹介(パート$2$で詳しく)

背景

Laplace方程式と呼ばれる$2$階線形楕円形偏微分方程式の境界値問題は，数学，物理，工学の観点から典型的で重要な問題として知られています．
本記事では，その問題を数値計算で近似的に解くための BIE (Boundary Integral Equation) Method と呼ばれる手法を紹介します．
数学的にはかなり大雑把にその枠組みに触れるだけにとどめ，後に紹介するfmm2d(Fortranライブラリ)がどのように役立つかを紹介します．

以下のように表される$2$次元の境界値問題を考え，その解を$u(x)$と表記します．

\begin{align}
\Delta u(x) &= 0, \, x = (x_1, x_2) \in \Omega,\\ 
u(x) &= f(x), \, x \in \Gamma.
\end{align}

ここで$\Delta = \frac{\partial^2}{\partial x_1^2} + \frac{\partial^2}{\partial x_2^2}$はLaplace作用素，$\Omega$は$2$次元空間内の有界な領域, $\Gamma = \partial \Omega$はその境界，$f(x)$は$\Gamma$上で定義された関数とします．($\Gamma$や$f(x)$の滑らかさ等に立ち入ると複雑になるので，"十分滑らか"と仮定して雑に進めます.)本記事では以下の図のような$\Omega$を単位円とする単純な例を考え，$\Gamma$上の$f(x)$の値を既知とします．

この問題の解(方程式と境界条件を満たす関数)$u(x)$が，$\Omega$内部の$x$でとる値を知りたいという状況です．

BIE (Boundary Integral Equation) Method

偏微分方程式の数値としては，FDM(有限差分法)，FEM(有限要素法)が有名どころですが，BIE method と呼ばれる手法も知られていて，FDM や FEM に比べると(数値解法としては)比較的新しい手法ではあるものの，近年活発に$2$次元，$3$次元の問題が研究されています．

BIE method は いわゆる Potential Theory に基づいており，標準的な資料として，Kress (2014, 参考文献1)が挙げられます．
まずLaplace方程式に対応する Green's function $G(x,y)$というものを導入します．

\begin{equation}
    G(x,y) = -\frac{1}{2 \pi} \ln |x-y|
\end{equation}

ここで$x$を$\Omega$内，$y$を$\Gamma$上の相異なる$2$点を表すベクトルとします．すると$G$の $y$における normal derivative $\frac{\partial G}{\partial n_y}(x,y)$を具体的に計算することが出来て，

\begin{equation}
    \frac{\partial G}{\partial n_y}(x,y)
    = \frac{1}{2\pi r^2}r \cdot \hat{n}_{y}
\end{equation}

と表されます．ここで$\hat{n}_{y}$は境界$\Gamma$上の$y$における外向き単位法線ベクトルです．

境界$\Gamma$や境界値を与える$f(x)$が何らかの意味である程度滑らかであれば，$\Gamma$上で定義されるある関数$\rho(y)$(密度関数と呼ばれます)が存在して，$\Omega$内の解$u(x)$は，$\frac{\partial G}{\partial n_y}(x,y)$と$\rho(y)$の convolution (畳み込み) と等しいことが知られています:

\begin{equation}
    u(x) = D[\rho](x) 
    = \int_{\Gamma}
    \frac{\partial G}{\partial n_y}(x,y) \rho(y) ds_y
\end{equation}

しかしながら$\rho$はまだ未知の関数です．
$u$を数値計算で近似的に求めるために，右辺の積分を，離散化した数値積分で十分良く近似したいので，少なくとも$\Gamma$上で離散化のために選んだ標本点における$\rho$の値が必要です．

上記$D[\rho](x)$の$x$は$\Omega$内の点ですが，$x$を$\Gamma$上のある点に，その点における境界に対する単位法線ベクトルと平行な向きで近づける極限をとると，以下の積分方程式(boundary integral equation)が(厳密には Principal Value, あるSchwartz distribution の意味で)成り立つことが知られています．

\begin{equation}
    \lim_{x \to y} u(x) = -\frac{1}{2}\rho(y) + \lim_{x \to y}D[\rho](x)
\end{equation}

左辺は境界値$f(y)$に収束します．ので，少々雑な表記をすると，

\begin{equation}
    f(x) = -\frac{1}{2}\rho(x) + D[\rho](x)
    = -\frac{1}{2}\rho(x)
      + \int_{\Gamma} \frac{\partial G}{\partial n_y}(x,y) \rho(y) ds_y
\end{equation}

が$\Gamma$に近づける極限の意味での各$x$で成り立つとみなすことが出来ます．

$D[\rho]$は積分で表せるので，数値計算のために，$\Gamma$上でサイズ$N_b$の標本点集合を選び，$(y_j)$ $(1 \le j \le N_b)$と表します．更に$(y_j)$と同一の標本点集合を用意し，異なる文字を用いて$(x_i)_{i=1}^{N_b}$を用意し，上記$x$と関連付けます．$(y_j)$の構成は，境界$\Gamma$を曲線でパラメータ表示し，標準区間[-1,1]上のGauss-Legendre数値積分の標本点を用いて設定する方法が標準的です．

$D[\rho]$の数値積分は，各$x_i$に対して，

\begin{equation}
    D[\rho](x_i)
    = \int_{\Gamma} \frac{\partial G}{\partial n_y}(x,y) \rho(y) ds_y
    \simeq \sum_{j=1}^{N_b} \frac{\partial G}{\partial n_{y_j}}(x_i,y_j) \rho(y_j) w_j
\end{equation}

と離散化して近似出来ます．$(w_j)_{j=1}^{N_b}$は上述のGauss-Legendre数値積分の標本点に対応する既知のウェイトと，曲線のパラメータ表示を標準区間[-1,1]上で表示するための変数変換によって現れる$1$階の導関数から既知の値として構成されます．

この近似により，未知の$\rho$を含む BIE は，各$x_i$に対して，

\begin{equation}
    \left[
        1 - 2\sum_{j=1}^{N_b} \frac{\partial G}{\partial n_{y_j}}(x_i,y_j) w_j
    \right]
    \rho(y_j)
    \simeq -2f(x_i)
\end{equation}

となり，近似的に解くために$\simeq$を$=$で置き換えて方程式を得ます．
全ての$i \ (1 \le i \le N_b)$についてそれらの方程式をまとめると，

\begin{equation}
    (I - K) \rho = b
\end{equation}

という形の線形方程式を得られます．
ここで$I$は$N_b \times N_b$単位行列，$K = (K_{ij}) = \left(\frac{\partial G}{\partial n_{y_j}}(x_i,y_j) w_j\right)$，$b = (b_i) = -2(f(x_i))$と表せます $(1 \le i, j \le N_b)$.

数値線形代数のアルゴリズムを用いて，この未知の$\rho$ベクトルの値を得られれば，その値を用いて$u(x) = D[\rho](x) = \int_{\Gamma} \frac{\partial G}{\partial n_y}(x,y) \rho(y) ds_y$の値を$\Omega$内の$x$においても同様に数値積分で近似的に計算することが出来ます．
という訳で，$\rho$のために上記線型方程式を解く必要がありますが，行列$K$は密行列で対称性も一般的には無く，反復法，典型的にはGMRES(参考文献2)がよく用いられます．GMRES等の反復法では，$K$と各反復ステップで得られる近似解ベクトルの行列-ベクトル積を計算します．$\rho$の十分な近似解を得るために必要な反復回数を$N_s$とおくと，上記の線型方程式を数値的に解くために必要な時間計算量は$O(N_b^2 N_s)$で，$N_s$の理論上の上限は$N_b$であるものの，実用的には$N_b$に比べて無視できるほど小さいことが知られています．
したがって，実用的な時間計算量は$O(N_b^2)$とみなすことが出来ます．
また，ここで紹介されていることと同様に， BIE method では$d$次元の領域における問題を，$d-1$次元の領域の境界上の密度関数を求めることで問題を近似的に解くため，
FDM や FEM といった他の代表的な手法に比べて問題を離散化するために必要な標本点の数を少なくすることが出来ます．($N_b$は FDM や FEM の離散化による標本点より次元が落ちる分少なくなる．)それでもなお，$O(N_b^2)$という密行列を伴う反復法の計算コストは FDM や FEM といった粗行列を扱う他の代表的な手法に比べて優れている訳では無く，時間・空間計算量の効率化が必要です．
ここで FMM (Fast Multipole Method) という近似アルゴリズムを適用することにより，$O(N_b^2)$ではなく$O(N_b)$で線型方程式を解くことが出来ます．

FMM は一言で大雑把にまとめてしまうと，Green's function やその方向微分，Hessian によって構成される$N \times M$行列と，密度関数のような$M \times 1$ベクトルの積(行列-ベクトル積)を，$O(NM)$ではなく$O(N + M)$で近似的に求めるアルゴリズムです．
FMM の仕組みについては本記事では深入りしませんが，$2$次元，$3$次元，対象となる偏微分方程式に対応する Green's function によって layer potentials が異なり，かなり多くの亜種が存在し，それに伴い論文や紹介記事も存在します．$2$次元の問題を対象とした最も古い文献は Greengard and Rokhlin (1987, 参考文献3)で，Chen (2015, 参考文献4)は入門の参考資料として良いかもしれません．

さらに，先述の線型方程式の$\rho$ベクトルを求めた後，最終目標である$u(x)$の$\Omega$内の近似解を求める段階でも，

\begin{align}
    u(x) &= D[\rho](x) 
    = \int_{\Gamma}
    \frac{\partial G}{\partial n_y}(x,y) \rho(y) ds_y \\
    &\simeq \sum_{j=1}^{N_b} \frac{\partial G}{\partial n_{y_j}}(x,y_j) \rho(y_j) w_j 
\end{align}

と同様な数値積分が先述の線型方程式の$1$行と似た形のベクトルと，$\rho$ベクトルの内積で得られます．領域$\Omega$内で計算したい$x$が$N_e$箇所あるとすると，全ての$(x_i)_{i=1}^{N_e}$についてこれらの内積をまとめて，線型方程式を解いたときと同様に FMM を適用出来ます．つまり，FMM によって時間計算量を$O(N_e N_b)$から$O(N_e + N_b)$に削減出来ます．

以下の２つの図は 単位円内のサンプル問題の数値解を求めるための

1. 線型方程式の$\rho$を求めるために要した時間
2. $\rho$の値を用いて数値解$u(x)$を求めるために要した時間
   をそれぞれ異なるインプットデータサイズで計算し，結果を描画したものです．
   詳細については本記事の続きのパート$2$で説明します．

これらの数値実験のために，以下で紹介する fmm2d という$2$次元FMMのライブラリを用いました．

fmm2d: ２次元問題用Fortranライブラリ

fmm2d は以下のリポジトリで公開されています．

本記事で具体的に使用したライブラリ内の関数は，
fmm2d/src/laplace/lfmm2dwrap.f
で実装されている

1. 'lfmm2d_s_d_p'
2. 'lfmm2d_t_d_p'
   の$2$つです．１つ目を BIE の線型方程式を解くために，２つ目を$u(x)$の数値解を計算するために用いました．また，　データ入力の都合上，それぞれに以下のラッパー関数を実装しました．

fmm2d_laplace_wrap.f90

module lfmm2d_wrap
  use iso_c_binding
  use, intrinsic :: iso_fortran_env
  implicit none
  contains

  subroutine lfmm2d_t_d_p_wrap(eps_fmm, n_src, src_in, dip_str_in, dip_vec_in, &
    n_trg, trg_in, pot_trg_out, flg_err) bind(c)
    implicit none
    real(c_double), value :: eps_fmm
    integer(c_int), value :: n_src, n_trg
    integer(c_int), value :: flg_err
    ! sizes of *_in below: n_src*2 all
    real(c_double), intent(in) :: src_in(*), dip_str_in(*), dip_vec_in(*)
    ! ((i-1)*2 + j) -> (j, i) (size (n*2,1) -> (2, n))
    double precision, allocatable :: src(:,:)
    double complex, allocatable :: dip_str(:)
    double precision, allocatable :: dip_vec(:,:)
    ! sizes of *_in below: n_trg*2 all
    real(c_double), intent(in) :: trg_in(*)
    real(c_double), intent(inout) :: pot_trg_out(*)
    ! (i-1)*2 + j) -> (j, i) (size (n*2,1) -> (2, n))
    double precision, allocatable :: trg(:,:)
    double complex, allocatable :: pot_trg(:)
    integer :: i, j
    double complex :: iu

    iu = (0.0d0, 1.0d0)

    allocate(src(2, n_src), dip_str(n_src), dip_vec(2, n_src))
    allocate(trg(2, n_trg), pot_trg(n_trg))

    do i = 1, n_src
      do j = 1, 2
        src(j, i) = src_in((i-1)*2 + j)
        dip_vec(j, i) = dip_vec_in((i-1)*2 + j)
      end do
      dip_str(i) = dip_str_in((i-1)*2 + 1) + dip_str_in((i-1)*2 + 2)*iu
    end do
    do i = 1, n_trg
      do j = 1, 2
        trg(j, i) = trg_in((i-1)*2 + j)
      end do
    end do
    
    flg_err = 0
    call lfmm2d_t_d_p(eps_fmm, n_src, src, &
      dip_str, dip_vec, &
      n_trg, trg, pot_trg, flg_err)
    do i = 1, n_trg
      pot_trg_out((i-1)*2 + 1) = real(pot_trg(i))
      pot_trg_out((i-1)*2 + 2) = imag(pot_trg(i))
    end do

    deallocate(src, dip_str, dip_vec)
    deallocate(trg, pot_trg)
  end subroutine lfmm2d_t_d_p_wrap

  subroutine lfmm2d_s_d_p_wrap(eps_fmm, n_src, src_in, dip_str_in, dip_vec_in, &
     pot_out, flg_err) bind(c)
    implicit none
    real(c_double), value :: eps_fmm
    integer(c_int), value :: n_src
    integer(c_int), value :: flg_err
    ! sizes of *_in/out below: n_src*2 all
    real(c_double), intent(in) :: src_in(*), dip_str_in(*), dip_vec_in(*)
    real(c_double), intent(inout) :: pot_out(*)
    ! ((i-1)*2 + j) -> (j, i) (size (n*2,1) -> (2, n))
    double precision, allocatable :: src(:,:)
    double complex, allocatable :: dip_str(:)
    double precision, allocatable :: dip_vec(:,:)
    double complex, allocatable :: pot_src(:)
    integer :: i, j
    double complex :: iu

    iu = (0.0d0, 1.0d0)

    allocate(src(2, n_src), dip_str(n_src), dip_vec(2, n_src), pot_src(n_src))

    do i = 1, n_src
      do j = 1, 2
        src(j, i) = src_in((i-1)*2 + j)
        dip_vec(j, i) = dip_vec_in((i-1)*2 + j)
      end do
      dip_str(i) = dip_str_in((i-1)*2 + 1) + dip_str_in((i-1)*2 + 2)*iu
    end do

    flg_err = 0
    call lfmm2d_s_d_p(eps_fmm, n_src, src, dip_str, dip_vec, pot_src, flg_err)

    do i = 1, n_src
      pot_out((i-1)*2 + 1) = real(pot_src(i))
      pot_out((i-1)*2 + 2) = imag(pot_src(i))
    end do

    deallocate(src, dip_str, dip_vec, pot_src)
  end subroutine lfmm2d_s_d_p_wrap
end module lfmm2d_wrap

fmm2d_laplace_wrap.h

void lfmm2d_t_d_p_wrap(double eps_fmm, int n_src, double *src_in,
  double *dip_str_in, double *dip_vec_in,
  int n_trg, double *trg_in, double *pot_trg_out,
  int flg_err);

void lfmm2d_s_d_p_wrap(double eps_fmm, int n_src, double *src_in,
  double *dip_str_in, double *dip_vec_in,
  double *pot_out,
  int flg_err);

これらのコードと行った数値実験について，パート$2$で詳しく紹介する予定です．

参考文献

1. R. Kress, Linear Integral Equations, Third Edition, Springer, 2014.
2. Y. Saad and M.H. Schultz, GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems, J. Sci. Statist. Comput., 7 (1986), pp. 856-869.
3. L. Greengard and V. Rokhlin, A fast algorithm for particle simulations, J. Comput. Phys., 73(2), 1987
4. L. Chen, Introduction To Fast Multipole Methods, https://www.math.uci.edu/~chenlong/MathPKU/FMMsimple.pdf