FFT アルゴリズムの探索メソッド比較 （Estimate/Measure/Patient）

# やったこと

`fft` 関数[^1]を実行すると、まず FFT に使用するアルゴリズムを決定する処理が実行されます。

```matlab
method = fftw('planner')
```

で表示されるアレです。場合によっては結構時間を食うこのステップをスキップする方法について紹介します。

計算処理への影響は、ほとんどの場合**極小**かと思いますが、単発で大きな配列に `fft` を実行する場合は、知っておいて損はしないかもしれません。

[^1]: see https://jp.mathworks.com/help/matlab/ref/fft.html


### 環境

MATLAB R2019b


# FFTW と planner と `fftw`

`fft` は使ったことがある方も多いと思いますが、これはもちろん、あの FFTW (the Fastest Fourier Transform in the West) に「基づいた」実装になっています[^1]。FFTW は以下の引用にある通り、ハードウェアの状況に合わせて最適なアルゴリズムを決めていて、それは MATLAB で `fft` をした場合も同じです。

>FFTW3 の doc から引用（http://www.fftw.org/fftw3_doc/Introduction.html#Introduction）
>> FFTW does not use a fixed algorithm for computing the transform, but instead it adapts the DFT algorithm to details of the underlying hardware in order to maximize performance. Hence, the computation of the transform is split into two phases. First, FFTW’s planner “learns” the fastest way to compute the transform on your machine.

決定方法には 4 つのメソッドが用意されており、ESTIMATE の場合は配列サイズから推定しますが、それ以外の場合は実際にフーリエ変換を実行した上で決定するなど、下に行くにつれて最適なアルゴリズムの検索に気合が入ってきます（＝時間がかかります）。

- FFTW_ESTIMATE
- FFTW_MEASURE
- FFTW_PATIENT
- FFTW_EXHAUSTIVE

これら Planner と呼ばれるものの MATLAB 側での入り口になるのが [`fftw`](https://jp.mathworks.com/help/matlab/ref/fftw.html?s_eid=PSM_29435) です。

```matlab
method = fftw('planner')
```

で確認できます。ちなみに MATLAB だと ESTIMATE がデフォルトになっていますが、変更する場合は

```matlab
fftw('planner','estimate');
fftw('planner','measure');
fftw('planner','patient');
```

などなど。


# 検証１：Plan の作成にかかる時間

以下の通りメソッド間で大きな違いが見て取れます。

MATLAB では plan の作成だけを行う処理はなく、`fft` の初回実行時に plan 作成されます。ですので、以下の時間には plan 作成時間に加えて `fft` の実行にかかった時間も入ります。

```
経過時間は 0.000670 秒です。 % estimate plan 作成 + fft
経過時間は 0.000051 秒です。 % 2回目（fft のみ）

経過時間は 0.021619 秒です。 % measure plan 作成 + fft
経過時間は 0.000062 秒です。 % 2回目（fft のみ）

経過時間は 0.117527 秒です。 % patient plan 作成 + fft
経過時間は 0.000064 秒です。 % 2回目（fft のみ）
```

Patient では Estimate の200倍くらい時間かけています。それぞれ2回目の fft 実行はかなり速くなっていますが、これは plan の作成処理が無いことに加えて、JIT が効いているからかなと。


以下のコードを実行しました。


```matlab
X = rand(503,1);
fftw('dwisdom',[]); % 初期化

fftw('planner','estimate');
tic; fft(X); toc;
tic; fft(X); toc;

fftw('planner','measure');
tic; fft(X); toc;
tic; fft(X); toc;

fftw('planner','patient');
tic; fft(X); toc;
tic; fft(X); toc;
```



# 検証２：Plan の fft 処理時間への影響

上で見た通り、初回の `fft` の実行は plan の作成処理も入ってしまいますので、初回以降の `fft` 実行の処理時間を比較してみます。一番重要なところですね。今回は複数回実行して中央値を返す `timeit` を使用します。

```
経過時間は 0.000023 秒です。 % estimate
経過時間は 0.000021 秒です。 % measure
経過時間は 0.000019 秒です。 % patient
```

多少（10%）早くなっている？という程度です。もっと大きな配列を使う場合で、総計算時間に FFT が占める割合が大きければ、patient/measure の使用も考慮してもいいかもしれません。

以下のコードを実行しました。

```matlab
fftw('dwisdom',[]); % 初期化

X = rand(503,1);
f = @() fft(X);

fftw('planner','estimate');
t = timeit(f);
disp("経過時間は "+ sprintf('%.6f',t) + " 秒です。")

fftw('planner','measure');
t = timeit(f);
disp("経過時間は "+ sprintf('%.6f',t) + " 秒です。")

fftw('planner','patient');
t = timeit(f);
disp("経過時間は "+ sprintf('%.6f',t) + " 秒です。")
```



# 検証３：事前に作成した plan の使用

さて plan の作成時間が `fft` の実行に比べてかなり時間がかかっていることが分かりました。そうなると、事前に plan だけ作成しておいて再利用することが有効な場合もあります。

同じ配列サイズであれば、plan の作成は初回実行のみ。なので、実際問題で影響は少ないはずですが、2回目以降の `fft` 処理速度を初回から発揮させたい、などの貪欲な要望に応えることができます。

まず plan 作成して mat ファイルに保存します。

```matlab

X = rand(503,1);
fftw('dwisdom',[]); % 初期化

fftw('planner','estimate');
fft(X); % plan 作成
fftinfo = fftw('dwisdom');
save('fftinfo_estimate.mat','fftinfo'); % 保存

fftw('planner','measure');
fft(X); % plan 作成
fftinfo = fftw('dwisdom');
save('fftinfo_measure.mat','fftinfo'); % 保存

fftw('planner','patient');
fft(X); % plan 作成
fftinfo = fftw('dwisdom');
save('fftinfo_patient.mat','fftinfo'); % 保存
```

それぞれの plan (`fftinfo`) が保存できました。`fftinfo` は文字列で定義されています。

```
(FFTW-3.3.3 fftw_wisdom #x3c273403 #x192df114 #x4d08727c #xe98e9b9d)
```

---

さて、次にそれぞれ個別に読み込んだうえで `fftw('dwisdom',fftinfo);` と設定し `fft` 実行にかかる時間を調べてみます。
 

```
経過時間は 0.000341 秒です。 % 初回(estimate)
経過時間は 0.000082 秒です。 % 2回目

経過時間は 0.000307 秒です。 % 初回(measure)
経過時間は 0.000088 秒です。 % 2回目

経過時間は 0.000413 秒です。 % 初回(patient)
経過時間は 0.000052 秒です。 % 2回目
```

plan の作成を行っていないので初回実行もかなり速くなっています（注：patient plan 作成は 0.12 秒かかっていました）し、各メソッドの違いによる影響は（ほぼ）見られません。2回目が速いのは JIT でしょう。たぶん。

以下のコードを実行しました。


```matlab

fftw('dwisdom',[]); % 初期化
X = rand(503,1);

load('fftinfo_estimate.mat','fftinfo');
fftw('dwisdom',fftinfo);
tic; fft(X); toc;
tic; fft(X); toc;

load('fftinfo_measure.mat','fftinfo');
fftw('dwisdom',fftinfo);
tic; fft(X); toc;
tic; fft(X); toc;

load('fftinfo_patient.mat','fftinfo');
fftw('dwisdom',fftinfo);
tic; fft(X); toc;
tic; fft(X); toc;
```

# まとめ

生の FFTW のライブラリを使ったことがある方はご存知、 planner のお話でした。MATLAB の `fft` は FFTW が基になっているので plan を活用しているんですが、基本的には気にしなくて済む実装になっていますね。

単発で配列サイズが大きい `fft` を実行する場合など、限られた状況において `fftw` 関数を使って plan をいじることで、計算速度の面でメリットが出る可能性はありますので、もし機会があったら試してみてください。

こんなケースだと大きな違いがでる、などありましたらコメントもお待ちしています。