FFTWを使った畳み込み積分の高速化

目的

無限区間で定義された次の畳み込み積分

を数値的に計算することを考えよう。関数 f(x), g(x) は x=0 の近くに局在している関数として、計算区間としては対称的な区間 -W/2<x<W/2 を取りたいとする。
最も単純な方法として、各 y に対して上の数値積分をしてしまうと、分割点数 N に対して O(N²) の計算コストがかかってしまい、N が大きくなると計算のネックとなりかねない。 そこで高速フーリエ変換を利用して上の畳み込みを実行することで、計算コストのオーダーを軽減したい。FFTの計算コストはO(N log N)、また変換後の畳み込みは（単に N 回のかけ算になり）O(N)、全体で O(N log N) に軽減される。
高速フーリエ変換ライブラリFFTW3を利用することにする。

FFTができること

次のN個の複素数のデータ

に対して、FFT は次の N 個の複素数 F_t を高速に計算してくれる。

つまり、離散フーリエ変換に特化したアルゴリズムとなっている。畳み込み積分は（連続）フーリエ変換によって、もとの関数の単純な積の形にする事が出来ることはよく知られているが、離散化した場合も"ほぼ"同様のことがいえる。ただ連続積分を離散化する時の対応関係を正しく理解していないと、位相がずれて見当違いな結果が出てしまうので注意したい。

畳み込み積分の変形

先の畳み込み積分において、まず有限区間 [-W/2,W/2] に制限して離散化する。

次に畳み込み積分を離散和で近似する。

ここで、先の離散化の約束に従うと

となるから、次が得られる。

ここで g の引数に N/2 の項が出てきたが、これは区間の始点を -W/2 に取ったから生じたもので忘れてはならない。もし区間が [0,W] であればこのような項は生じない。 さて上式の i について離散フーリエ変換を行うと、単純な計算ののち、FFT版の畳み込み積分の公式が得られる。

FFT版畳み込み積分の公式と実装

先の公式を利用して実装したものをgithubにアップロードした。 FFTライブラリを利用したコアの部分を抜粋すると次の通り。特に(-1)^tの位相因子と規格化因子に注意。

積分幅依存性と注意点

例として、以下２つのケースを例にとって計算を実行してみよう。

例１：y 軸対称な配置での畳み込み

の二つを畳み込んだ結果が以下になる。区間幅 W=20,40,80 に対してプロットしてあり、上が全体の様子、下が厳密解との差分。

狭い区間でも十分な結果が得られているように見える。実は対称な配置では、周期 W の関数と見たときに両端が偶然連続的に繋がるため、たまたま良い結果が得られるだけである。

例１：y 軸非対称な配置での畳み込み

では次にx軸に非対称な配置でやってみよう。先ほどとの大きな違いは、厳密な答えが区間の両端で異なる値を取っている点であり、周期関数としては不連続点が存在してしまう。 グラフを見れば明らかなように、区間端点での不連続性を反映して、その付近での誤差がかなり大きくなる。 このような場合は、両端での値の差が必要な精度より小さくなるよう W を大きく取らなければならない。