高速フーリエ変換（FFT）の転置写像を観察することで時間間引き・周波数間引きのアルゴリズムが相互に導出できることを紹介する。時間間引き・周波数間引きの具体的なアルゴリズムについては知っていることを前提にする（sugarknriさんのFFTの解説）。
クロネッカー積
\(m \times n\) 行列 \(A\) と任意のサイズの行列 \(B\) のクロネッカー積を
$$
A\otimes B := \begin{bmatrix}
a_{11} B & \cdots & a_{1n} B \\
\vdots & \ddots & \vdots \\
a_{m1} B & \cdots & a_{mn} B
\end{bmatrix}$$で定義する。
解説
数列 \((x_i)_{1\leq i\leq n}\) の離散フーリエ変換（DFT）\((\bar{X}_i)_{1 \leq i\leq n}\) を原始 \(n\) 乗根 \(\zeta\) を用いて$$\bar{X_k}:=\sum_{m=1}^n x_m\zeta^{km}$$と定義する。この変換は行列 \(W=(\zeta^{ij})_{1\leq i,j\leq n}\) を用いて
$$
\begin{pmatrix}
\zeta^0 & \zeta^0 & \zeta^0 & \cdots & \zeta^0 \\
\zeta^0 & \zeta^1 & \zeta^2 & \cdots & \zeta^{n-1} \\
\zeta^0 & \zeta^2 & \zeta^4 & \cdots & \zeta^{2(n-1)} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
\zeta^0 & \zeta^{n-1} & \zeta^{2(n-1)} & \cdots & \zeta^{(n-1)(n-1)}
\end{pmatrix} \begin{pmatrix}
x_0 \\ x_1 \\ x_2 \\ \vdots \\ x_{n-1}
\end{pmatrix} = \begin{pmatrix}
\bar X_0 \\\bar X_1 \\ \bar X_2 \\ \vdots \\ \bar X_{n-1}
\end{pmatrix}
$$
と表せる。
\(I_i\) を \(i\) 次の単位行列、\(Z_i:=\mathrm{diag} (\zeta^{\frac{n}{i}0},\zeta^{\frac{n}{i}1},\ldots,\zeta^{\frac{n}{i}{(i-1)}})\)、\(P\) を \(n\) 次の bit-reversal permutation を表す行列、
$$T_i:=\begin{bmatrix}
I_{2^{i-1}} & Z_{2^{i-1}} \\ I_{2^{i-1}} & -Z_{2^{i-1}}
\end{bmatrix} \otimes I_{n/2^{i-1}}
$$とすると時間間引きの FFT は \(W=T_k \ldots T_2T_1P\) という分解で表される。\(W\) は対称行列であるから転置を取ると \(W=W^\top =P^\top T_1^\top T_2^\top \ldots T_k^\top\) 、
$$T^\top=\begin{bmatrix}
I_{2^{i-1}} & I_{2^{i-1}} \\ Z_{2^{i-1}} & -Z_{2^{i-1}}
\end{bmatrix} \otimes I_{n/2^{i-1}}
$$
である。\(P\) は disjoint な互換の積だから \(P=P^\top=P^{-1}\) 。これはよく見ると周波数間引きの FFT そのものである。つまり、転置を取ることで時間間引きのアルゴリズムから周波数間引きのアルゴリズムが得られた。式の形から、FFT に周波数間引き、inverse-FFT に時間間引きを使えば、\(P\) は打ち消し合って計算の必要がないことも分かる。
実装
転置を取ることで時間間引きの実装から周波数間引きの実装を機械的に書くことができる。下の2つのコードはそれぞれ時間間引き・周波数間引きの C++ のコードである。まず、転置を取ると \((AB)^\top=B^\top A^\top\) のようになるから、配列 a に関数する操作の順序を全て逆にする。これは for 文を逆向きに回すだけで良い。3重 for 文 の内部は
$$\begin{bmatrix}
a[\mathrm{pos}] \\ a[\mathrm{pos+len}] \end{bmatrix} \leftarrow \begin{bmatrix}
1 & x \\ 1 & -x
\end{bmatrix} \begin{bmatrix}
a[\mathrm{pos}] \\ a[ \mathrm{pos+len} ] \end{bmatrix} $$
と表せるから、作用している行列を転置して
$$\begin{bmatrix}
a[\mathrm{pos}] \\ a[\mathrm{pos+len}] \end{bmatrix} \leftarrow \begin{bmatrix}
1 & 1 \\ x & -x
\end{bmatrix} \begin{bmatrix}
a[\mathrm{pos}] \\ a[ \mathrm{pos+len} ] \end{bmatrix} $$
とすればよい。実際に下の周波数間引きのコードは、時間間引きのコードから機械的に生成しており、その動作原理は全く考えずに書いた（が、正しい結果を返す）。
時間間引きのコード
constexpr long long p = 998244353;
constexpr long long g = 3;
void fft(std::vector &a, long long g) {
int n = a.size();
g = pow_mod(g, (p - 1) / n);
{
int i = 0;
for (int j = 0; j < n; ++j) {
if (i > j) std::swap(a[i], a[j]);
for (int k = n / 2; k > (i ^= k); k /= 2);
}
}
for (int len = 1; len <= n / 2; len *= 2) {
long mul = pow_mod(g, n / (2 * len));
for (int src = 0; src < n; src += 2 * len) {
long long x = 1;
for (int pos = src; pos < src + len; ++pos) {
long long s = a[pos];
long long t = a[pos + len] * x % p;
a[pos] = (s + t) % p;
a[pos + len] = (s + p - t) % p;
x = x * mul % p;
}
}
}
} 周波数間引きのコード
constexpr long long p = 998244353;
constexpr long long g = 3;
void fft(std::vector &a, long long g) {
int n = a.size();
g = pow_mod(g, (p - 1) / n);
for (int len = n / 2; len >=1; len /= 2) {
long mul = inv(pow_mod(g, n / (2 * len)));
long src_x = pow_mod(mul, len - 1);
for (int src = n - 2 * len; src >=0; src -= 2 * len) {
long x = src_x;
for (int pos = src + len - 1; pos >= src; --pos) {
long long u = a[pos];
long long v = a[pos + len];
a[pos] = (u + v) % p;
a[pos + len] = (u + p - v) % p * x % p;
x = x * mul % p;
}
}
}
{
int i = 0;
for (int j = 0; j < n; ++j) {
if (i > j) std::swap(a[i], a[j]);
for (int k = n / 2; k > (i ^= k); k /= 2);
}
}
} 余談
高速フーリエ変換では転置写像が元の写像と一致したたため、転置写像が元の問題の非自明なアルゴリズムを与えた。しかし、転置写像が元の写像に一致せずとも、転置写像を観察することで元の問題の非自明なアルゴリズムを導けることがある。例えば、A. Bostan らは転置写像を経由することで multipoint evaluation や interpolation の計算量を定数倍改善した。転置写像と元の写像に本質的な違いはないはずなのに、転置を取ることで高速な解法の思いつきやすくなるという事に驚いた。転置を取るだけで元の問題を非自明な角度から見ることができるというのは意外で面白いと思う。
参考文献
- Bostan, Alin, Grégoire Lecerf, and Éric Schost. “Tellegen’s principle into practice.” Proceedings of the 2003 international symposium on Symbolic and algebraic computation. 2003.
- 転置写像を考察することで multipoint evalution や interpolation の定数倍高速化をしています。multipoint evalution については競プロLTで tko919 さんが解説してくれました。
- Charles Van Loan 氏の授業スライド
Comments