DFTを行列の対角化の視点から捉える

はじめに

行列の \(n\) 乗の計算をする時,対角化すると,固有値を \(n\) 乗するだけになるので計算が楽になります (競プロだと繰り返し二乗法で十分速いのであまり使い所がないですが,手計算でやるときは対角化をすると圧倒的に楽になります).

ところで,数列の \(n\) 乗を計算する時,数列を DFT してから,各要素を \(n\) 乗し,それをIDFT すれば \(O(N\log N)\) で高速に求まります.

これって行列の対角化に似てませんか?実は DFT は行列の対角化として捉えることができるのです.

DFT のおさらい

以下の記事が非常によくまとまっているのでそちらを見るのがおすすめですが,一応こちらにもかんたんに定義や基本的な性質を書いておきます.

qiita.com


数列 \((a_n)_{n=0,\dots,N-1}\) の離散 Fourier 変換 (DFT) \((A_n)_{n=0,\dots,N-1}\) は,\(W_N = \exp\left(-i\frac{2\pi}{N} \right)\) (1の\(N\)乗根) とすると,次の式で定義されます.

$$
A_n = \sum_{k=0}^{N-1} a_k W_N^{k}
$$

以降もとの数列を小文字で,それの DFT を大文字で書きます.また数列の長さは全部 \(N\) とします.

また,通常の畳み込みの代わりに巡回畳み込みを考えることにします.これは,添字を \(\mathrm{mod}\,N\) で考えて,畳み込み後の数列の長さが変わらないようにするだけです.つまり,数列 \(a, b\) の巡回畳み込み \(c\) は次式で定義されます.

$$
c_n = \sum_{k=0}^{N-1} a_{n-k\,\mathrm{mod}\,N} b_k \quad \text{for } n = 0,\dots, N-1
$$

通常の畳み込みを計算したいときは,適当に数列を延長して0を詰めて巡回畳み込みを計算すれば良いです.

DFT の著しい性質は,巡回畳み込みが要素ごとの積に置き換わることです.つまり,\(C_n = A_n B_n \) が成り立ちます.実際,

$$
C_n = \sum_{k=0}^{N-1} c_k W_N^{nk} \\
= \sum_{k=0}^{N-1} \left( \sum_{l=0}^{N-1} a_{k-l} b_l \right) W_N^{nk} \\
= \sum_{l=0}^{N-1} \left( \sum_{k=0}^{N-1} a_{k-l}W_N^{n(k-l)} \right) b_l W_N^{nl} \\
= \sum_{l=0}^{N-1} \left( \sum_{k=0}^{N-1} a_k W_N^{n(k-l)} \right) b_l W_N^{nl} \\
= \left( \sum_{k=0}^{N-1} a_k W_N^{nk} \right) \left( \sum_{l=0}^{N-1} b_l W_N^{nl} \right) \\
= A_n B_n
$$

3行目の変形では,和を取る順番を入れ替えて,因子 \(W_N^{nk}\) を分解しました.4行目の変形では,\(k-l \rightarrow k\) と置き換えました (添字を \(\mathrm{mod}\,N\) で考えており,さらに \(W_N^N=1\) なのでこれができる)

要素ごとの積は \(O(N)\) で求まるので,DFT が高速に計算できれば数列の畳み込みが高速に計算できます.DFT を定義どおりに計算すると \(O(N^2)\) となり,愚直に畳み込みを計算するのと変わりませんが,DFT を分割統治で高速に \(O(N\log N)\) で計算する高速 Fourier 変換 (FFT) というアルゴリズムが知られており,それを用いると畳み込みが \(O(N\log N)\) で求まります.

数列を線形写像として表す

DFTを行列の対角化として捉えたいので,まず数列を行列で表すことから始めます.数列を\(N\) 次元ベクトルとして捉えます.数列 \(x \in \mathbb{R}^N \) に,数列 \(a \in \mathbb{R}^N \) との巡回畳込み \(ax \in \mathbb{R}^N\) を対応させる写像を \(f_a(x)\) とします. これは線形写像です.

\(f_a\) の表現行列 \(M_a\) は以下のように表せます.

$$
M_a = \begin{pmatrix}
a_0 & a_{N-1} & \dots & a_1 \\
a_1 & a_0 & \dots & a_2 \\
\vdots & \vdots & \ddots & \vdots \\
a_{N-1} & a_{N-2} & \dots & a_{0}
\end{pmatrix}
$$

\(M_{ab}=M_aM_b, \quad M_{a^n} = M_a^n\) です.これで,数列の畳み込みが行列の掛け算として表現できました.

表現行列の対角化

\(M_a\) を対角化してみましょう.よく見ると \(M_a\) は巡回行列になっているので,この構造を利用すると,\(M_a\) の固有ベクトルは \(\boldsymbol{v}_n = \left( W_N^{n0}, W_N^{n1}, \dots, W_N^{n(N-1)} \right)^\intercal \text{for } n=0,\dots, N-1 \) となり,\(v_n\) に対応する固有値は \(A_n = \sum_{k=0}^{N-1} a_k W_N^{nk}\) になっていることがわかります.

これはDFTになっていますね?

よって,
$$
W = \begin{pmatrix}
W_N^{00} & W_N^{10} & \dots & W_N^{(N-1)0} \\
W_N^{01} & W_N^{11} & \dots & W_N^{(N-1)1} \\
\vdots & \vdots & \ddots & \vdots \\
W_N^{0(N-1)} & W_N^{1(N-1)} & \dots & W_N^{(N-1)(N-1)}
\end{pmatrix}
$$
とすれば,

$$
M_a = W\mathrm{diag}\left( A_0, A_1, \dots, A_{N-1} \right) W^{-1}
$$

と対角化できるわけです.

畳み込みの高速計算

\(a^n\) の表現行列 \(M_{a^n}\) は \(M_a^n\) として表せました.\(M_a\) が対角化できたので,その \(n\) 乗は固有値を \(n\) 乗するだけで簡単に求まります.

$$
M_{a^n} = M_a^n = W\mathrm{diag}\left( A_0^n, A_1^n, \dots, A_{N-1}^n \right)W^{-1}
$$

これは,DFTして,各要素を \(n\) 乗して,IDFTするという操作に対応しています.これで,DFTを用いて数列の \(n\) 乗を計算するアルゴリズムが,行列を対角化して固有値を \(n\) 乗するというアルゴリズムの特別な場合であることがわかりました.

ところで,対角化に用いた行列 \(W\) は,数列 \(a\) に依存しない形になっています.つまり,\(M_x\) は \(x\) が何であろうと全部 \(W\) で対角化できるということです.これすごいですね.

\(M_a = W\mathrm{diag}\left( A_0, A_1, \dots, A_{N-1} \right)W^{-1},\,M_b = W\mathrm{diag}\left( B_0, B_1, \dots, B_{N-1} \right) W^{-1}\) ですから,

$$
M_{ab} = M_aM_b = W\mathrm{diag}\left( A_0B_0, A_1B_1, \dots, A_{N-1}B_{N-1} \right) W^{-1}
$$

と計算できます.これは,DFT して要素ごとの積をとってIDFTするという操作に対応しています.一般の行列では,同じ行列で対角化できるとは限りませんから,行列積も対角化すると高速に計算できるというのは巡回行列の間で成り立つ特別な性質です.

おわりに

DFTすると畳み込みが積になるというのは,式は追えるんだけどあまり直観的な理解ができないでいました.今回,行列のDFTは対角化に対応しているという解釈を知って,少しは見通しが良くなった感覚があります.まあ直観的な理解に至ったかと言われるとまだ微妙ですが.

Hadamard 変換もこんな感じで解釈できるのでしょうか.