TFHEでの多項式乗算におけるFFTの利用

2020/9/27 式が一部間違っていることに気づいたので修正

# TFHEでの多項式乗算
TFHEにおけるFFTの使用目的は、Bootstrapping処理において発生するTRGSWとTRLWEの積の計算に必要になる多項式剰余環上の乗算を計算することです。

```math
\begin{aligned}
&x:多項式の変数\\
&a_n,b_n \in \mathbb{R}:多項式の係数\\
&N \in \mathbb{N}:多項式の係数の個数\\
&A = \sum^{N-1}_{n = 0} a_n x^n\\
&B = \sum^{N-1}_{n = 0} b_n x^n
\end{aligned}
```

と定義したとき、目的は以下のような演算を行うことです。

```math
\begin{aligned}
c_n &\in \mathbb{R}\\
C &= \sum^{N-1}_{n = 0} c_n x^n = A * B \pmod{x^N+1} = \sum^{N-1}_{n = 0}\sum^{N-1}_{m = 0} a_nb_m x^{n+m} \pmod{x^N+1} = \sum^{N-1}_{l = 0}(\sum_{m+n = l}a_nb_m - \sum_{m+n = l+N}a_nb_m) x^l
\end{aligned}
```

この式で重要なことは、一番最後で$\pmod{x^N+1}$が消えることです。つまり、剰余を計算しないで済みます。

## FFTでの乗算

前提知識として、DFT(実際にはFFTで計算します)を用いた乗算のアルゴリズムを確認します。多項式を係数の配列としてみなします。
$A$と$B$の係数をDFTにそれぞれかけて、かけ合わせ、IDFTをしたものが、$A * B \pmod{X^N - 1}$と一致することを確認します。

```math
\sum^{N-1}_{k = 0}\frac{1}{N}[\sum^{N-1}_{l=0}(\sum^{N-1}_{n=0}a_ne^{-i\frac{2\pi n l}{N}} \cdot \sum^{N-1}_{m=0}b_me^{-i\frac{2\pi m l}{N}})e^{i\frac{2\pi k l}{N}}]x^k
```

ここで、

```math 
\sum^{N-1}_{n=0} e^{\frac{2\pi l n}{N}} = \begin{cases}N\ if\ l\equiv 0\bmod{N}\\ 0\ otherwise\end{cases}
```
であることを用いると(逆離散フーリエ変換をすれば容易に確かめらます)、上の式は以下のようになります。

```math
\sum^{N-1}_{k = 0}(\sum_{n+m = k}a_nb_m + \sum_{n+m = k+N}a_nb_m) x^k = A*B\pmod{x^N-1}
```

一般的な多倍長整数の乗算での使用の場合などでは、$A$と$B$の下半分にのみ要素を詰め、上の要素を0で埋めるのが最も原理的な使い方です。

## TFHEでSPQLIOSを用いているときのアルゴリズム

SPQLIOSの場合は単純なDFTではなく、一般にDWT(WはWighted、日本語では離散荷重変換)と呼ばれる処理を用いています。また、実数FFTを用いるのではなく入力を複素数として虚数部にも要素を詰めるというアプローチをとっています。端的に述べれば以下のような計算を行います。

```math
\begin{aligned}
\sum^{\frac{N}{2}-1}_{l = 0}&\frac{1}{\frac{N}{2}}e^{-i\frac{2\pi l}{2N}}\{\sum^{\frac{N}{2}-1}_{m=0}[\sum^{\frac{N}{2}-1}_{n=0}e^{i\frac{2\pi n}{2N}}(a_n+ia_{\frac{N}{2}+n})e^{-i\frac{2\pi l n}{\frac{N}{2}}} \cdot \sum^{\frac{N}{2}-1}_{n=0}e^{i\frac{2\pi n}{2N}}(b_n+ib_{\frac{N}{2}+n})e^{-i\frac{2\pi l n}{\frac{N}{2}}}]e^{i\frac{2\pi l m}{\frac{N}{2}}}\}x^l\\
&= \sum^{\frac{N}{2}-1}_{l = 0}[e^{-i\frac{2\pi l}{2N}}\sum_{n+m=l}e^{i\frac{2\pi l}{2N}}(a_n+ia_{\frac{N}{2}+n})(b_m+ib_{\frac{N}{2}+m}) + \sum_{n+m=l+\frac{N}{2}}e^{i\frac{2\pi (l+\frac{N}{2})}{2N}}(a_n+ia_{\frac{N}{2}+n})(b_m+ib_{\frac{N}{2}+m})]x^l\\
&= \sum^{\frac{N}{2}-1}_{l = 0}[\sum_{n+m=l}(a_n+ia_{\frac{N}{2}+n})(b_m+ib_{\frac{N}{2}+m}) + \sum_{n+m=l+\frac{N}{2}}i(a_n+ia_{\frac{N}{2}+n})(b_m+ib_{\frac{N}{2}+m})]x^l \\
&= \sum^{\frac{N}{2}-1}_{l = 0}\{[\sum_{n+m=l}(a_nb_m-a_{\frac{N}{2}+n}b_{\frac{N}{2}+m}) - \sum_{n+m=l+\frac{N}{2}}(a_nb_{\frac{N}{2}+m}+a_{\frac{N}{2}+n}b_m)]+i[\sum_{n+m=l}(a_nb_{\frac{N}{2}+m}+a_{\frac{N}{2}+n}b_n) + \sum_{n+m=l+\frac{N}{2}}(a_nb_m-a_{\frac{N}{2}+n}b_{\frac{N}{2}+m})]\}x^l
\end{aligned}
```

よって、計算結果の実数部が$A*B\pmod{x^N-1}$の下半分に、虚数部が上半分になっています。
この方式は$\frac{N}{2}$要素のFFTしか必要としません。

## TFHEでFFTWを用いているときのアルゴリズム

おまけで書きますが、この方式は非効率です。原論文者がなぜ次のSPQLIOSを用いた場合のアルゴリズムと同じになるような形でFFTWを用いなかったのかは謎です。
まず、多項式乗算自体は以下のように、多項式の係数の配列を2N個に拡張して、下半分に多項式の係数の配列をそのまま詰め、上半分に多項式の係数の符号を反転したものを入れて考えます。

```math
\begin{aligned}
f_n &= \begin{cases}a_n\ if\ n< N\\ -a_{n-N}\ otherwise\end{cases} \\
g_n &= \begin{cases}b_n\ if\ n< N\\ -b_{n-N}\ otherwise\end{cases} \\
F*G&=\sum^{2N-1}_{l = 0}\frac{1}{N}[\sum^{2N-1}_{m=0}(\sum^{2N-1}_{n=0}f_ne^{-i\frac{2\pi l n}{2N}} \cdot \sum^{2N-1}_{n=0}g_ne^{-i\frac{2\pi l n}{2N}})e^{i\frac{2\pi l m}{2N}}]x^l= \sum^{2N-1}_{l = 0}(\sum_{n+m = l}f_ng_m + \sum_{n+m = l+N}f_ng_m) x^l
\end{aligned}
```

$l\leq N-1$の場合にのみ興味があるので、その場合に限ります。

```math
\sum_{n+m = l}f_ng_m + \sum_{n+m = l+N}f_ng_m = \sum_{n=0}^{l}a_nb_{l-n} - \sum_{n=l+1}^{N-1}a_nb_{l-n+N} + \sum_{n=N}^{N+l}a_{n-N}b_{l-n+N} - \sum_{n=N+l+1}^{2N-1}a_{n-N}b_{l-n+2N} = 2(\sum_{n+m=l}a_nb_m - \sum_{n+m=l+N}a_nb_m)
```

よって、$F*G$の下半分の成分は$A*B\pmod{x^N+1}$の2倍に一致することがわかります。実際の実装では$f_n$と$g_n$を2で先に割っています。
TFHEでは上記に加えて対称性を利用した２つの最適化を施しています。1つ目は実数FFTに一般的なもので、実数FFTの結果は対象性が存在するため、実数の個数の半分の計算結果の複素数を計算すれば十分であるというものです。もう1つは$F$と$G$は上部と下部が符号が違うだけなためにFFT結果の偶数成分が全て0になるため要素ごとの乗算のときに計算しなくてよいというものです。
後者のみ示します。$l=2p$とします。

```math
\sum^{2N-1}_{n=0}f_ne^{-i\frac{2\pi l n}{2N}} = \sum^{2N-1}_{n=0}f_ne^{-i\frac{2\pi p n}{N}} = \sum^{N-1}_{n=0}(f_n + f_{N+n})e^{-i\frac{2\pi p n}{N}} = 0
```

この方法の非効率な点は、最も負荷の大きいFFTの要素数が先の方法に比べて倍になっていることです。先に述べたように実数FFTでは半分の個数のFFTの結果を計算すれば良いので概ね$N$要素のFFTと等価ですが、2つ目の最適化がいうようにその半分は自明に0であるので、理想的には$\frac{N}{2}$要素のFFTを計算すれば十分なわけで、実際SPQLIOSはそうなっています。

#参考文献
https://math.stackexchange.com/questions/1435448/negacyclic-fft-multiplication
https://nufhe.readthedocs.io/en/latest/implementation_details.html
多数桁計算における高速アルゴリズムの研究, 後保範, 早稲田大学大学院 博士論文