今回はあまり知られていない機械学習の手法である教師なしカーネル回帰（Unsupervised Kernel Regression: UKR）¹という手法についてご紹介したいと思います。このUKRですが、実は昨今流行の変分自己符号化器（Variational Autoencoder: VAE)や敵対的生成ネットワーク（Generative Adversarial Nets: GAN)といったいわゆる深層生成モデル（deep generative models）の仲間であると解釈することができます。UKRは非常にシンプルな仕組みをしていて、VAEやGANとはまた別の良さがある手法なので、ぜひこちらの記事を読んで興味を持っていただければと思います。

暗黙的生成モデルとは？

UKRがどのような問題を解いているか説明するために、Implicit generative models²（以降、暗黙的生成モデルと呼ぶ）について説明します。大まかに言うと暗黙的生成モデルは確率分布の表現アプローチの1つで、確率分布をあからさまに指定することを避けることで表現力の向上を狙っています。深層生成モデルと呼ばれるVAEやGAN、そして本記事で注目するUKRはこの暗黙的生成モデルを採用し、それをデータに当てはめるアルゴリズムを示しています。まずは基本的な確率分布の推定問題から説明していきます。

統計的推論とモデルの確率密度関数表現

観測されたデータセット$X=\{x_i\}, x_i\in\mathcal{X}=\mathbb{R}^D$からそれを生成した真の確率分布$\tilde{p}(x)$を知りたいということがあります。ですが、当然のことながら真の$\tilde{p}(x)$を知る事はできません。ですので、機械学習や統計では、様々な確率分布を表現できるモデル$p(x)$をユーザが仮定し、$p(x)$を観測データ$X$に当てはめる³事で$\tilde{p}(x)$に近づけるということをやるわけです。このような問題は統計的推論と呼ばれます。
このときにまず考えなければならないのは、どんなモデル$p(x)$を仮定するか？ということです。よく用いられるのはデータとは独立したパラメータ$\theta$を持つ（パラメトリックな）確率密度関数を指定するアプローチです。例えば正規分布を用いるのであれば
$$p(x)=\mathcal{N}(x\mid\mu, \Sigma)$$
となります。この場合パラメータは平均と共分散行列、つまり$\theta=(\mu, \Sigma)$となります。あとはデータ$X$にうまく当てはまるような$\theta$を推定することになります。
また別のやり方としてはデータそのもので確率密度関数を表現する方法（ノンパラメトリック）があります。代表的な手法としては以下のようなカーネル密度推定が挙げられます⁴ 。
$$p(x)=\frac{1}{n}\sum_i \mathcal{N}(x \mid x_i, \Sigma)$$
パラメトリックに表現する場合、最終的に問題をパラメータの$\theta$の推定問題に落とし込めるので扱いやすくはなりますが、表現できる$p(x)$がシンプルすぎるのが悩ましいところです。逆にノンパラメトリックの場合はデータを用いて表現するので柔軟性があり、データが多いほど表現力は増します。ですが$x$が高次元になると莫大な数のデータを用意しない限りうまく働きません⁵。

暗黙的生成モデル

パラメトリックにしろノンパラメトリックにしろ、ここまで紹介したのは確率密度関数を直接指定するアプローチです。このアプローチだと高次元なデータに対して柔軟性の高い確率モデルを構築することは難しいかもしれません。そこで近年注目されているのが暗黙的生成モデルです。

暗黙的生成モデルでは、$p(x)$を定義するのに確率密度関数といった解析的な形で指定する（明示的）のではなく、$p(x)$から得られるサンプル$x$の生成過程を定義する（暗黙的）ということを考えます。その生成過程を図で表現すると以下のような感じです。

確率変数は円で囲んでいます。観測できない変数はグレーの円、観測できる変数は白い円で示しています。$z\in\mathcal{Z}=\mathbb{R}^L$は観測できない変数なので潜在変数と呼ぶことにします。この潜在変数$z$は単純な確率分布$p(z)$に従って生成されるとします。例えば
$$p(z)=\mathcal{N}(z\mid 0, I_L)$$
のように平均$0$、共分散行列が単位行列$I_L$となるような正規分布です。
この時、$x$は生成された$z$と潜在変数の空間$\mathcal{Z}$からデータの空間$\mathcal{X}$への写像$f:\mathcal{Z}\rightarrow\mathcal{X}$により次のように生成されるとします。
$$x=f(z)$$

この時、$z$が確率変数なので$x$も確率変数になります。つまり
$$x\sim p(x)$$
と$x$は何らかの確率分布$p(x)$に従っています。この場合、具体的な$p(x)$の確率密度を知る事はできませんが、$p(z)$から$z$を生成し、それを$f$で写像することで、この$p(x)$からのサンプルは得ることができます⁶。つまり、$f$と$p(z)$によって$p(x)$が間接的に表現されているので"implicit"（暗黙的）と呼んでいるわけですね。この$f$をデータから獲得しようというのが暗黙的生成モデルのデータへの当てはめです⁷。$f$を深層ニューラルネットワークで表現するものは深層生成モデルと呼ばれます。万能関数近似能力を持つニューラルネットをモデルの構築に導入することで、確率密度函数を指定するような明示的なモデル化を超える表現力を持つことを狙っているわけです。代表的な手法としてはVAEやGANがあります。これらの学習の目的は、ニューラルネットのパラメータ$\theta$を推定することです。

教師なしカーネル回帰

本記事で取り扱うUKRは、深層ではない暗黙的生成モデルです。どのような手法かこれから紹介していきます。

写像をどう表現するか

暗黙的生成モデルでは$f$をどう表現するかということが肝になってきます。VAEやGANではこの$f$をdeep neural networkで構築し、データ$X$から$f$のパラメータ$\theta$を最適化します。これは$f$をパラメトリックに表現しているわけですが、UKRでは$f$をデータを用いてノンパラメトリックに表現します。

どうするかというと、潜在変数を考えます。暗黙的生成モデルでは潜在変数からデータが生成されると仮定したので、1つの観測データ$x_i$に対して、生成する元となった$z_i$というものを考えることができます。全てのデータ$X$についての潜在変数をまとめて$Z=\{z_i\}$とします。この時$f$を、この潜在変数集合$Z$とデータ$X$を用いて次のように表現します⁸。
$$f(z\mid Z)=\frac{1}{K(z)}\sum_i k(z,z_i)x_i$$
ただし$K(z)=\sum_i k(z, z_i), k(z,z')=\exp(-\frac{1}{2\sigma^2}|z-z'|^2)$と定義しています。また$\sigma$は写像$f$の滑らかさを調整するハイパーパラメータになります。つまり$f$を表現するパラメータとして、データに対応する潜在変数$Z$を導入しているわけです。データに対して適応的に$f$を表現しているので、ノンパラメトリックであると言えます。

UKRの学習

最適化する目的関数（学習の結論）

$f$をどう表現するかを定義し終えたので、これにより$q(x)$が表現できる確率分布の集合が暗黙に定義されたことになります。ここからは実際にデータに当てはめていくフェーズになります。これを便宜上、学習と呼ぶことにします。暗黙的生成モデルでは$f$のパラメータを学習するわけですが、UKRでは$f$がデータに対応する潜在変数$Z$をパラメータとして持つので、この$Z$をデータから推定するのがUKRの学習になるわけです。

UKRの学習は具体的にどのように行われるか、結論から言うと、$Z$についての目的関数$L(Z)$を以下のように定義し、これを最小化するような$Z$を勾配法によって推定します。
$$L(Z)=\frac{1}{n}\sum_i \left|x_i-f(z_i\mid Z)\right|^2+\lambda \sum_i|z_i|^2\ \ \ (1)$$
目的関数の中身を見てみると、一項目は二乗誤差です。VAEでも出てくるような、潜在変数からデータを復元した時の誤差を表現しています。二項目は正則化項です。これを入れずに誤差だけを小さくしようとすると、回帰における過学習のような状態になります。$\lambda$はその正則化の強さをコントロールするハイパーパラメータです。

経験分布とモデルのKL divergence最小化から目的関数を導出する

と、ここまでは原論文¹に忠実なnotationでのミニマムな説明でして、基本的にはこれで以上です。暗黙的生成モデルでは最終的に$f$を獲得できればOKなのでこれで十分です。ですが決定論的に問題を捉えているので、確率モデルを考えているというのがいまいちピンとこないと思います。ここからは確率的な枠組みでこの目的関数を導出してみます。統計プロパーな方には自明な話かも知れません。また、込み入った話が続くのでUKRの特徴だけ知りたいという方は次節「UKRの特徴」にスキップしていただいても構いません。

ここで推定したいのは$f$のパラメータである潜在変数$Z$です。$Z$によってモデルが表現する確率分布が決まるので、ここではその確率分布を$p_Z(x)$と書くことにします。この$p_Z(x)$が観測されたデータから構成される経験分布
$$p_{\mathrm{data}}(x)=\frac{1}{n}\sum_i \delta(x-x_i)$$
にできるだけ近づけるようにするのが$Z$の推定の方針です。ここで$\delta(x)$はディラックのデルタ関数です。$\int\delta(x)dx=1$でなおかつ$x\neq 0$のとき$\delta(x)=0$となるような関数です。また任意の関数$f(x)$と一緒に積分すると$\int f(x)\delta(x)dx=f(0)$となる性質があります。この後の計算でこの性質は使います。

ここで必要ないくつかの仮定を導入します。観測データ$X$に対する$Z$の尤度$p(X\mid Z)$を次のように定義します。

\begin{align}
p(X\mid Z)&=\prod_i p_Z(x_i\mid z_i) \\
&=\prod_i \mathcal{N}\left(x_i\mid f(z_i\mid Z),\beta^{-1}I_D\right) 
\end{align}

ここで$\beta$は精度パラメータで既知としておきます。$I_D$は$D\times D$の単位行列です。要するにデータ$x_i$は平均を$f(z_i\mid Z)$とする等方な正規分布に従って生成されるということです。最初の変形から$X$はi.i.dであるように見えますが、厳密には潜在変数$Z$が与えられた上での条件付き独立という形になっています。データ$x_i$に対応する$z_i$は$f$のパラメータという形で$x_j$にも影響を与える形になっています。一方で$Z$に対する事前分布$p(Z)$を次のように仮定します。
$$p(Z)=\prod_i \mathcal{N}(z_i\mid 0, \alpha^{-1}I_L)$$
$\alpha$は事前分布の精度パラメータ$I_L$は$L\times L$の単位行列です。

このように仮定すると、UKRが表現する確率分布$p_Z(x)$は
$$p_Z(x)=\int p_Z(x\mid z)p(z)dz$$
となります。ただし$z$と$x$の間には非線形性のある$f$が入っているのでこの積分計算は解析的に解けません。この$p_Z(x)$を経験分布$p_{\mathrm{data}}(x)$に近づけることを考えます。近づき度合いの指標としてここではKL divergenceを用いて
$$D_{KL}\left[p_{\mathrm{data}}(x)||p_Z(x)\right]=-\int p_{\mathrm{data}}(x)\ln\frac{p_Z(x)}{p_{\mathrm{data}}(x)}dx$$
を最小化することを考えます。これを変形すると

\begin{align}
D_{KL}\left[p_{\mathrm{data}}(x)||p_Z(x)\right]&=-\int p_{\mathrm{data}}(x)\ln\frac{p_Z(x)}{p_{\mathrm{data}}(x)}dx\\
&=-\int p_{\mathrm{data}}(x)\ln p_Z(x)dx + C\\
&=-\frac{1}{n}\sum_i \ln p_Z(x_i)+C\\
&=-\frac{1}{n}\ln p_Z(X)+C
\end{align}

となります。$Z$に関係しない項は定数として$C$にまとめています。従って、このKL divergenceの最小化は対数周辺尤度$\ln p_Z(X)=\ln\prod_i p_Z(x_i)$の最大化と等価になります。ここからの流れはVAEと似たような感じで⁹、$p_Z(X)$は解析的な形を得ることはできないので、この下限を取り、それを最小化することを考えます。対数周辺尤度$\ln p_Z(X)$を整理すると

\begin{align}
\ln p_Z(X)&=\ln p_Z(X)\int q(Z\mid \hat{Z})dZ\\
&=\int q(Z\mid \hat{Z})\ln\frac{p_Z(X)p(Z\mid X)}{p(Z\mid X)}dZ\\
&=\int q(Z\mid \hat{Z})\left(\ln\frac{q(Z\mid \hat{Z})}{p(Z|X)}+\ln \frac{p(Z,X)}{q(Z\mid \hat{Z})}\right)dZ\\
&=D_{KL}\left[q(Z\mid \hat{Z})||p(Z|X)\right]+L(\hat{Z})
\end{align}

ここで$p(Z\mid X)$は厳密な事後分布です。これは解析的に求めることができないので、近似する分布$q(Z\mid \hat{Z})$を導入しています。$q(Z\mid \hat{Z})$は以下のようなパラメータ$\hat{Z}=\{\hat{z_i}\}$を持つデルタ関数として定義します¹⁰。

\begin{equation}
q(Z\mid\hat{Z})=\prod_i\delta(z_i-\hat{z}_{i})
\end{equation}

整理した周辺尤度の第一項$D_{KL}\left[q(Z\mid \hat{Z})||p(Z|X)\right]$は非負なのでそれ以外をまとめて下限$L(\hat{Z})$としています。第一項が$0$、つまり事後分布$p(Z\mid X)$とそれを近似する$q(Z\mid \hat{Z})$が一致すれば周辺尤度と下限が一致します。変分下限を最大化するようにこの近似事後分布のパラメータ$\hat{Z}$を求めることになります。

さらに下限を整理すると

\begin{align}
L(\hat{Z})&=\int q(Z\mid\hat{Z})\ln \frac{p(Z, X)}{q(Z\mid\hat{Z})}dZ\\
&=\int q(Z\mid\hat{Z})\left(\ln p(X\mid Z)+\ln\frac{p(Z)}{q(Z\mid\hat{Z})}\right)dZ\\
&=\int q(Z\mid\hat{Z})\ln p(X\mid Z)dZ - D_{KL}\left[q(Z\mid\hat{Z})\|p(Z)\right]\\
&=\ln p(X\mid \hat{Z})+\ln p(\hat{Z})\\
&=-\frac{\beta}{2}\sum_i \left\|x_i-f(z_i\mid \hat{Z})\right\|^2-\frac{\alpha}{2}\left\|\hat{z}_i\right\|^2
\end{align}

となり、(1)と比較するとデルタ関数の近似事後分布を導入しているため$Z$が$\hat{Z}$という表記になりますが、$\beta=2/N, \alpha=2/\lambda$とおいて符号を反転させれば一致することが分かります。つまり元の目的関数(1)の最小化は周辺尤度の下限を最大化していることに相当することが分かりました。

また余談ですが、この目的関数の導出の流れはVAEとほぼ同じとなっています。VAEをご存知の方向けにVAEとUKRの違いを以下に列挙しておきます。

写像の表現方法。VAEはニューラルネット。UKRはNadaraya-Watson kernel estimator。
近似事後分布に仮定するパラメトリックな分布。VAEは正規分布でUKRはデルタ関数。言い換えればVAEは潜在変数を分布推定するのに対してUKRは点推定する。
パラメータを推定するアルゴリズム。VAEの場合は確率的勾配法を用いてエンコーダやデコーダのパラメータを推定することが前提になっている（オンライン学習もしくはセミバッチ学習が可能）。UKRは$f$をデータ集合とそれに対応する潜在変数集合$Z$でパラメタライズ（ノンパラメトリックに表現）するため、基本的に確率的でない勾配法を用いることになる（原則としてバッチ型学習になる）。
両者ともに、あるデータ$x_i$が別のデータ$x_j$の潜在変数$z_j$の推論に影響を与えるが、その影響の与え方が異なる。VAEの場合は個々のデータに対応する潜在変数$Z=\{z_i\}$を推定するのに確率的なエンコーダ$q_{\phi}(z\mid x)$を導入しており、全データに共通のエンコーダのパラメータ$\phi$を推定することで、$x_i$が$x_j$の潜在変数$z_j$の推定に影響を与えることになる。一方でUKRでは$x_j$に対する$z_j$を推定する際に用いる$f$のパラメータに潜在変数$z_i$が含まれるため、結果的に$x_i$から$z_i$を推定する際には$z_j$にも影響を与えることになる。

UKRの特徴

確率モデルの表現方法としては冒頭でも述べたように明示的に密度関数を定義するパターンと潜在変数の確率密度関数と写像で定義する暗黙的生成モデルがあるわけですが、それぞれパラメトリックとノンパラメトリックなものがあります。これを整理すると次のようになります。

	明示的に密度関数を定義	暗黙的生成モデル
パラメトリック	ガウス分布など	VAE, GAN
ノンパラメトリック	カーネル密度推定	UKR

ということで、UKRは暗黙的生成モデルでありながらノンパラメトリックという珍しい立ち位置の手法になります¹¹。
パラメトリックな暗黙的生成モデルと比べたUKRの特徴は、アルゴリズムが非常にシンプルになるということです。例えばVAEの場合は潜在変数$Z$¹²とエンコーダのパラメータと写像$f$を表現するデコーダのパラメータの両方を推定する必要があります。これに対してUKRの場合は潜在変数$Z$が決まると写像$f$が決まるので、アルゴリズム上は$Z$を推定すれば良いということになります。
どちらが良いかというのは場合によりけりだと思いますが¹³、後者の良さとしては最終的に解く問題が非常にシンプルなので、フルスクラッチでも容易に実装できるというのは利点としてあるかと思います。ですので、簡単なトイデータを学習させてみることで暗黙的生成モデルの学習という問題の癖や勘所を体験する、という入門的な使い方にはかなり良いのではないかと思います。
またこれもノンパラメトリックに由来する特徴として、ユーザが置かなければいけないハイパーパラメータや仮定が少ないという点があります¹⁴。UKRの場合はハイパーパラメータは正則化の強さ$\lambda$や潜在変数の次元$L$などかなり少数です。一方でVAEやGANの場合は$f$を表現するneural networkのアーキテクチャは任意なので、これはユーザが指定する必要があります。

実装例

有志の皆さんがPythonやJuliaで実装したコードを公開し、トイデータでのお試しまで行っていますので本記事でUKRに興味を持った方はぜひ読んでみてください。
Variants of UKRを読んで実装したよ
 Juliaで Unsupervised Kernel Regression 実装
 Juliaに入門して教師なしカーネル回帰を実装してみる

最後に

個人的にはこのUKRという手法のシンプルさはとても気に入っていて、その良さを言語化するのが難しいのではありますが、皆さんにご紹介してみた次第です。本記事で取り上げた暗黙的生成モデルの学習という問題は実はそう簡単ではなくて、潜在変数の次元$L$によって結果が大きく変わったり、$f$が複雑すぎると過学習¹⁵になったりします（この辺りもどこかでまたご紹介できると良いのですが）。その難しさを体験するのにUKRというのはピッタリな題材であると私は思います。フルスクラッチでの実装コードも充実していますので、皆さんぜひ実装して色々試してみてもらえたらとても嬉しいです。

Principal Surfaces from Unsupervised Kernel Regression https://ieeexplore.ieee.org/document/1471705 ↩
Learning in Implicit Generative Models https://arxiv.org/abs/1610.03483 ↩
もう少し厳密に言うと、$p(x)$が表現しうる確率分布の中で$X$を生成した「尤もらしさ」が高い特定の確率分布$\hat{p}(x)$を求めることになります。 ↩
厳密に言えばこの場合、正規分布の共分散行列は$\Sigma$はユーザが与えることになります。 ↩
ものすごく雑にうまく働かない理由を説明すると、次元が上がると空間を満たすために必要なデータ数が指数的に増加するため、高次元だとデータがある程度あってもスカスカになってしまいます。するとカーネル密度推定をしてもカーネル関数（例えば正規分布）とカーネル関数が重ならず、ほとんど経験分布と変わらない確率分布を推定してしまうのだと考えられます。ざっとググった限りだとこの論文はその辺りをちゃんと数学的に示してくれていそうなので興味のある方は読んでみてください。 ↩
$p(x)$の具体的な確率密度が求められない理由の１つは、VAEやGANでは潜在空間の次元をデータ空間の次元より小さく取ることです。深層生成モデルの一種であるnormalizing flowでは、潜在空間をデータ空間の次元と同一とし、加えて写像$f$は可逆であるという制約を加えることで、表現力を向上させつつ確率密度を求めることができます。 ↩
潜在変数からデータを生成する写像$f$をデータセットから推定しようという問題は実は古くからあります。主曲線（principal curve）の推定問題や次元削減の一部の手法ではかなり前から取り組まれていました。 ↩
$Z$と$X$が既知の場合の$f$の推定は回帰問題であると捉えることができます。この観点から言えば、$f$をこのように表現するという事は、Nadaraya-Watson kernel estimatorもしくはKernel Regressionと呼ばれる方法で回帰を解いていることになります。 ↩
下限の導出は色んな文献で行われていますが、今回は深層生成モデルを巡る旅(2): VAEを参考にさせていただきました。 ↩
近似分布をデルタ関数に指定しているので結局のところMAP推定と同じことですが、今回は経験分布を近似していることが分かりやすいようにこのような導出を行いました。 ↩
他にノンパラな暗黙的生成モデルと解釈できる手法としてはガウス過程潜在変数モデル（Gaussian Process Latent Variable Models: GPLVM）もあります。 ↩
厳密に言うと潜在変数そのものを推定しているのではなく、データから潜在変数を確率的に推定するエンコーダを導入し、そのパラメータを推定しています。逆にGANの場合は、学習の中でデータに対応する潜在変数集合$Z$は登場せず、潜在変数$z$は写像$f$に入力するただの乱数という扱いで、あくまで$f$のパラメータ$\theta$を獲得することに主眼が置かれています。完全に脱線ですが、このデータに対する$Z$を考慮しないのがモード崩壊の一因なのではないか…？などと素人ながら思っているのですが、この辺り詳しい方がいたらぜひ教えてください。 ↩
潜在変数と写像の両方を推定するタイプと潜在変数のみを推定するタイプ、この両者の違いがもたらす影響を厳密に検証して良し悪しを議論するのは難しいので、潜在変数のみ推定＝利点ということではありません。 ↩
本文では利点のように書いてしまっていますが、これも利点と捉えるか欠点と捉えるかは場面によりけりだと思います。明確な性能評価の指標があり、それに向けてゴリゴリとチューニングできるのであればパラメトリックなVAEやGANの方が向いていると考えられます。UKRは良くも悪くもユーザの恣意性を入れる余地がかなり少ないです。 ↩
経験分布を完全にトレースしてしまい、観測データしか生成しないような状況をここでは指しています。 ↩

深層（でない）生成モデルとしての教師なしカーネル回帰