円周率.jp > 多数桁計算 > FFT > 2次元化FFT

2次元化FFT

多数桁乗算に用いられる FFT で計算対象になる数値は, 一般に桁数 n の 1 次元配列でもって表される. FFT の演算ではほぼランダムアクセスに近いメモリアクセスをするため, n がキャッシュメモリより大きくなると極端に演算速度が遅くなる.

これを回避するため,npq と因数分解できると仮定し,FFT の計算式を次のように書き換える.

kk1pk2jj1qj2     (0≦k1j2q, 0≦k2j1p
Ak1p+k2 \sum_{j_1=0}%5e{p-1} \sum_{j_2=0}%5e{q-1} aj1q+j2 ω(j1q+j2)(k1p+k2)
\sum_{j_1=0}%5e{p-1} \sum_{j_2=0}%5e{q-1} aj1q+j2 ωj1k2q+j2k1p+j2k2
\sum_{j_2=0}%5e{q-1} \Bigl( \Bigl( \sum_{j_1=0}%5e{p-1} aj1q+j2 ωpj1k2 \Bigr) ω j2k2 \Bigr) ωqj2k1

このように変形することで,pq 程度の長さの FFT を繰り返すことになるので,pq 次第ではキャッシュメモリに収まるので演算全体の高速化が期待できる. 基本的には max (pq) を最小化するため, pq\sqrt{n} にする.

2 次元化 (four-step FFT)

上記の式では若干わかりにくいので,入力となる 1 次元データ列を 2 次元に配置し再解釈することで分かりやすくなる. 大きく 4 つの計算ステップを踏むことから four-step FFT という.

a0[j1j2] aj1qj2    (配列の 2 次元再解釈)
1st STEP   a1[k2j2] \sum_{j_1=0}%5e{p-1} a0[j1j2] ωpj1k2    (for 0≦j2q
2nd STEP   a2[k2j2] a1[k2j2] ω j2k2    (for 0≦k2p,0≦j2q
3rd STEP   a3[j2k2] a2[k2j2]    (for 0≦k2p,0≦j2q
4th STEP   a4[k1k2] \sum_{j_2=0}%5e{q-1} a3[j2k2] ωqj2k1    (for 0≦k2p
Ak1p+k2 a4[k1k2]    (配列の 1 次元再解釈)

six-step FFT

four-step FFT では各 FFT を行う際に用いるデータが不連続に並んでいるため, 計算効率が非常に悪い. それを改善するために,最初と最後に行列を転置することで FFT でアクセスするメモリが連続になるようにしたものを six-step FFT という. 実際には STEP 3 のひねり係数をかける作業と STEP 4 の行列転置を同時に行うので, 実質的には five-step FFT と考えることもできる.

  1. 行列を転置する
  2. 1 次元 FFT を行う (four-step の 1st STEP 相当)
  3. ひねり係数 ωjk をかける (four-step の 2nd STEP 相当)
  4. 行列を転置する (four-step の 3rd STEP 相当)
  5. 1 次元 FFT を行う (four-step の 4th STEP 相当)
  6. 行列を転置する

3 次元化 FFT (nine-step FFT)

同様に 3 次元化も可能である.2 次元化で使う 1 次元 FFT もキャッシュサイズに載らない場合に使うことで効果が上がる. また,分散メモリ型並列計算にも応用させやすい. 式変形は面倒だが,1 度やっておくと理解が深まる. 変換前後のデータ ajAk の 3 軸をどう取るかについては, 変換後の式を考えて ω の項が少なくなるように取ると良い. 1 つの分割,計算例を以下に書く.

npqrkk1pqk2pk3jj1qrj2rj3, ( 0≦k3j1p, 0≦k2j2q, 0≦k1j3r
Ak1pqk2pk3 \sum_{j_1=0}%5e{p-1} \sum_{j_2=0}%5e{q-1} \sum_{j_3=0}%5e{r-1} aj1qrj2rj3 ω(j1qrj2rj3)(k1pqk2pk3)
\sum_{j_1=0}%5e{p-1} \sum_{j_2=0}%5e{q-1} \sum_{j_3=0}%5e{r-1} aj1qrj2rj3 ωj1k3qrj2k2prj3k1pqj2k3rj3k2pj3k3
\sum_{j_3=0}%5e{r-1} \Bigl( \Bigl( \sum_{j_2=0}%5e{q-1} \Bigl( \Bigl( \sum_{j_1=0}%5e{p-1} aj1qrj2rj3 ωpj1k3 \Bigr) ωpqj2k3 ωj3k3 \Bigr) ωqj2k2 \Bigr) ωqrj3k2 \Bigr) ωrj3k1

行列転置を含めて再解釈

3 次元化 FFT では最初と最後を含めて 4 回の行列転置が必要になる. また,分散メモリ型計算機を用いてこの計算を行う場合はそのうち 2 回((1)と(9))が全対全のデータ通信を用いての転置となるため, 通信速度がネックになりやすい.

a0[j1j2j3] aj1qrj2rj3    (配列の 3 次元再解釈)
(1) a1[j3j2j1] a0[j1j2j3]    (for 0≦j1p,0≦j2q,0≦j3r
(2) a2[j3j2k3] \sum_{j_1=0}%5e{p-1} a1[j3j2j1] ωpj1k3    (for 0≦j2q, 0≦j3r
(3) a3[j3j2k3] a2[j3j2k3] ωpqj2k3 ω j3k3    (for 0≦k3p,0≦j2q,0≦j3r
(4) a4[j3k3j2] a3[j3j2k3]    (for 0≦k3p,0≦j2q,0≦j3r
(5) a5[j3k3k2] \sum_{j_2=0}%5e{q-1} a4[j3k3j2] ωqj2k2    (for 0≦k3p,0≦j3r
(6) a6[j3k3k2] a5[j3k3k2] ωqrj3k2    (for 0≦k3p,0≦k2q,0≦j3r
(7) a7[k3k2j3] a6[j3k3k2]    (for 0≦k3p,0≦k2q,0≦j3r
(8) a8[k3k2k1] \sum_{j_3=0}%5e{r-1} a7[k3k2j3] ωrj3k1    (for 0≦k3p,0≦k2q
(9) a9[k1k2k3] a8[k3k2k1]    (for 0≦k3p,0≦k2q,0≦k1r
Ak1pq+k2p+k3 a9[k1k2k3]    (配列の 1 次元再解釈)