周期カーネルに対するFeature Map

ガウス過程回帰の課題のひとつとしてスケーラビリティがあります。
データ数をNとしてfitting の時間計算量 O(N^3)、predict の時間計算量が O(N^2)であることから大量のデータを扱いたいときは何らかの近似が必要です。

そのような近似手法のひとつとして、カーネルを有限次元の特徴ベクトルの内積で近似してベイズ線形回帰に帰着するものがあります。
特徴量ベクトルの次元をDとして、ベイズ線形回帰のfitting の時間計算量は O(D^3) 、predict の時間計算量は O(D)なので特徴量ベクトルの次元がデータ数よりも小さければより高速に計算ができることになります。

与えられたカーネルに対応する特徴量ベクトルを作る方法としてRandom Fourier Features
があります。*1
Gaussian カーネルやLaplacian カーネルについてはRandom Feature Map で特徴量ベクトルを得られますが、周期カーネルにこの方法を適用するのは難しそうでした。

そこでいろいろ調べたところ、周期カーネルを近似する論文*2を見つけたので概要と実装を備忘録としてまとめておきます。

(追記: 一部の数式が上手く表示されないようなので画像で置き換えました。何もわからない。)

Random Fourier Features

そもそもRandom Fourier Features がどのようなアルゴリズムなのか簡単に述べておきます。
 \mathbb{R}^{D}上のstationaryなカーネル k(\mathbf{x}, \mathbf{y}) = k(\mathbf{x} - \mathbf{y})がpositive definite であることと、ある有限非負測度のフーリエ変換で書けることが同値である(Bochner の定理)ことを使います。

Bocher の定理より、適当に定数倍したpositive definite なカーネルはある確率密度関数 p(\mathbf{\omega})フーリエ変換で書けます。
\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}
最後の式において、 \mathbf{\omega} p(\mathbf{\omega})に従い、bは [0, 2\pi]の一様分布に従います。

 \mathbf{\omega}と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}
ここで\omega_{0} = \frac{2\pi}{T}で、このカーネルは周期Tの関数を説明したい時などに使われます。

Random Fourier Features で特徴量ベクトルを作るなら、周期カーネルを逆フーリエ変換して得られる確率密度関数に従って周波数をサンプリングすることになりますが、どうみても簡単にサンプリングができる確率分布を得られそうにありません。
そこで、周期カーネルフーリエ級数展開して近似することを考えます。

周期カーネルフーリエ級数展開を次のように書きます。
f:id:lan496:20180716013606p:plain

k_{perSE}(\tau)は偶関数であることと、第一種変形ベッセル関数の積分表示I_{n}(x) = \frac{1}{\pi} \int_{0}^{\pi} e^{x \cos \theta} \cos n \theta\, d\theta \, (n \in \mathbb{Z})
*5を使うと、 c_{k}は次のように計算できます。
f:id:lan496:20180716013619p:plain

フーリエ級数の第K項以降を切り捨てて近似します。
f:id:lan496:20180716013630p:plain
よって、周期カーネルを2K次元の特徴量ベクトル \mathbf{\phi}(t)で近似できることがわかりました。

実装とテスト

論文の公式実装が見当たらなかったので自分で書きました。*6

論文ではランダムに2000個の入力に対して、周期T=2の周期カーネルでグラム行列 \mathbf{K}を計算したものと、第C項までを使って近似で得られる行列 \tilde{\mathbf{K}}のフロベニウスノルムでの誤差 \frac{||  \tilde{\mathbf{K}} -
 \mathbf{K}||_{F}}{|| \mathbf{K} ||_{F}} を測っていたので、自分の実装で検証してみました。

特徴量ベクトルの次元とハイパーパラメターlength_scale(l)を変えた時の誤差は下図のとおりです。
図の左上部分の誤差が元論文と異なりますが、傾向としては同じ結果が得られました。
f:id:lan496:20180716011413p:plain

下図はlength_scale(l)=0.2の場合にカーネルの値をプロットしたものです。
フーリエ級数展開の16次まで取り込めば、元の周期カーネルをほぼ再現できていることがわかります。
f:id:lan496:20180716012425p:plain