やったこと
fft 関数1を実行すると、まず FFT に使用するアルゴリズムを決定する処理が実行されます。
method = fftw('planner')
で表示されるアレです。場合によっては結構時間を食うこのステップをスキップする方法について紹介します。
計算処理への影響は、ほとんどの場合極小かと思いますが、単発で大きな配列に fft を実行する場合は、知っておいて損はしないかもしれません。
環境
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 です。
method = fftw('planner')
で確認できます。ちなみに MATLAB だと ESTIMATE がデフォルトになっていますが、変更する場合は
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 が効いているからかなと。
以下のコードを実行しました。
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 の使用も考慮してもいいかもしれません。
以下のコードを実行しました。
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 ファイルに保存します。
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 でしょう。たぶん。
以下のコードを実行しました。
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 をいじることで、計算速度の面でメリットが出る可能性はありますので、もし機会があったら試してみてください。
こんなケースだと大きな違いがでる、などありましたらコメントもお待ちしています。