2024年04月01日(月) [長年日記]
■ [c++][dev][math] 高速フーリエ変換
高速フーリエ変換についてメモ。
高速フーリエ変換とは、離散フーリエ変換を高速に計算するアルゴリズムのこと。
離散フーリエ変換
離散フーリエ変換とは、項数 N の複素数列 \((a_n)_{n=0}^{N-1}\) から、次で定める項数 N の複素数列 \((b_n)_{n=0}^{N-1}\) を求める操作のこととする。
\[ \ b_n = \sum_{x=0}^{N-1}{a_x \exp(-i \frac{2 \pi nx}{N})} \]
\(W_N=\exp\left(-i\frac{2 \pi}{N}\right)\) とおくと、この式は次のように書ける。
\[ \ b_n = \sum_{x=0}^{N-1}{W_N^{nx} a_x} \]
この式の掛け算の回数は N 回で、すべての 0≦n≦N に対して計算するには N×N 回の掛け算が必要になる。
高速フーリエ変換の解説は N を2のべき乗に制限するものが多いが、ここでは非素数の N に対して（素朴に計算するのに比べて）高速に計算できるアルゴリズムを考える。なお、以下の話は N が素数であっても適用できる（R=N、S=1 とすればよい）が、その場合は計算量が減らない。
計算量を減らすための式変形
R を N の約数、S を N/R とする。求めたい複素数列 \((b_n)_{n=0}^{N-1}\) を項数 S の R 個の複素数列 \((b_{Rm})_{m=0}^{S-1}, (b_{Rm+1})_{m=0}^{S-1}, ..., (b_{Rm+R-1})_{m=0}^{S-1} \) に分割する。数列 \((b_{Rm+k})_{m=0}^{S-1}\) の一般項は次のように表せる（この式変形の導出は本稿の最後に示す）。
\[ \ b_{Rm+k} = \sum_{y=0}^{S-1}{W_S^{my} \left( W_N^{ky} \sum_{x=0}^{R-1}{W_R^{kx} a_{Sx+y}} \right)} \]
この式の掛け算の回数を考えると、
- 括弧内の計算に R+1 回の掛け算が必要で、これをすべての 0≦y≦S-1 と 0≦k≦R-1 に対して計算するので全体では (R+1)×S×R 回
- \(W_S^{my}\)との積が S 回でこれをすべての 0≦m≦S-1 と 0≦k≦R-1 に対して計算するので全体では S×S×R 回
になる（と思う）。合わせて (R+1)×S×R + S×S×R = (R+S+1)×N 回の掛け算が必要で、これはほとんどのケースで元の式の N×N 回より少ない回数になる。括弧内の計算がmに依存しておらず計算結果を再利用できることが効いている。ただし、元の式には無かった \(W_S^{my}\) と \(W_R^{kx}\) の計算が新たに発生していることには注意。また、\((b_n)_{n=0}^{N-1}\) を分割する際に並び順が変わっているので、並び替えのコストも発生する。
ここで、\(b_{Rm+k}\) の式が項数 S の複素数列 \(\left( W_N^{ky} \sum_{x=0}^{R-1}{W_R^{kx} a_{Sx+y}} \right)_{y=0}^{S-1}\) に対する離散フーリエ変換の形になっていることに注目してほしい。これにより、R個の数列 \((b_{Rm+k})_{m=0}^{S-1}\) それぞれの計算に対しても上記と同じ式変形で計算量を減らすことができ、この操作は数列の要素数が1になるまで再帰的に行えることになる。
高速フーリエ変換の実装
上のようにして計算量を減らしたコードと減らさず素朴に計算するコードの比較を行うプログラムを書いた。200行くらいあるので全体はGistに置きここにコードを抜粋する。コンパイルにはC++20が必要だが、もっと古いC++でコンパイルできるように直すのは難しくない。
高速フーリエ変換の計算をする関数 fft() は次の通り。
// inputの値を離散フーリエ変換した結果を返す
template<typename T>
std::vector<std::complex<T>> fft(
std::vector<std::complex<T>> input,
const twiddle_factor_map<T>& tf_map)
{
std::vector<std::complex<double>> work(input.size());
detail::fft(
input.begin(),
input.end(),
work.begin(),
detail::find_divisor(input.size()),
tf_map);
return input;
}
- 入力をvectorで与え結果をvectorで受け取る。
- \( W_N \) のべき乗については事前に計算したものを tf_map というマップに設定して渡すようにした。
- 計算の作業領域用に入力と同じサイズのvectorを使用している。
上で呼び出している実際の処理を行う関数 detail::fft() は次の通り。
// 範囲[first, last)の値を離散フーリエ変換した結果を範囲[first, last)に上書きする。
// work_firstは[first, last)と同じサイズの作業用範囲の先頭。
// radixには範囲[first, last)の要素数の約数を指定する。
template<typename T, typename Iter>
void fft(
Iter first,
Iter last,
Iter work_first,
std::size_t radix,
const twiddle_factor_map<T>& tf_map)
{
const auto n = static_cast<std::size_t>(std::distance(first, last));
if (n <= 1) return;
const auto& tfs_r = tf_map.at(radix);
const auto& tfs_n = tf_map.at(n);
const auto s = n / radix;
for (std::size_t i = 0; i < s; ++i) {
auto it = work_first + i;
for (std::size_t j = 0; j < radix; ++j) {
*it = 0;
for (std::size_t k = 0; k < radix; ++k) {
*it += *(first + i + s * k) * tfs_r[j * k % radix];
}
*it *= tfs_n[i * j];
it += s;
}
}
for (std::size_t i = 0; i < radix; ++i) {
const auto it = work_first + s * i;
fft(it, it + s, first, find_divisor(s), tf_map);
}
for (std::size_t i = 0; i < s; ++i) {
auto it = first + radix * i;
for (std::size_t j = 0; j < radix; ++j) {
*it = *(work_first + i + s * j);
++it;
}
}
}
- detail::fft() は自分自身を呼び出している。再帰呼び出しの中でvectorを作り直さなくていいように、イテレータ範囲でデータを渡すようにしている。
- 作業領域内にR個の複素数列 \(\left( W_N^{ky} \sum_{x=0}^{R-1}{W_R^{kx} a_{Sx+y}} \right)_{y=0}^{S-1}\) を設定している。
- 再帰呼び出しでは作業領域を入力データ、入力データ範囲を作業領域として渡している。
- 再帰呼び出しが終わったら作業領域内のデータを並び替えたものを計算結果として入力データ範囲に設定している。
実行結果
Gistに置いたファイルの実行結果は以下のような感じ。素朴に計算した場合の実行時間（dft）を高速に実行した場合の実行時間（fft）で割った値を ratio として出力している。delta は計算結果の差の絶対値の最大値。
なお、NumPyのFFTルーチンと比べると遅い（項数 129140163 = 3の17乗 のときに上のコードは23秒、NumPyは10秒で計算）ので、上のコードもしくは本アルゴリズムは改善の余地が大いにある。
入力の項数が少ないうちは素朴に計算する方が速い。
size 2 ratio 0.379 dft 0.000 fft 0.000 delta 0.000e+00 size 3 ratio 0.858 dft 0.000 fft 0.000 delta 0.000e+00 size 4 ratio 0.897 dft 0.000 fft 0.000 delta 0.000e+00 size 5 ratio 0.950 dft 0.000 fft 0.000 delta 0.000e+00 size 6 ratio 0.400 dft 0.000 fft 0.000 delta 7.589e-15 size 7 ratio 0.907 dft 0.000 fft 0.000 delta 0.000e+00 size 8 ratio 0.640 dft 0.000 fft 0.000 delta 3.794e-15 size 9 ratio 2.072 dft 0.000 fft 0.000 delta 0.000e+00 size 10 ratio 0.621 dft 0.000 fft 0.000 delta 7.944e-15 size 11 ratio 0.966 dft 0.000 fft 0.000 delta 0.000e+00 size 12 ratio 0.600 dft 0.000 fft 0.000 delta 2.197e-14 size 13 ratio 1.021 dft 0.000 fft 0.000 delta 0.000e+00 size 14 ratio 0.704 dft 0.000 fft 0.000 delta 1.507e-14 size 15 ratio 0.524 dft 0.000 fft 0.000 delta 1.888e-14
項数2000くらいで、大きくて40倍くらいの差が出てくる。項数が素数の場合は速くならない。
size 1987 ratio 0.998 dft 0.010 fft 0.010 delta 0.000e+00 size 1988 ratio 22.035 dft 0.010 fft 0.000 delta 3.122e-12 size 1989 ratio 46.875 dft 0.010 fft 0.000 delta 2.890e-12 size 1990 ratio 8.610 dft 0.010 fft 0.001 delta 2.621e-12 size 1991 ratio 10.168 dft 0.009 fft 0.001 delta 2.390e-12 size 1992 ratio 19.410 dft 0.010 fft 0.000 delta 2.912e-12 size 1993 ratio 1.000 dft 0.009 fft 0.009 delta 0.000e+00 size 1994 ratio 1.997 dft 0.009 fft 0.005 delta 3.429e-12 size 1995 ratio 46.528 dft 0.010 fft 0.000 delta 2.430e-12 size 1996 ratio 3.932 dft 0.010 fft 0.002 delta 4.590e-12 size 1997 ratio 1.039 dft 0.010 fft 0.009 delta 0.000e+00 size 1998 ratio 33.008 dft 0.010 fft 0.000 delta 2.993e-12 size 1999 ratio 0.892 dft 0.010 fft 0.011 delta 0.000e+00
項数が2のべき乗の場合。数が大きくなるとどんどん差がつく。
--- 2のべき乗 --- size 2 ratio 1.911 dft 0.000 fft 0.000 delta 0.000e+00 size 4 ratio 0.982 dft 0.000 fft 0.000 delta 0.000e+00 size 8 ratio 0.741 dft 0.000 fft 0.000 delta 7.105e-15 size 16 ratio 0.908 dft 0.000 fft 0.000 delta 1.137e-14 size 32 ratio 1.188 dft 0.000 fft 0.000 delta 6.071e-14 size 64 ratio 2.473 dft 0.000 fft 0.000 delta 2.050e-13 size 128 ratio 4.226 dft 0.000 fft 0.000 delta 3.422e-13 size 256 ratio 8.059 dft 0.000 fft 0.000 delta 1.833e-12 size 512 ratio 15.466 dft 0.001 fft 0.000 delta 2.854e-11 size 1024 ratio 28.213 dft 0.003 fft 0.000 delta 8.084e-11 size 2048 ratio 56.931 dft 0.012 fft 0.000 delta 3.748e-10 size 4096 ratio 104.531 dft 0.044 fft 0.000 delta 7.867e-10 size 8192 ratio 204.091 dft 0.172 fft 0.001 delta 1.531e-09 size 16384 ratio 361.594 dft 0.644 fft 0.002 delta 2.899e-08 size 32768 ratio 715.049 dft 2.648 fft 0.004 delta 6.994e-08 size 65536 ratio 1349.311 dft 11.721 fft 0.009 delta 3.444e-07
項数が3のべき乗の場合。
--- 3のべき乗 --- size 3 ratio 1.306 dft 0.000 fft 0.000 delta 0.000e+00 size 9 ratio 1.336 dft 0.000 fft 0.000 delta 0.000e+00 size 27 ratio 1.078 dft 0.000 fft 0.000 delta 4.550e-14 size 81 ratio 2.865 dft 0.000 fft 0.000 delta 2.558e-13 size 243 ratio 8.203 dft 0.000 fft 0.000 delta 2.752e-12 size 729 ratio 23.335 dft 0.001 fft 0.000 delta 1.347e-11 size 2187 ratio 57.941 dft 0.013 fft 0.000 delta 1.113e-10 size 6561 ratio 169.265 dft 0.116 fft 0.001 delta 2.806e-09 size 19683 ratio 466.513 dft 1.041 fft 0.002 delta 4.331e-08 size 59049 ratio 1307.727 dft 9.614 fft 0.007 delta 8.864e-08
項数の素因数が異なる数になる場合。
--- 素数の積 --- size 2 ratio 1.211 dft 0.000 fft 0.000 delta 0.000e+00 size 6 ratio 0.452 dft 0.000 fft 0.000 delta 7.119e-15 size 30 ratio 1.100 dft 0.000 fft 0.000 delta 3.197e-14 size 210 ratio 6.740 dft 0.000 fft 0.000 delta 2.975e-12 size 2310 ratio 61.381 dft 0.014 fft 0.000 delta 6.370e-11 size 30030 ratio 576.535 dft 2.444 fft 0.004 delta 1.515e-07
項数が素数の場合は同じくらいの速さ。
--- 素数 --- size 2 ratio 0.943 dft 0.000 fft 0.000 delta 0.000e+00 size 3 ratio 0.777 dft 0.000 fft 0.000 delta 0.000e+00 size 5 ratio 0.982 dft 0.000 fft 0.000 delta 0.000e+00 size 7 ratio 1.054 dft 0.000 fft 0.000 delta 0.000e+00 size 11 ratio 1.069 dft 0.000 fft 0.000 delta 0.000e+00 size 13 ratio 1.436 dft 0.000 fft 0.000 delta 0.000e+00 size 1009 ratio 0.999 dft 0.003 fft 0.003 delta 0.000e+00 size 10009 ratio 0.967 dft 0.269 fft 0.278 delta 0.000e+00 size 100003 ratio 0.988 dft 28.695 fft 29.054 delta 0.000e+00
参考
次のページが参考になった。
- FFT (高速フーリエ・コサイン・サイン変換) の概略と設計法 （京都大学数理解析研究所）
上で行った変換は混合基数アルゴリズムに該当する。
式変形の詳細
最後に、先に示した \((b_{Rm+k})_{m=0}^{S-1}\) の一般項の導出をメモしておく。
計算の過程で次を使う。xは整数、Nは自然数、RはNの約数、Sは N/R とする。
\[ \ W_N^{Nx} = 1, \quad \ W_N^{Rx} = W_S^x, \quad \ \sum_{x=0}^{N-1}{f(x)} = \sum_{y=0}^{S-1}{\sum_{x=0}^{R-1}{f(Sx+y)}} \]
これらを使うことで次のように変形できる。
\begin{align} \ b_{Rm+k} &= \sum_{x=0}^{N-1}{W_N^{(Rm+k)x} a_x} \\ \ &= \sum_{y=0}^{S-1}{\sum_{x=0}^{R-1}{W_N^{(Rm+k)(Sx+y)} a_{Sx+y}}} \\ \ &= \sum_{y=0}^{S-1}{\sum_{x=0}^{R-1}{W_N^{Nmx} W_N^{Rmy} W_N^{kSx} W_N^{ky} a_{Sx+y}}} \\ \ &= \sum_{y=0}^{S-1}{\sum_{x=0}^{R-1}{W_S^{my} W_R^{kx} W_N^{ky} a_{Sx+y}}} \\ \ &= \sum_{y=0}^{S-1}{W_S^{my} W_N^{ky} \sum_{x=0}^{R-1}{W_R^{kx} a_{Sx+y}}} \\ \end{align}