IMQ kernelの可視化をする

最近、[A Visual Exploration of Gaussian Processes](https://www.jgoertler.com/visual-exploration-gaussian-processes/)というサイトを見ました。今まであまりよく理解できていなかったカーネル法やガウシアンRBFカーネルを視覚的に見ることができて理解が深まった気がします(理解したとはいっていない)。

これを真似すれば、[以前書いたWAEの記事](https://qiita.com/knok/items/ed6e829cf35aeae2d68b)での[IMQ関数](https://qiita.com/knok/items/ed6e829cf35aeae2d68b#inverse-multiquadric-kernel-imq)も理解が深まるのではないかと思い、やってみました。

# TL;DR

jupyter notebookで以下のような感じのグラフ可視化ができました。

![wae-imq-inter](https://knok.github.io/images/wae-imq-inter.gif)

コード全体は[github](https://github.com/knok/viz-kernel)に掲載してあります。

## RBFについて

自分はこれまでRadial basis kernel(RBF)について誤解していたのですが、RBFにも種類が複数あり一般にSVMなどで単に「RBFカーネル」として紹介されるのはGaussian kernelであるようです。

定義

```math
k(x, y) = \exp \left(- \frac{ ||x - y||^2}{2\sigma^2} \right)
```

他の書き方もあるようです。

```math
k(x, y) = \exp \left(- \gamma ||x - y||^2 \right)
```

これは単に $\gamma = \frac{1}{2\sigma^2}$ と表現しているだけで同じものだと理解できました。

* ref: http://crsouza.com/2010/03/17/kernel-functions-for-machine-learning-applications/#gaussian

## IMQについて

Inverse Multiquadric(IMQ)はRBFカーネルの一種で日本語では[逆多重二乗RBF](https://ja.wikipedia.org/wiki/%E6%94%BE%E5%B0%84%E5%9F%BA%E5%BA%95%E9%96%A2%E6%95%B0)という名称がついているようです。

こちらも異なる式での解説が複数見受けられます。Wikipediaでの記述に習いつつ、表記を合わせると以下になるでしょう。

```math
k(x, y) = \frac {1} {\sqrt{1 + \epsilon ||x - y||^2} }
```

先の[自分の記事から貼ったリンク(PDF)](https://qiita.com/knok/items/ed6e829cf35aeae2d68b#inverse-multiquadric-kernel-imq)でも同じ式が示されています。

一方で、[César Souzaの解説記事](http://crsouza.com/2010/03/17/kernel-functions-for-machine-learning-applications/#inverse_multiquadric)ではちょっと違った式になっています。

```math
k(x, y) = \frac{1}{\sqrt{||x - y||^2 + c^2 }}
```

$||x-y||$は頻出するので以下$r = ||x - y||$とします。分母の平方根の中身に着目すると、

```math
r^2 + c^2 = 1 + \epsilon r^2 \\
c^2 = 1 + \epsilon r^2 - r^2 \\
c^2 = 1 + ( \epsilon -1) ||x - y||^2
```

と考えれば内容は同じだと考えられます。ただし$c^2 > 0$なのでこの場合$\epsilon \geq 1$となるような気がします(間違っていたらご指摘ください)。

### WAE-MMDでの式

[Wasserstein Auto-Encoders](https://arxiv.org/abs/1711.01558)の論文内ではまた違う式で表現しています。

```math
k(x, y) = \frac{C}{C+||x - y||^2}
```

分母に平方根がありませんね…[Wikipediaの解説](https://ja.wikipedia.org/w/index.php?title=%E6%94%BE%E5%B0%84%E5%9F%BA%E5%BA%95%E9%96%A2%E6%95%B0&oldid=65814482#RBF%E3%81%AE%E7%A8%AE%E9%A1%9E)を見る限り、逆二乗カーネル(Inverse quadratic kernel)に見えます。
式変形をしてみましょう。

```math
\frac{C}{C+r^2} = \frac{1}{1+\frac{r^2}{C}} \\
( \epsilon r)^2 = \frac{r^2}{C} \\
\epsilon^2 r^2 = \frac{r^2}{C} \\
\epsilon^2 = \frac{1}{C} \\
\epsilon = \sqrt{\frac{1}{C}}
```

$r^2 = ||x - y||^2$なので

```math
k(x, y) = \frac{1}{1+(\epsilon r)^2}
```

と変換でき、よってこのカーネル関数はInverse Multiquadric kernelではなくInverse quadratic kernelが正しい気がします。

<!---
```math
\frac{C}{C + r^2} = \frac{1}{\sqrt{1 + \epsilon r^2}} \\
\frac{1}{1 + r^2/C} = \frac{1}{\sqrt{1 + \epsilon r^2}} \\
r^2/C = \sqrt{1 + \epsilon r^2} \\
(r^2/C)^2 = 1 + \epsilon r^2 \\
(r^2/C)^2 - \epsilon r^2 = 1 \\
(r^2)^2 - C^2 \epsilon r^2 = C^2 \\
C^2 + C^2 \epsilon r^2 = (r^2)^2 \\
C^2(1 + \epsilon r^2) = (r^2)^2 \\
C^2 = \frac{r^4}{1 + \epsilon r^2} \\
C = \frac{r^2}{ \sqrt{1 + \epsilon r^2}}
```

元の式も分子・分母をCで除算して

```math
\frac{1}{1 - r^2/C} = \frac{1}{1 - \frac{r^2}{\frac{r^2}{ \sqrt{1 + \epsilon r^2}}}} \\
= \frac{1}{1-\sqrt{1 + \epsilon r^2}}
 ```

…若干定数は異なりますが、性質としては近いように思います。ただし$\sqrt{1 + \epsilon r^2} < 1$を満たす必要がありそうです。
--->

なお、WAE論文の原文では"In all experiments we used $C = 2_z\sigma^2_d$, which is the expected squared distance between two multivariate Gaussian vectors drawn from $P_Z$"とあります。
Cは確率分布$P_Z$由来の多変量ガウスベクトル2点間の距離の2乗、ということのようです。

## ガウシアンカーネルの可視化

まずは"A Visual Exploration of Gaussian Process"のRBFガウシアンカーネルの可視化と同じような出力を出すことを考えてみましょう。
リアルタイムにインタラクティブな操作ができるところまでいきなり到達するのは難しそうなので、まずはjupyter notebookとmatplotlibで似たような図を出すことを考えます。

$k(x, y)$の値をそれぞれ-5 ~ 5の間で計算し、その値をヒートマップとして表示しているようなので、matplotlibのpcolormeshを使うことができそうです。

* 参考: [Python + matplotlib によるヒートマップ - Qiita](https://qiita.com/ynakayama/items/7dc01f45caf6d87a981b)

ガウシアンカーネルを関数として定義します。一般的には入力x, yはそれぞれベクトルですが、今回は可視化のためスカラ値を入力し、定義に沿って出力します。計算式は"A Visual Exploration of Gaussian Processes"に合わせてみます。

```python
def rbfgaussian(x, y, length=0.8, sigma=0.8):
    z = (sigma)**2 * np.exp(-((x - y)** 2/ (2 * length**2)))
    return z
```

単純にX, Y = [0, 1, 2, 3, 4, 5]としたときの図を書いてみます。

```python
X = list(range(5+1))
data = np.zeros((5+1, 5+1))
for x in X:
    for y in X:
        data[x, y] = rbfgaussian(x, y)
plt.pcolormesh(X, X, data, cmap='Blues')
plt.colorbar()
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/31714/0af111cc-e32c-24e3-24c9-30d75be495a1.png)

書きたい図と比べるとちょっと異なります。これはおそらく軸の方向がA Visual Exploration of Gaussian Processesとは逆になっているのではないかと思います。y軸を反転させてみましょう。

```python
X = list(range(5+1))
data = np.zeros((5+1, 5+1))
for x in X:
    for y in X:
        data[x, y] = rbfgaussian(x, y)
plt.gca().invert_yaxis() # Y軸反転
    plt.pcolormesh(X, X, data, cmap='Blues', vmin=0, vmax=1)
plt.colorbar()
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/31714/bdb95217-0f09-ee5b-812d-36e1c03be6ed.png)

それっぽくなってきました。
粒度を上げて、値の範囲も-5~5に変更し、描画するまでの処理を関数にしてみます。

```python
def drawheatmap(max_xy, step, kfunc, **kwargs):
    X = np.arange(-max_xy, max_xy+1, step)
    chunks = X.shape[0]
    data = np.zeros((chunks, chunks))
    X = np.arange(-max_xy, max_xy+1, step)
    for i, x in enumerate(X):
        for j, y in enumerate(X):
            data[i, j] = kfunc(x, y, **kwargs)
    plt.gca().invert_yaxis()
    plt.pcolormesh(X, X, data, cmap='Blues', vmin=0, vmax=1)
    plt.colorbar()
```

描画させると以下のようになります。

```python
drawheatmap(5, 0.5, rbfgaussian)
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/31714/679480a0-babf-40cf-0f92-c8c1e64ef4a5.png)

パラメータも変えてみましょう。sigmaを大きくすれば濃くなるはずです。

```python
drawheatmap(5, 0.5, rbfgaussian, sigma=1)
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/31714/3e85be49-262a-b15d-f6eb-fb6f1cf6ecb3.png)

## インタラクティブに表示させる

もとのページではパラメータ$\sigma, l$がスライダーで変更できるようになっていたので、同じことをやらせます。ipywidgetsを使えばできそうです。

[こちらの記事](https://qiita.com/tanemaki/items/e965cd9c25acbeed860c)を参考にしました。ただ、この記事では古いモジュール名(IPython.html.widgets)を使っているので、ipywidgetsに読み替えます。

```python
@interact(s=(0.1, 1), l=(0.1, 1))
def inter_rbf(s, l):
    drawheatmap(5, 0.5, rbfgaussian, sigma=s, length=l)
```

うっかり0除算しないように範囲を決める必要があります。

![rbfg-inter.gif](https://qiita-image-store.s3.amazonaws.com/0/31714/1fea3701-c4de-d010-e0ee-f10b0551951b.gif)

計算は基本jupyter notebookが動作している側でやることになりますが、件のページとだいたい同じようなことが再現できました。

## IMQ kernelの可視化

既に計算を行う部分を関数化して定義しているので、その中身をIMQに置き換えれば、同じことができます。

```python
def rbfimq(x, y, c=1):
    z = 1 / np.sqrt((x-y)**2 + c**2)
    return z
@interact(c=(0.01, 5))
def inter_rbf(c):
    drawheatmap(5, 0.5, rbfimq, c=c)
```

![imq-inter.gif](https://qiita-image-store.s3.amazonaws.com/0/31714/3ad9e17a-a577-2933-06a6-eb66bc6674aa.gif)

ガウシアンカーネルと比較すると、色の濃い(値の大きい)領域がより太くなっている傾向が見られます。また、Cを1より大きくすると全体的に薄くなっていくことが確認できます。

## WAE-MMDでの式

おそらくはIMQではなくInverse quadraticだと思われる、WAE論文の式を可視化してみます。

```python
def rbfwaeimq(x, y, c=1):
    z = c / (c + (x-y)**2)
    return z
@interact(c=(0.01, 5))
def inter_rbf(c):
    drawheatmap(5, 0.5, rbfwaeimq, c=c)
```

![wae-imq-inter](https://knok.github.io/images/wae-imq-inter.gif)

$x = y$の時は常に$C/C = 1$となるので、その領域は色が濃くなっています。Cを1より大きくすると色の濃いエリアが広がる傾向が見られます。

やはりWAEのカーネル関数はIMQではないのでは…との思いが強くなります。

## 線形カーネル

せっかくなので線形カーネルでもやってみましょう。

```python
def linear(x, y, c = 0, sigma=0.3, sigmab=0.8):
    z = sigmab**2 + sigma **2 * (x - c) * (y - c)
    return z

@interact(c=(-1.5, 1.5), sigma=(0.01, 1), sigmab=(0.1, 1))
def inter_rbf(c, sigma=0.34, sigmab=0.8):
    drawheatmap(5, 0.5, linear, c=c, sigma=sigma, sigmab=sigmab)
```

![linear-inter](https://knok.github.io/images/linear-inter.gif)

## 最後に

可視化をしてみることで理解を深めることができました。世の中の解説は微妙に式が違うことがあり、式変形を自分でやってみることも重要である、という気づきを個人的には得られました。

しかしどうせなら、A Visual Exploration of Gaussian Processesと同じようにすべてクライアントサイドで完結するような方法でも実装してみたいものです。

## 参考

* [A Visual Exploration of Gaussian Processes](https://www.jgoertler.com/visual-exploration-gaussian-processes/)
* [Kernel Functions for Machine Learning Applications – César Souza](http://crsouza.com/2010/03/17/kernel-functions-for-machine-learning-applications/)
* [[1711.01558] Wasserstein Auto-Encoders](https://arxiv.org/abs/1711.01558)
* [Python + matplotlib によるヒートマップ - Qiita](https://qiita.com/ynakayama/items/7dc01f45caf6d87a981b)
* [IPythonでクリクリできる図を作ろう - Qiita](https://qiita.com/tanemaki/items/e965cd9c25acbeed860c)
* [WAE-MMDをChainerで実装 - Qiita](https://qiita.com/knok/items/ed6e829cf35aeae2d68b#inverse-multiquadric-kernel-imq)