- 高次の多項式を学校で習う方法で展開すると、すべての項の組を掛け合わせる必要があるため、O(n²) のコストがすぐにボトルネックになる
- 多項式の係数ベクトルの乗算は離散信号の 畳み込み と同じであり、
[2, 3, 4] と [5, 6, 7] の結果は [10, 27, 52, 45, 28] になる
- DFT は離散信号を周波数領域へ移し、FFT は同じ変換を O(n log n) で計算するため、大きな入力で差が出る
- 時間領域での畳み込みは周波数領域での要素ごとの乗算に変わるため、FFT で変換して掛け合わせた後に IFFT で戻せば、多項式の乗算をより高速に処理できる
- 低次では FFT/IFFT の往復コストが利点を相殺することもあるが、次数が大きくなるほど FFT 方式のほうが効率的である
多項式の乗算が遅くなる理由
- 多項式
P(x) は係数 a_k と変数 x のべき乗項を足し合わせた形で表される
- 例
P(x)=5x²+2x+9 は次数 2 の多項式である
- 係数ベクトルは表記方法によって
[5, 2, 9] または [9, 2, 5] のように表せる
- 加算と減算は同じ次数の項どうしを足したり引いたりすればよいため、比較的単純である
- Python では
zip(p, q) で各係数を走査しながら a + b または a - b を計算できる
- 次数が異なる場合は
zip_longest を使える
- 乗算は各項を互いに掛けた後、同じ次数の項を再び合計しなければならないため、計算量が大きくなる
(2x²+3x+4) × (5x²+6x+7) の結果は 10x⁴+27x³+52x²+45x+28 である- この方法の計算量は O(n²) であり、次数が大きくなるほど必要な乗算回数が増える
係数ベクトルと畳み込み
- 離散領域で 2 つの信号
p と q の畳み込みは y[n]=Σ p[k]·q[n-k] と定義される
- 計算は
q を反転したうえで p の上を左から右へ移動させ、重なった要素の積を足し合わせる方式である
- 例の信号は次のとおり
p = [2, 3, 4]q = [5, 6, 7]
q を反転して移動すると、各出力係数は次の順で作られる
2×5 = 102×6 + 3×5 = 272×7 + 3×6 + 4×5 = 523×7 + 4×6 = 454×7 = 28
- 畳み込み結果は
y = [10, 27, 52, 45, 28] である
- これは多項式の乗算で得られる
10x⁴+27x³+52x²+45x+28 の係数と同じである
- したがって 多項式の乗算 は係数ベクトルの畳み込みとして見ることができる
フーリエ変換と FFT
- フーリエ変換は信号を 時間領域 から 周波数領域 へ変換する
- 時間の観点では、信号を特定時点の値として見る
- 周波数の観点では、信号を異なる振動周波数の和として解釈する
- 振動周波数はサイン波とコサイン波で表され、それぞれ係数と位相を持つ
- 5Hz の純粋なサイン波に FFT を適用すると、周波数領域では 5Hz の位置にデルタのように現れる
- これは時間領域のサイン波が 5Hz のサイン 1 つで表現できることを示している
- 関連用語は次のように区別される
- Fourier Transform(FT): 連続領域で定義されたフーリエ変換
- Discrete Fourier Transform(DFT): 離散信号に対して定義されたフーリエ変換
- Fast Fourier Transform(FFT): DFT を O(n²) の代わりに O(n log n) で計算するアルゴリズム
- DFT は離散時間信号
x[n] を周波数領域の X[k] に変換する
- 各
X[k] は、入力サンプルを特定の周波数を表す複素数と掛け合わせて足し合わせることで計算される
周波数領域で乗算に変える
- DFT と周波数領域の最大の利点は、畳み込みを要素ごとの乗算に変えられる 点にある
- 時間領域で 2 つの信号を畳み込むことは、周波数領域で 2 つの信号を掛け合わせることと同じである
- 乗算は畳み込みよりも高速に計算できる
- 多項式の乗算を高速に行う手順は次のとおり
- 多項式を FFT で周波数領域へ変換する: O(n log n)
- 周波数領域で要素ごとに掛ける: O(n)
- 結果を IFFT で再び時間領域へ戻す: O(n log n)
- 全体として、FFT を使えば多項式の乗算を O(n log n) の計算量で実行できる
- 大きな多項式では、学校で習う O(n²) の乗算より高速である
Python 実装とベンチマーク
multiply_naive は二重ループですべての係数の組を掛け合わせ、結果位置 i + j に加算する
- 結果の長さは
len(p) + len(q) - 1 である
- 計算量は O(n²) である
multiply_fft は FFT/IFFT ベースで係数の乗算を行う
- 結果の長さを収められるように、
len(p) + len(q) - 1 以上の 2 のべき乗長を計算する
np.pad で 2 つの入力をパディングするnp.fft.fft で変換した値を要素ごとに掛けるnp.fft.ifft で戻した後、実部を丸めて整数係数へ変換する
- 例の入力
p = [2, 3, 4], q = [5, 6, 7] では、どちらの方法も [10, 27, 52, 45, 28] を返す
- ベンチマークでは
multiply_naive の代わりに、np.convolve を使う multiply_convolve と FFT 方式を比較する
multiply_naive は Python のループが遅く、np.fft.fft を使う FFT 方式と直接比較しにくいためであるnp.convolve は同じ演算を低水準の C コードで実行する
- 次数は
range(1, 30000, 1000) の範囲で増やし、各次数で 1 から 999999 の間のランダム係数を持つ 2 つの多項式を生成する
- 各方式は
n_runs = 5 で平均時間を測定する
- 低次では FFT/IFFT の往復変換コストのため、FFT 方式が有利でないこともある
- 次数が大きくなると、FFT 方式がはるかに効率的な結果を示す
1件のコメント
Hacker News のコメント
こういう説明でいつも気になるのは、たいてい数値誤差を忘れている点。
係数の乗算を単に「定数時間」として抽象化することはできない。そうするなら、そもそも乗算全体を抽象化しても同じことになる。数値精度を考慮すると、O(n (log n)^3) により近い [1]
[1]: http://numbers.computation.free.fr/Constants/Algorithms/fft....
[1] One-Dimensional Quaternion Discrete Fourier Transform and an Approach to Its Fast Computation:
https://www.mdpi.com/2079-9292/12/24/4974
[2] Convolution Theorems for Quaternion Fourier Transform: Properties and Applications:
https://onlinelibrary.wiley.com/doi/10.1155/2013/162769
[3] On the Matrix Form of the Quaternion Fourier Transform and Quaternion Convolution:
https://arxiv.org/abs/2307.01836
この方法で長い数同士を掛け算できる。要点は、多項式の乗算が、桁上がり（carry）をしない通常の長い数の乗算と同じだということ。
たとえば 1000 桁の数があるなら、各桁を 1000 個の要素を持つ多項式の係数とみなす。次に、記事で説明されている FFT の方法でこれらの多項式を掛け合わせる。結果を数に戻すには、桁上がりを処理する必要がある。ある要素が 10 より大きければ、超過分を次の桁へ渡し、係数を数に変換すればよい。
基本的な考え方はこうで、桁上がりに必要な精度や、FFT の結果を最も近い整数に丸めても正しいことを保証する部分には微妙な点がある。この方式は、この分野の代表的なライブラリである GMP が大きな数の乗算を行う方法である
実際に意味が出るには数がどれくらい大きい必要があるのか、そしてどこで使われるのかが気になる
まだ見ていないなら、この動画を見るとよい
https://youtu.be/h7apO7q16V0?si=bmgUEMTQSqU3flIv
多項式の乗算から FFT アルゴリズムを導出していて、本当に素晴らしい。半年に一度くらい見返している
FFT の「畳み込みは点ごとの乗算」という性質は、任意の巡回乗法群でも成り立つ。より代数的な導出については https://www.sciencedirect.com/science/article/pii/S002200007... を参照。
これを「調和 FFT」と呼ぶ場合もあり、非調和 FFT もある：GF(2^n) 上の [LCH14] “additive NTT”、有限体の単位円 X^2+Y^2=1 上の [HLP24] circle FFT、楕円曲線の同種写像列上の [BCKL21] ecfft
[LCH14]: https://arxiv.org/abs/1404.3458
[HLP24]: https://eprint.iacr.org/2024/278
[BCKL21]: https://arxiv.org/pdf/2107.08473
より高速な多項式乗算に FFT を使うことを最初に提案したのは誰だったのだろう？
最近気になって調べてみたが、引用をうまくたどれたわけではないものの、David Eppstein の 1995 年の論文 [0] まではさかのぼれた。ここでは、漸進的な更新後の部分和問題を効率よく解くためにこれを使っている。きっと Knuth の TAOCP にはもっと早い時期に載っていそうだ
FFT 多項式乗算によって、重複を許す厳密な部分和問題を準指数時間で解けるという事実もかなり衝撃的だった [1]。重要なのは、このアルゴリズムは O(N log N) だが、ここでの N は集合のサイズではなく最大要素なので、P ≠ NP に対する反例のようなものではないという点だ
[0] https://escholarship.org/content/qt6sd695gn/qt6sd695gn.pdf
[1] https://x.com/festivitymn/status/1788362552998580473?s=46&t=...
Strassen は 1968 年に Pollard 風のアプローチを発見したとされているが、文書としての記録はない。また、FFT そのものの誕生ではないにしても、Cooley-Tukey の 1965 年の論文 [4] が FFT とその応用研究を本格的に促した点も見るべきだ。その数年後の出来事である
[1] https://doi.org/10.1090/S0025-5718-1971-0301966-0
[2] https://doi.org/10.1016/S0022-0000(71)80014-4
[3] https://doi.org/10.1007/BF02242355
[4] https://doi.org/10.1090/S0025-5718-1965-0178586-1
https://www.cis.rit.edu/class/simg716/FFT_Fun_Profit.pdf
すべての機械学習は畳み込み方程式を解くことだと思っている
この論文は強化学習の文脈で扱っているが https://arxiv.org/abs/1712.06115、ほとんどのアプローチはこのパラダイムに当てはまる
たった今、FFT を使って時系列の部分列の大きな集合に対する内積を計算するアルゴリズム（matrix profile）を実装した。時系列の長さ n は数億規模まであり得る
FFT を使った高速な畳み込み計算により、計算時間は O(n) から O(log n) に減り、この規模では速度向上がすさまじい。GPU まで使えば、ノート PC で 1,000 万個のデータポイントを 0.1 秒で処理する、といった具合にさらに速くなる
この演算の核心となる「トリック」は、この気づきのように見える:
「時間領域で 2 つの信号の畳み込みを行うことは、周波数領域で 2 つの信号を掛け合わせることと同じ」という性質があり、FFT は時間領域から周波数領域へ変換できるようにしてくれる。だから多項式を FFT で周波数領域へ移した後、その領域では乗算するだけでよい。畳み込みより速い。この説明で抜けていたステップがはっきりするか気になるし、抜けている部分があれば記事を更新できる
だとすると、整数因数分解は離散デコンボリューションなのだろうか？ FFT 表現、つまり点ごとの乗算の逆演算と、tableax、つまり通常の筆算による乗算／繰り上がり加算を並べてみると対称性が破れて、高速アルゴリズムに十分な情報が得られるのか気になる
もちろん素朴な多項式乗算は、多項式の次数に対して遅い。とはいえ実際に100次多項式を2つ扱う必要があるのはどんな時だろうか？
こうした理由から、計算機代数システムではこのような方法は使われないという印象がある
https://www.youtube.com/watch?v=CcZf_7Fb4Us
https://en.wikipedia.org/wiki/Reed%E2%80%93Solomon_error_cor... が一例
[1]: https://github.com/8051enthusiast/delsum
Hazy Research の 2020〜2023 年のブログ記事には、このアプローチに関する情報が多い
「（…）物理学研究で、長さがほぼ 1 テラバイト、項が 1 億個を超える式を扱ったことがある」