周期カーネルに対するFeature Map
ガウス過程回帰の課題のひとつとしてスケーラビリティがあります。
データ数をNとしてfitting の時間計算量、predict の時間計算量がであることから大量のデータを扱いたいときは何らかの近似が必要です。
そのような近似手法のひとつとして、カーネルを有限次元の特徴ベクトルの内積で近似してベイズ線形回帰に帰着するものがあります。
特徴量ベクトルの次元をDとして、ベイズ線形回帰のfitting の時間計算量は 、predict の時間計算量はなので特徴量ベクトルの次元がデータ数よりも小さければより高速に計算ができることになります。
与えられたカーネルに対応する特徴量ベクトルを作る方法としてRandom Fourier Features
があります。*1
Gaussian カーネルやLaplacian カーネルについてはRandom Feature Map で特徴量ベクトルを得られますが、周期カーネルにこの方法を適用するのは難しそうでした。
そこでいろいろ調べたところ、周期カーネルを近似する論文*2を見つけたので概要と実装を備忘録としてまとめておきます。
(追記: 一部の数式が上手く表示されないようなので画像で置き換えました。何もわからない。)
Random Fourier Features
そもそもRandom Fourier Features がどのようなアルゴリズムなのか簡単に述べておきます。
上のstationaryなカーネルがpositive definite であることと、ある有限非負測度のフーリエ変換で書けることが同値である(Bochner の定理)ことを使います。
Bocher の定理より、適当に定数倍したpositive definite なカーネルはある確率密度関数のフーリエ変換で書けます。
\begin{align}
k(\mathbf{x} - \mathbf{y})
&= \int_{\mathbb{R}^{D}} p(\mathbf{\omega}) e^{i (\mathbf{x} - \mathbf{y})^{T} \mathbf{\omega}} d\mathbf{\omega} \\
&= \mathbb{E}_{\mathbf{\omega}, b} \left[ z_{\mathbf{\omega}}(\mathbf{x}) z_{\mathbf{\omega}}(\mathbf{y})\right]
\quad (\mbox{where} \quad z_{\mathbf{\omega}}(\mathbf{x}) := \sqrt{2} \cos(\mathbf{x}^T \mathbf{\omega} + b))
\end{align}
最後の式において、はに従い、bは]の一様分布に従います。
とbをそれぞれD個ずつサンプリングして
\begin{align}
\mathbf{z}(\mathbf{x}) := \sqrt{\frac{2}{D}} \left( \cos(\mathbf{x}^T \mathbf{\omega}_{1} + b_{1}), \cdots, \cos(\mathbf{x}^T \mathbf{\omega}_{D} + b_{D})\right)^{T}
\end{align}
とすると、カーネルをモンテカルロ積分の要領で近似できます。
\begin{align}
k(\mathbf{x} - \mathbf{y})
&\simeq \frac{1}{D} \sum_{j = 1}^{D} z_{\mathbf{\omega}_{j}}(\mathbf{x}) z_{\mathbf{\omega}_{j}}(\mathbf{y}) \\
&= \mathbf{z}(\mathbf{x})^{T} \mathbf{z}(\mathbf{y})
\end{align}
例えばGaussian カーネルは正規分布をフーリエ変換で書けるのでこの方法で対応する特徴量ベクトルを得ることができます。
これは既にsklearn に実装されています。*3
Fourier Features for Periodic Kernels
次に周期カーネルを考えます。
具体的には次のようなカーネルです。*4
\begin{align}
k_{perSE}(t, t') = k_{perSE}(t, t')
&= \exp \left( -2 \left( \frac{\sin (\omega_{0} (t - t'))}{l} \right)^{2} \right) \\
&= \exp \left( \frac{\cos \omega_{0} (t - t') - 1}{l^{2}} \right)
\end{align}
ここでで、このカーネルは周期Tの関数を説明したい時などに使われます。
Random Fourier Features で特徴量ベクトルを作るなら、周期カーネルを逆フーリエ変換して得られる確率密度関数に従って周波数をサンプリングすることになりますが、どうみても簡単にサンプリングができる確率分布を得られそうにありません。
そこで、周期カーネルをフーリエ級数展開して近似することを考えます。
は偶関数であることと、第一種変形ベッセル関数の積分表示
*5を使うと、は次のように計算できます。
フーリエ級数の第K項以降を切り捨てて近似します。
よって、周期カーネルを2K次元の特徴量ベクトルで近似できることがわかりました。
実装とテスト
論文の公式実装が見当たらなかったので自分で書きました。*6
論文ではランダムに2000個の入力に対して、周期T=2の周期カーネルでグラム行列を計算したものと、第C項までを使って近似で得られる行列のフロベニウスノルムでの誤差を測っていたので、自分の実装で検証してみました。
特徴量ベクトルの次元とハイパーパラメターlength_scale(l)を変えた時の誤差は下図のとおりです。
図の左上部分の誤差が元論文と異なりますが、傾向としては同じ結果が得られました。
下図はlength_scale(l)=0.2の場合にカーネルの値をプロットしたものです。
フーリエ級数展開の16次まで取り込めば、元の周期カーネルをほぼ再現できていることがわかります。
*1:https://people.eecs.berkeley.edu/~brecht/papers/07.rah.rec.nips.pdf
*2:http://web.it.usyd.edu.au/~framos/Publications_files/Tompkins-Ramos.pdf
*3:sklearn.kernel_approximation.RBFSampler — scikit-learn 0.19.1 documentation
*4:sklearn.gaussian_process.kernels.ExpSineSquared — scikit-learn 0.19.1 documentation
*5:Modified Bessel Function of the First Kind -- from Wolfram MathWorld