ガウス過程と機械学習

本記事は書きかけです。最近は主に Zenn で記事を書いており、本記事の内容の拡大板を書籍として公開しているので続きはそちらをご覧ください。

まえがき

この記事は2019年3/11日(日)に大阪大学中之島センターで開講された『機械学習・データ科学スプリングキャンプ初日「ガウス過程と機械学習入門」』の備忘録になります。当日は持橋先生と大羽先生が代わる代わる刺激的な講演をしてくださいました。GPLVMの計算量削減エゲツない。

学生なのでいまTwitterで話題のこちらの書籍¹を無料でいただきました。なんと太っ腹。

注意

本記事ではもともと講義で聴講した内容をなるべくそのまま残そうとしましたが、初稿から記事修正までに時間が経ってしまったこと、解説記事では講義と同じテンポでは説明できないことを踏まえて、私なりに再構成しています。

数式の表記は基本的にガウス過程本に倣います。ベクトル${\bf x}$は縦ベクトルを表します。特に断りのない限り、行列$X$といえばそれは縦ベクトル${\bf x}$を複数個並べたものを表します。同様に$X^T$といえばそれは横ベクトル${\bf x}^T$を複数個並べたものを表します。

また、${\bf y}=f({\bf x})$と書いたときは
$$
\left\{
\begin{eqnarray}
y^{(1)} &= f(x_1^{(1)}, x_2^{(1)}, \ldots, x_n^{(1)}) \\
y^{(2)} &= f(x_1^{(2)}, x_2^{(2)}, \ldots, x_n^{(2)}) \\
\vdots & \\
y^{(m)} &= f(x_1^{(m)}, x_2^{(m)}, \ldots, x_n^{(m)})
\end{eqnarray}
\right. \tag{1}
$$を表しています。

ガウス過程とは（イントロ）

ガウス過程とは一言で言えば「関数をランダムサンプリングする箱」である。

ランダムサンプリングというのは確率変数に対して試行を行い実現値を得る操作のことです。簡単に言えば「サイコロを振って出目を確認する」ということなので、ガウス過程は「出目が関数になっているサイコロを振って目を確認することで関数を取り出す」ようなモデルと言えます。

ガウス過程（Gaussian Process）

関数$x \mapsto y$を生成する確率分布のこと。数学的には関数解析で扱われる対象²で、時間$t$の関数$y=f(t)$として見れば「軌跡」になるが、必ずしも時間の関数でなくてもよい（空間の関数$y=f(x)$でもよい）。

あとで『ガウス過程と機械学習¹』で確認したところ111ページの定義4.1.と定義4.2.の付近にもう少し詳しく書かれていました。これは講義3限目の内容になるので、また別の記事で言及すると思います。

ガウス過程回帰（Gaussian Process Regression）

$(x,y)$に関するデータを集めて既知の集合とし、未知の$x^\ast$ に対して$y^\ast$を予測する、という回帰問題をガウス過程で行う。ガウス過程で回帰を行うことの利点のひとつは、カーネルさえ定義できればどんな特徴量空間$\mathcal{X}$に対してもガウス過程回帰を定義できることである（例：グラフ、化学構造、文字列など）。

講義の後半で明らかにされることですが、ガウス過程はカーネル法で用いられるグラム行列を共分散行列とする多変量ガウス分布から関数を生成します。つまりカーネル法で培われてきた、様々なデータ構造に対するカーネルをそのまま流用できるということです。

出力値$y\in \mathbb{R}$は連続値である³。

パッと見気になるのは「出力値$y\in \mathbb{R}$は」という部分です。当然、多変数$y \in \mathbb{R}^d$のときにどう拡張したらいいかということが気になるわけですが、『Gaussian Process for Machine Learning』にも詳しいことは書いていなかったので探すのはもう少し時間がかかりそうです。

回帰とは

講義では「ガウス過程回帰とは」を説明するまえにまずは「回帰とは」の説明から入り、線形回帰について説明して1限を締めくくりました。ガウス過程やガウス過程回帰は線形回帰をベイズ推論によって拡張したものなので、ここをきっちり抑えておきたいところです。

でもまぁいうて線形回帰は他に無限にいい記事や本が存在するので、この記事では適当にハショります。

『ガウス過程と機械学習』のChapter1に詳しい解説が載っていますし、なんなら私も『線形な手法とカーネル法（1. 回帰分析）』で似たような解説をしているので、そちらもどうぞ。

ただ、ガウス過程の前提になる線形回帰は、以前書いた記事『線形な手法とカーネル法（1. 回帰分析）』に登場したカーネル回帰とは微妙に異なるので、誤解を避けるためもう一度ひと通りの説明はしておきます。

回帰とは

回帰とは入力${\bf x} \in \mathcal{X}$（$\mathcal{X}$は任意の集合）から出力$y \in \mathbb{R}$を予測する問題⁴のことで、入力$x$と出力$y$の関係を
$$
y = f(x)
$$という素直な式で表現したモデルです。回帰問題を解くとは、観測された有限の$N$個のデータ$\mathcal{D}=\{(x^{(1)}, y^{(1)}),(x^{(2)}, y^{(2)}),\ldots, (x^{(N)}, y^{(N)}) \}$を用いてこの未知の関数$f$を推定することを言います。

未知の入力$x^\ast$があっても、それに対応する出力$y^\ast$は求まった関数$f$を用いれば
$$
y^\ast = f(x^\ast)
$$として推定することができます。

単回帰

単回帰は1次元の入力$x \in \mathbb{R}$の値をそのまま特徴量として用いて、直線近似による回帰を行うものです。出力$y$との関係を以下の式

$$
y = a x + b
$$と仮定し、未知の入力$x^\ast$に対応する出力$y^\ast$は
$$
y^\ast = a x^\ast + b
$$と表すことができます。

与えられたデータ$\mathcal{D}=\{(x^{(1)}, y^{(1)}), (x^{(2)}, y^{(2)}), \ldots, (x^{(N)}, y^{(N)})\}$についての入力と出力の関係式を$(1)$式のように連立方程式で表記すると、

$$
\left\{
\begin{eqnarray}
y^{(1)} &= ax^{(1)} + b \\
y^{(2)} &= ax^{(2)} + b \\
\vdots & \\
y^{(N)} &= ax^{(N)} + b
\end{eqnarray}
\right.
$$となります。この連立方程式は行列形式で

$$
\left(
\begin{matrix}
y^{(1)} \\
y^{(2)} \\
\vdots \\
y^{(N)}
\end{matrix}
\right)=\left(
\begin{matrix}
x^{(1)} & 1\\
x^{(2)} & 1\\
\vdots \\
x^{(N)} & 1
\end{matrix}
\right)\left(
\begin{matrix}
a \\
b
\end{matrix}
\right)
$$と書くことができるので、ベクトルと行列の表記に直せば
$$
\mathbf{y} = X\mathbf{w}
$$と書けます。この方程式は$X$が正則であれば逆行列$X^{-1}$を用いて
$$
\mathbf{w} = X^{-1}\mathbf{y}
$$と解くことができます（これができる状況は、与えられた2点を通る直線を求めることに相当します）が、実際は観測されるデータ数は2個よりは多い（$N > 2$）ですから、逆行列は存在せず方程式を正確に満たす$\mathbf{w}=(a,b)^\mathrm{T}$は存在しません。

そこで「観測されたすべての点について、直線との二乗誤差の和をとった値$L$」という「尤もらしい（もっともらしい）指標」を用意してこの値$L$を最小化しよう、という方法で直線を求めるのが「最小二乗法」です。

機械学習においては$L$のような「その値が最小化（または最大化）されたときに、モデルはデータの有様をもっともうまく表現しているであろう」という評価尺度となる関数のことを**損失関数(loss function)**といいます。

二乗誤差を損失関数として用いる場合、損失関数$L$は$\mathcal{D}=(X,\mathbf{y})$と$\mathbf{w}$の関数で、
$$
L(\mathcal{D},\mathbf{w}) = \|\mathbf{y} - X\mathbf{w} \|_2^2
$$と表すことができます。データ$\mathcal{D}$は観測されることで定数となりますので、データから$\mathbf{w}$を求めるときは

$$
L(\mathbf{w}) _ \mathcal{D} = \|\mathbf{y} - X\mathbf{w} \|_2^2 \tag{2}
$$と$\mathbf{w}$のみの関数になります。$L(\mathbf{w}) _ \mathcal{D}$は単に$L(\mathbf{w})$と略記することが多いです。念のために説明しておくと、この$\| \cdot \| _2$はベクトルのユークリッドノルムであり、$\| \cdot \| _ 2 = (\| \cdot \| _ 2)^2$です。

さて、$(2)$式を最小化する方法について詳しくは参考¹や参考²を参照していただきたいですが、条件だけ言うと$L(\mathbf{w})$を$\mathbf{w}$で微分して$0$になる点が最適解です。すなわち
$$
\frac{\partial L(\mathbf{w})}{\partial \mathbf{w}} = -2X^{\mathrm T}(\mathbf{y} -X \mathbf{w}) = 0
$$を満たす$\mathbf{w}$が最適解となります。したがって単回帰モデルの解は、正規方程式
$$
X^\mathrm{T} X w = X^\mathrm{T} y
$$の解として、$X^T X$が逆行列を持つときには
$$
w = (X^\mathrm{T} X)^{-1} X^\mathrm{T} y
$$で表すことができます。

重回帰

重回帰は単回帰で1次元だった入力を多次元$x \in \mathbb{R}^n$に拡張したときのモデルであり、入力$x=(x _ 1, x _ 2, \ldots, x _ n)$の各元$x_i$をそのまま特徴量として、その線型結合で回帰を行うモデルです。
$$
y = \mathbf{w}^\mathrm{T} \mathbf{x} = w_0 \cdot 1 + w_1 x_1 + w_2 x_2 + \cdots + w_n x_n
$$$w_0 \cdot 1$は切片であり、単回帰でいうところの$b$に相当します。

データについての連立方程式は
$$
\left(
\begin{matrix}
y^{(1)} \\
y^{(2)} \\
\vdots \\
y^{(N)}
\end{matrix}
\right)=\left(
\begin{matrix}
1 & x _1^{(1)} & x _2^{(1)} & \cdots & x _n^{(1)}\\
1 & x _1^{(2)} & x _2^{(2)} & \cdots & x _n^{(2)}\\
\vdots \\
1 & x _1^{(N)} & x _2^{(N)} & \cdots & x _n^{(N)}
\end{matrix}
\right)\left(
\begin{matrix}
w _0 \\
w _1 \\
\vdots \\
w _n
\end{matrix}
\right)
$$となるので、ベクトルと行列の表記に直せば単回帰とまったく同じように
$$
\mathbf{y} = X\mathbf{w}
$$と書けますから、その解もまったく同様に
$$
w = (X^\mathrm{T} X)^{-1} X^\mathrm{T} y
$$で表すことができます。

線形回帰

線形回帰は入力$x \in \mathcal{X}$に特徴量写像$\phi \colon \mathcal{X} \mapsto \mathbb{R}^n$を施したあとのベクトル$\phi(x)=(\phi _ 1(x), \phi _ 2 (x), \ldots, \phi _ n (x))^\mathrm{T}$の各元$\phi _ i (x)$を特徴量として、その線型結合で回帰を行うモデルです。$\phi _ i(x)$をこの線形回帰の**基底関数(basis function)**といいます。
$$
y = \mathbf{w}^\mathrm{T} \phi(x) = w _ 1 \phi _ 1(x) + w _ 2 \phi _ 2(x) + \cdots + w_n \phi_n(x) \tag{3}
$$$\mathcal{X}$は任意の集合（文字列やグラフなどでもよい）であり、$\mathcal{X}$から$\mathbb{R}$上への写像が定義できればなんでも構いません。また、参考¹に登場する範囲では常に$\phi _ 0 (x) = 1$として切片に割り当ててあります。
$$
y = w _ 0 \cdot 1 + w _ 1 \phi _ 1(x) + w _ 2 \phi _ 2(x) + \cdots + w_n \phi_n(x) \tag{3}
$$「なんでも」と言ってしまうとあまりに抽象的でイメージがつかなくなるかと思いますので、具体例を挙げながら説明します。入力$x \in \mathbb{R}$として、

$$
\left\{
\begin{eqnarray}
\phi _ 0(x) &= 1 \\
\phi _ 1(x) &= x \\
\phi _ 2(x) &= x^2 \\
\phi _ 3(x) &= x^3 \\
\phi _ 4(x) &= x^4 \\
\phi _ 5(x) &= x^5 \\
\end{eqnarray}
\right.
$$とおいて、$(3)$式に代入すると
$$
y = w _ 0 + w _ 1 x + w _ 2 x^2 + w _ 3 x^3 + w _ 4 x^4 + w _ 5 x^5
$$となりますので、これは$x$に関する5次式になっています。つまり「線形」と言いつつこのモデルで表現できる関数は非線形関数になります。「非線形関数の線形和」になっているわけですね。重回帰は$\phi _i(x) = x _i$であるような、線形回帰の特殊なモデルだったことになります。

この場合もデータに関する連立方程式は

$$
\left(
\begin{matrix}
y^{(1)} \\
y^{(2)} \\
\vdots \\
y^{(N)}
\end{matrix}
\right)=\left(
\begin{matrix}
1 & \phi _1(x)^{(1)} & \phi _2(x)^{(1)} & \cdots & \phi _n(x)^{(1)}\\
1 & \phi _1(x)^{(2)} & \phi _2(x)^{(2)} & \cdots & \phi _n(x)^{(2)}\\
\vdots \\
1 & \phi _1(x)^{(N)} & \phi _2(x)^{(N)} & \cdots & \phi _n(x)^{(N)}
\end{matrix}
\right)\left(
\begin{matrix}
w _0 \\
w _1 \\
w _2 \\
\vdots \\
w _n
\end{matrix}
\right)
$$となるので、ベクトルと行列の表記に直せばやはり
$$
\mathbf{y} = \Phi\mathbf{w}
$$と書けますから、その解もまったく同様に
$$
w = (\Phi^ \mathrm{T} \Phi)^{-1} \Phi^\mathrm{T} y
$$で表すことができます。

まだよく調べてないですがこのときの$\Phi$を計画行列(design matrix)というそうです。

一般化線形モデル

統計の分野で出てくる一般化線形モデル（GLM: Generalized Linear Model）は
$$
y = f(w^T \phi(x))
$$のことを指すので注意とのこと。このときの$f$をリンク関数というそうです。

リッジ回帰

さて、今までのモデルでは$\Phi ^ \mathrm{T} \Phi$に逆行列が存在していることを仮定していましたが、必ずしも存在するとは限りません。

もっとも起こりうるパターンは、データ$\mathcal{D}$の中で同じ値を持つデータが2回以上観測されているケースです。この場合$\Phi ^ \mathrm{T} \Phi$にはランク落ちが発生しており逆行列が存在しません。また、まったく同じでなくても「非常に近い値のデータ」が観測されていると逆行列は観測ノイズの影響をシビアに受けることになります¹。

そこで$w$に対してペナルティを課して
$$
\|y - \Phi w \| _2^2 + \alpha \|w\| _2^2
$$という目的関数を最小化します。このとき$w$は
$$
w = (\Phi^\mathrm{T} \Phi + \alpha I)^{-1} \Phi^ \mathrm{T} \mathbf{y}
$$と求まり、$\Phi^\mathrm{T} \Phi$がランク落ちしていても$\alpha I$が加算されてフルランクになるので常に逆行列が計算できます。

おまけ：カーネル回帰

カーネル回帰は線形回帰の一種ですが「どんな特徴量写像を定義していいかわからないので、観測されたデータ$x'$を特徴量写像$\phi(x)$に変換するカーネル関数$k(x,x')$を定義して、観測されたデータから特徴量写像を作ってしまおう」というモデルです。

カーネル関数$k(x,x')$は正定値性さえ満たせばなんでもよいのですが、たとえばガウスカーネルと呼ばれるものは以下の式で定義されます。
$$
k(x,x')= exp(-\beta \|x - x'\| _2^2)
$$このとき、観測されたデータ$\mathcal{D}=\{(x^{(1)}, y^{(1)}), (x^{(2)}, y^{(2)}), \ldots, (x^{(N)}, y^{(N)})\}$の$x^{(i)}$に対応した基底関数は
$$
\phi _i(x) = k(x, x^{(i)}) = exp(-\beta \|x - x^{(i)}\| _2^2)
$$となります。

カーネル回帰はこうして定義した特徴量写像を用いて線形回帰を行うモデルです。
$$
\begin{align}
y &= \mathbf{w}^\mathrm{T} \phi(x) = w _ 1 \phi _ 1(x) + w _ 2 \phi _ 2(x) + \cdots + w_n \phi_n(x) \\
&= w _ 1 k(x, x^{(1)}) + w _ 2 k(x, x^{(2)}) + \cdots + w_N k(x, x^{(N)})
\end{align}
$$観測されたデータについての連立方程式は
$$
\left(
\begin{matrix}
y^{(1)} \\
y^{(2)} \\
\vdots \\
y^{(N)}
\end{matrix}
\right)=\left(
\begin{matrix}
k(x^{(1)}, x^{(1)}) & k(x^{(1)}, x^{(2)}) & \cdots & k(x^{(1)}, x^{(N)})\\
k(x^{(2)}, x^{(1)}) & k(x^{(2)}, x^{(2)}) & \cdots & k(x^{(2)}, x^{(N)})\\
\vdots \\
k(x^{(N)}, x^{(1)}) & k(x^{(N)}, x^{(2)}) & \cdots & k(x^{(N)}, x^{(N)})
\end{matrix}
\right)\left(
\begin{matrix}
w _1 \\
w _2 \\
\vdots \\
w _N
\end{matrix}
\right)
$$と書けるので、ベクトルと行列の表記に直せば
$$
\mathbf{y} = K\mathbf{w}
$$となります。ここで注目していただきたいのは、行列$K$を計算するために必要なのはカーネル関数$k(x,x')$のみであり、特徴量写像$\phi(x)$を明示的に求める必要はないということです。これをカーネルトリックといいます。

カーネル回帰は最小二乗法で解くと過学習を起こすので、リッジ回帰による正則化を行い、
$$
\|\mathbf{y} - K \mathbf{w} \| _2^2 + \alpha \mathbf{w}^\mathrm{T}K\mathbf{w}
$$という損失関数を最小化します。

線形回帰と確率モデル

1限目の最後は、線形回帰と確率モデルの関連についての説明でした。

線形回帰モデルとガウス分布

線形回帰のモデルにノイズの混入を仮定し、
$$
y = f(x) + \varepsilon = w^T \phi(x)+\varepsilon
$$とします。ただし$\varepsilon$はホワイトノイズで平均$0$、分散$\sigma^2$の正規分布にしたがう確率変数であると仮定すると、
$$
\begin{align}
y - w^T \phi(x) \sim \mathcal{N}(0, \sigma^2) \\
\therefore y \sim \mathcal{N}(w^T \phi(x), \sigma^2)
\end{align}
$$により$y$も確率変数となります。

注：この時点で既に$x,y$の関係を確率的生成モデルで表現した、と考えることもできる。

$X = (x_1, x_2, \ldots, x_N)$が与えられた下での$y = (y_1, y_2, \ldots, y_N)$の確率（$y$の尤度）は、観測が独立だとすれば
$$
p(y|X) = \prod_{n=1}^N p(y_n|x_n) = \prod_{n=1}^N \mathcal{N}(w^T \phi(x)|\sigma^2)
$$となるから、両辺で対数を取れば
$$
\begin{align}
\log p(y|X)&=\sum_{n=1}^N \log \mathcal{N}(y_n | w^T \phi(x_n), \sigma^2) \\
&= - \|y- \Phi w\|^2 + const.
\end{align}
$$となり、最小二乗法の損失関数と同じ項が現れることがわかる。つまりこの対数尤度$\log p(y|X)$を最大化する問題（最尤推定）が最小二乗法と等価であることが見て取れるわけである。

線形回帰モデルとリッジ回帰

先ほどのモデルでは$w$についてなんら事前の情報を仮定しなかった。これだと$w$の値が$10^{100000}$のような「非常識な値」を取ったところで文句は言えない、ということになるので適当に$10^6$とかコンピュータで扱える程度の「常識的な値」に収まってほしいという事前情報をモデルに加えることにする。

$w$の各要素$w_i$がそれぞれガウス分布$\mathcal{N}(0,\sigma^2)$にしたがうとすれば、$y$と$w$の同時確率の対数は
$$
\begin{align}
\log p(y,w |X) &= \log p(y|w,X)p(w|X) \\
&= -(\|y-\Phi w\|^2 + \alpha \|w\|^2)
\end{align}
$$と表せるので、リッジ回帰と等価であることがわかる。

リッジ回帰の式を導出するにはMAP推定（事後確率の最大化）と絡めることが多いのだが、同時確率の最大化でも求めることができるのは驚きだ。どう違うのか（または本質的には同じことを言っているのか）あとで検算しておきたい。

1. 『ガウス過程と機械学習』持橋大地, 大羽成征, 講談社 ↩ ↩² ↩³ ↩⁴ ↩⁵

2. 『線形な手法とカーネル法（1. 回帰分析）』 ↩ ↩²

3. この主張は自明ではないです。「出力値が実数値を取る」ことと「出力値が連続になる」ことは厳密には等価ではなく、ここまでの議論では「出力値$y$は連続確率分布から生成された実数である」ことしか言っていないので「出力値が連続になる」ことは自明ではありません。この点に関しての詳細は講義のSlackで寺田さんという方が「そのうち, 関数データ解析のチュートリアルなど担当する機会があれば資料upする可能性あります．寺田 吉壱でググって頂くとそのうち資料など見れるかもしれません．」と書いていらっしゃったので、期待に胸を膨らませつつ気長に待とうと思います。この辺りの厳密な話はまだ調べている最中ですので、もしわかったら3限目の記事に書こうと思います。 ↩

4. 教師あり学習には大きく分けて回帰問題と分類問題の2種類があります。ガウス過程はどちらの問題にも応用できるモデルのようですが、講義では分類問題は扱いませんでした。 ↩