数百次元の手書き数字画像、数十関節からなる人間のモーションキャプチャデータ、数千遺伝子の発現量——こうした「高次元だが、実は少数の自由度で動いている」データは現実のあちこちに転がっています。手書きの「3」の画像は $8\times8=64$ 次元のベクトルですが、字の太さや傾き、丸みといった数個のパラメータが本質的に支配しているはずです。この「本質的な少数の自由度」を見つけ出すのが次元削減です。
線形の次元削減なら主成分分析（PCA）が定番です。しかしPCAは線形写像しか表現できないので、データが曲がった多様体（manifold）の上に乗っていると、その曲がりをほどけずに潰してしまいます。では、非線形の次元削減を、しかも確率モデルとしてきちんと定式化するにはどうすればよいのか。その美しい答えの一つが、Neil Lawrence が2005年に提案したガウス過程潜在変数モデル（Gaussian Process Latent Variable Model, GP-LVM） です。
GP-LVMのアイデアは驚くほどシンプルです。確率的PCA（PPCA）が「潜在変数 → 観測」を線形写像で結ぶのに対し、GP-LVMはその写像をガウス過程に置き換えます。たったこれだけで、線形だったPCAが非線形多様体を扱える強力な手法に化けるのです。この手法は、高次元データの可視化、人間の動作データのモデル化、欠損値の補完、新しいデータの生成など、幅広い場面で使われています。
本記事の内容
- PCAの限界と「非線形次元削減を確率モデルで書きたい」というモチベーション
- 確率的PCA（PPCA）の復習 — 潜在変数モデルとしての定式化
- PPCAの線形写像をガウス過程に置き換えるGP-LVMの定式化
- 対数周辺尤度と、潜在変数 $\boldsymbol{X}$ に関する勾配の導出
- GP-LVMとPCAの厳密な関係（線形カーネル ⇔ PPCA）
- ベイズGP-LVM（Titsias & Lawrence 2010）による潜在変数の積分消去とARD
- Pythonによるスクラッチ実装と、digitsデータの2次元可視化
- 動作データ・欠損補完・生成への応用
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
なぜ非線形次元削減が必要なのか
まず、PCAだけでは何が足りないのかを具体的に見ておきましょう。PCAは、データの分散が最大になる方向（主成分）に沿ってデータを射影します。これは「データを最もよく説明する直線的な部分空間を見つける」操作です。データが平たい板の上に分布しているなら、これで十分です。
ところが現実のデータは、しばしば曲がった面の上に乗っています。下の図は、3次元空間の中で「スイスロール」状に巻いた2次元の面の上にデータが分布している例です。本当はこの面を平らに伸ばせば（巻きをほどけば）2次元で完璧に表現できるのですが、PCAでこれを2次元に射影するとどうなるでしょうか。
左が3次元の元データ（色は多様体上の位置を表します）、右がPCAで2次元に射影した結果です。PCAは巻きをほどけず、本来は遠く離れているはずの場所（色の違う点）が右図では重なってしまっています。直線的な射影では、丸めた紙を上から押しつぶすように、構造が壊れてしまうのです。
私たちが欲しいのは、この「巻きをほどく」非線形な操作です。しかも単なる座標変換ではなく、確率モデルとして定式化したい。確率モデルにすれば、不確実性の評価、欠損値の補完、新しいサンプルの生成といった芸当が自然にできるようになるからです。この目標を達成するために、まずは線形の確率的次元削減である確率的PCAを復習し、それを足がかりにGP-LVMへと拡張していきます。
確率的PCA（PPCA）の復習
PCAを確率モデルとして書き直す
通常のPCAは分散最大化や再構成誤差最小化として説明されますが、Tipping & Bishop (1999) は、PCAを潜在変数を持つ確率モデルとして書けることを示しました。これが確率的PCA（Probabilistic PCA, PPCA）です。発想はこうです。観測できる高次元データ $\boldsymbol{y}_n \in \mathbb{R}^D$ の裏には、観測できない低次元の潜在変数 $\boldsymbol{x}_n \in \mathbb{R}^Q$（$Q \ll D$）が隠れていて、それが線形写像で高次元に持ち上げられ、ノイズが乗って観測される、と考えます。
図のように、各データ点 $n$ について、潜在変数 $\boldsymbol{x}_n$ から観測 $\boldsymbol{y}_n$ が生成されます。生成過程を数式で書くと次のようになります。
$$ \boldsymbol{x}_n \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{I}_Q), \qquad \boldsymbol{y}_n = \boldsymbol{W}\boldsymbol{x}_n + \boldsymbol{\mu} + \boldsymbol{\epsilon}_n, \qquad \boldsymbol{\epsilon}_n \sim \mathcal{N}(\boldsymbol{0}, \sigma^2 \boldsymbol{I}_D) $$
ここで $\boldsymbol{W} \in \mathbb{R}^{D \times Q}$ は潜在から観測へのローディング行列、$\boldsymbol{\mu}$ はデータ平均、$\boldsymbol{\epsilon}_n$ は等方的なガウスノイズです。潜在変数の事前分布を標準正規 $\mathcal{N}(\boldsymbol{0}, \boldsymbol{I})$ に固定するのがミソで、これによりモデルがすっきりします。
潜在変数を周辺化する
PPCAでは、潜在変数 $\boldsymbol{x}_n$ を積分で消去（周辺化）して、観測 $\boldsymbol{y}_n$ の分布を求めます。$\boldsymbol{y}_n$ はガウス分布する $\boldsymbol{x}_n$ の線形変換にガウスノイズを足したものなので、$\boldsymbol{y}_n$ 自身もガウス分布します。平均と共分散を計算しましょう。平均は
$$ \mathbb{E}[\boldsymbol{y}_n] = \boldsymbol{W}\,\mathbb{E}[\boldsymbol{x}_n] + \boldsymbol{\mu} = \boldsymbol{\mu} $$
です（$\mathbb{E}[\boldsymbol{x}_n]=\boldsymbol{0}$ を使いました）。共分散は、$\boldsymbol{x}_n$ と $\boldsymbol{\epsilon}_n$ が独立であることを使って
$$ \mathrm{Cov}[\boldsymbol{y}_n] = \boldsymbol{W}\,\mathbb{E}[\boldsymbol{x}_n\boldsymbol{x}_n^\top]\,\boldsymbol{W}^\top + \mathbb{E}[\boldsymbol{\epsilon}_n\boldsymbol{\epsilon}_n^\top] = \boldsymbol{W}\boldsymbol{W}^\top + \sigma^2 \boldsymbol{I}_D $$
と求まります。途中、$\mathbb{E}[\boldsymbol{x}_n\boldsymbol{x}_n^\top]=\boldsymbol{I}_Q$（事前分布が標準正規）を代入しました。まとめると、観測の周辺分布は
$$ \boldsymbol{y}_n \sim \mathcal{N}\big(\boldsymbol{\mu},\; \boldsymbol{C}\big), \qquad \boldsymbol{C} = \boldsymbol{W}\boldsymbol{W}^\top + \sigma^2 \boldsymbol{I}_D $$
となります。データ全体の対数尤度 $\log p(\boldsymbol{Y} \mid \boldsymbol{W}, \sigma^2) = \sum_n \log \mathcal{N}(\boldsymbol{y}_n \mid \boldsymbol{\mu}, \boldsymbol{C})$ を $\boldsymbol{W}$ について最大化すると、その最尤解はデータの共分散行列の上位 $Q$ 個の固有ベクトル、つまりPCAの主成分に一致することが示せます。PPCAは、PCAをきれいに確率モデル化したものなのです。
ここで重要な視点の切り替えをします。PPCAは「観測次元 $\boldsymbol{y}_n$ ごと」ではなく、見方を変えて「観測次元 $d$ ごと」に分解することもできます。次元 $d$ の観測値を $N$ 個並べたベクトルを $\boldsymbol{y}_{:,d} \in \mathbb{R}^N$ とすると、これは潜在 $\boldsymbol{X}$ の各行を入力とした $N$ 個の値の集まりです。この「次元ごとに見る」視点こそが、GP-LVMへの扉を開きます。
PPCAは線形写像 $\boldsymbol{W}$ に縛られているので、曲がった多様体は扱えません。では、この $\boldsymbol{W}\boldsymbol{x}$ という線形写像を、もっと柔軟な非線形写像に取り替えたらどうなるでしょうか。
GP-LVMの定式化
線形写像をガウス過程に置き換える
GP-LVMの核心的なアイデアは、PPCAの線形写像 $\boldsymbol{y} = \boldsymbol{W}\boldsymbol{x} + \boldsymbol{\epsilon}$ における写像を、ガウス過程（GP） に置き換えることです。観測の各次元 $d$ が、潜在変数 $\boldsymbol{x}$ を入力とする独立な関数 $f_d$ の出力だと考えます。
つまり、観測値を
$$ y_{nd} = f_d(\boldsymbol{x}_n) + \epsilon_{nd}, \qquad f_d \sim \mathcal{GP}\big(0, k(\cdot, \cdot)\big), \qquad \epsilon_{nd} \sim \mathcal{N}(0, \beta^{-1}) $$
と書きます。各 $f_d$ は平均ゼロ、共分散関数（カーネル）$k$ のガウス過程に従う未知の非線形関数です。$\beta^{-1}$ はノイズ分散です。カーネルとしてはRBF（ガウス）カーネルがよく使われます。
$$ k(\boldsymbol{x}_i, \boldsymbol{x}_j) = \alpha \exp\!\left(-\gamma \,\|\boldsymbol{x}_i – \boldsymbol{x}_j\|^2\right) $$
ここで $\alpha$ は信号分散、$\gamma$ は逆長さスケールです。下の図は、1次元の潜在変数 $x$ をGPで2次元の観測空間に写像した例です。潜在空間の滑らかな順序が、観測空間の滑らかな非線形曲線に対応していることがわかります。
左は1次元の潜在座標 $x$（色で表示）、右はそれをGP写像 $f=(f_1, f_2)$ で2次元の観測空間に持ち上げた結果です。色の連続性が観測空間でも保たれ、線形写像では決して作れない曲がった多様体が生成されています。GPを使うことで、PPCAの「平らな板」が「自在に曲がる面」になったわけです。
次元ごとに周辺化する
PPCAでは潜在変数 $\boldsymbol{x}$ を周辺化しました。GP-LVMでは逆に、写像 $f_d$ を周辺化します。ガウス過程の定義から、関数 $f_d$ を有限個の入力点 $\{\boldsymbol{x}_n\}_{n=1}^N$ で評価した値 $\boldsymbol{f}_d = (f_d(\boldsymbol{x}_1), \dots, f_d(\boldsymbol{x}_N))^\top$ は、多変量ガウス分布に従います。
$$ \boldsymbol{f}_d \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{K}) $$
ここで $\boldsymbol{K}$ は $K_{ij} = k(\boldsymbol{x}_i, \boldsymbol{x}_j)$ で定まる $N \times N$ のグラム行列（カーネル行列）です。観測 $\boldsymbol{y}_{:,d} = \boldsymbol{f}_d + \boldsymbol{\epsilon}_d$ はこれにノイズを足したものなので、
$$ \boldsymbol{y}_{:,d} \sim \mathcal{N}\big(\boldsymbol{0},\; \boldsymbol{K} + \beta^{-1}\boldsymbol{I}_N\big) $$
となります。ここからは記号を簡単にするため、ノイズを含めたカーネル行列を $\boldsymbol{K}$ と書き直します（$\boldsymbol{K} \leftarrow \boldsymbol{K} + \beta^{-1}\boldsymbol{I}$）。各観測次元 $d$ は独立で、しかも同じ潜在 $\boldsymbol{X}$ を共有する同じ共分散 $\boldsymbol{K}$ を持つ、というのが重要なポイントです。
ここで、PPCAとGP-LVMの双対性が見えてきます。PPCAは「データ点 $\boldsymbol{x}_n$ をガウス、写像 $\boldsymbol{W}$ をパラメータ」として扱いました。GP-LVMは逆に「写像 $f_d$ をガウス（過程）、潜在 $\boldsymbol{x}_n$ をパラメータ」として扱います。周辺化する対象が入れ替わっているのです。では、この定式化から尤度を書き下し、潜在 $\boldsymbol{X}$ をどう求めるかを見ていきましょう。
尤度と潜在変数Xの最適化
対数周辺尤度
各観測次元 $\boldsymbol{y}_{:,d}$ が独立にガウス分布するので、全データ $\boldsymbol{Y} \in \mathbb{R}^{N \times D}$（各列が $\boldsymbol{y}_{:,d}$）の対数周辺尤度は、$D$ 個のガウス尤度の和になります。
$$ \log p(\boldsymbol{Y} \mid \boldsymbol{X}) = \sum_{d=1}^{D} \log \mathcal{N}(\boldsymbol{y}_{:,d} \mid \boldsymbol{0}, \boldsymbol{K}) $$
多変量ガウスの対数密度 $\log \mathcal{N}(\boldsymbol{y} \mid \boldsymbol{0}, \boldsymbol{K}) = -\frac{N}{2}\log(2\pi) – \frac{1}{2}\log|\boldsymbol{K}| – \frac{1}{2}\boldsymbol{y}^\top \boldsymbol{K}^{-1}\boldsymbol{y}$ を代入します。すべての次元で共分散 $\boldsymbol{K}$ が共通なので、$\log(2\pi)$ 項と $\log|\boldsymbol{K}|$ 項は $D$ 倍され、二次形式の項は和としてまとまります。
$$ \log p(\boldsymbol{Y} \mid \boldsymbol{X}) = -\frac{DN}{2}\log(2\pi) – \frac{D}{2}\log|\boldsymbol{K}| – \frac{1}{2}\sum_{d=1}^{D} \boldsymbol{y}_{:,d}^\top \boldsymbol{K}^{-1} \boldsymbol{y}_{:,d} $$
最後の和はトレースを使ってきれいにまとめられます。$\sum_d \boldsymbol{y}_{:,d}^\top \boldsymbol{K}^{-1} \boldsymbol{y}_{:,d} = \mathrm{tr}\!\left(\boldsymbol{K}^{-1} \sum_d \boldsymbol{y}_{:,d}\boldsymbol{y}_{:,d}^\top\right) = \mathrm{tr}(\boldsymbol{K}^{-1} \boldsymbol{Y}\boldsymbol{Y}^\top)$ と変形できる（行列の積の和の定義そのもの）からです。したがって
$$ \log p(\boldsymbol{Y} \mid \boldsymbol{X}) = -\frac{DN}{2}\log(2\pi) – \frac{D}{2}\log|\boldsymbol{K}| – \frac{1}{2}\mathrm{tr}\!\left(\boldsymbol{K}^{-1}\boldsymbol{Y}\boldsymbol{Y}^\top\right) $$
が得られます。これがGP-LVMの目的関数です。ここで肝心なのは、$\boldsymbol{K}$ が潜在変数 $\boldsymbol{X}$ の関数だという点です。$\boldsymbol{X}$ を動かすと $\boldsymbol{K}$ が変わり、尤度が変わります。GP-LVMの学習とは、この尤度を最大にする潜在配置 $\boldsymbol{X}$（とカーネルハイパラ $\alpha, \gamma, \beta$）を求めることに他なりません。
潜在変数Xに関する勾配
$\boldsymbol{X}$ を勾配法で最適化するために、対数尤度 $\mathcal{L} = \log p(\boldsymbol{Y} \mid \boldsymbol{X})$ の勾配を求めます。$\boldsymbol{X}$ への依存はすべてカーネル行列 $\boldsymbol{K}$ を通じてなので、連鎖律を使います。まず $\mathcal{L}$ を $\boldsymbol{K}$ で微分しましょう。行列微分の公式 $\frac{\partial \log|\boldsymbol{K}|}{\partial \boldsymbol{K}} = \boldsymbol{K}^{-1}$ と $\frac{\partial\, \mathrm{tr}(\boldsymbol{K}^{-1}\boldsymbol{A})}{\partial \boldsymbol{K}} = -\boldsymbol{K}^{-1}\boldsymbol{A}\boldsymbol{K}^{-1}$ を用います。$\boldsymbol{A} = \boldsymbol{Y}\boldsymbol{Y}^\top$ とおくと、
$$ \frac{\partial \mathcal{L}}{\partial \boldsymbol{K}} = -\frac{D}{2}\boldsymbol{K}^{-1} + \frac{1}{2}\boldsymbol{K}^{-1}\boldsymbol{Y}\boldsymbol{Y}^\top\boldsymbol{K}^{-1} $$
となります。第1項が $\log|\boldsymbol{K}|$ 由来、第2項が二次形式由来です。次に、カーネル行列の各要素 $K_{ij}$ を潜在座標 $\boldsymbol{x}_n$ の各成分で微分し、連鎖律で組み合わせます。
$$ \frac{\partial \mathcal{L}}{\partial x_{nq}} = \sum_{i, j}\frac{\partial \mathcal{L}}{\partial K_{ij}}\frac{\partial K_{ij}}{\partial x_{nq}} $$
RBFカーネル $K_{ij} = \alpha \exp(-\gamma \|\boldsymbol{x}_i – \boldsymbol{x}_j\|^2)$ の場合、$K_{ij}$ を $x_{nq}$ で微分すると、$\boldsymbol{x}_i$ か $\boldsymbol{x}_j$ のどちらかが $\boldsymbol{x}_n$ のときだけ非ゼロになり、
$$ \frac{\partial K_{ij}}{\partial x_{nq}} = -2\gamma\, K_{ij}\,(x_{iq} – x_{jq})\,(\delta_{in} – \delta_{jn}) $$
という形になります（$\delta$ はクロネッカーのデルタ）。これらを代入すれば、潜在変数の全勾配が解析的に計算できます。あとはこの勾配を使って、L-BFGSなどの準ニュートン法で $\boldsymbol{X}$ を更新します。
最適化の実際
GP-LVMの最適化は非凸なので、初期値が重要です。標準的にはPCAの結果で $\boldsymbol{X}$ を初期化します。線形PCAが大まかな配置を与え、GP-LVMがそれを非線形に洗練する、という流れです。下の図は、合成データに対してL-BFGSで尤度を最適化したときの収束の様子です。
横軸が最適化のイテレーション、縦軸が対数周辺尤度です。イテレーションを重ねるごとに尤度が単調に増加し、やがて頭打ちになって収束しています。潜在配置 $\boldsymbol{X}$ とカーネルハイパラを同時に動かしながら、データを最もよく説明する低次元表現が自動的に見つかっていく様子が読み取れます。非凸ゆえ初期値依存はありますが、PCA初期化のおかげで安定して良い解にたどり着けます。
尤度の最適化で潜在 $\boldsymbol{X}$ が求まることがわかりました。ここで自然な疑問が湧きます。「カーネルを線形にしたら、GP-LVMはPCAに戻るのだろうか？」 実はその通りで、両者には厳密な対応関係があります。
PCAとの関係
線形カーネルのGP-LVMはPPCAになる
GP-LVMでカーネルを線形カーネル $k(\boldsymbol{x}_i, \boldsymbol{x}_j) = \boldsymbol{x}_i^\top \boldsymbol{x}_j$ に選ぶと、何が起きるか見てみましょう。このときグラム行列は $\boldsymbol{K} = \boldsymbol{X}\boldsymbol{X}^\top$（ノイズを含めると $\boldsymbol{X}\boldsymbol{X}^\top + \beta^{-1}\boldsymbol{I}$）になります。対数尤度は
$$ \log p(\boldsymbol{Y} \mid \boldsymbol{X}) = -\frac{DN}{2}\log(2\pi) – \frac{D}{2}\log|\boldsymbol{X}\boldsymbol{X}^\top + \beta^{-1}\boldsymbol{I}| – \frac{1}{2}\mathrm{tr}\!\left((\boldsymbol{X}\boldsymbol{X}^\top + \beta^{-1}\boldsymbol{I})^{-1}\boldsymbol{Y}\boldsymbol{Y}^\top\right) $$
となります。Lawrence (2005) は、この尤度を $\boldsymbol{X}$ について最大化した解が、PPCA（すなわちPCA）の解と一致することを示しました。直感的には、線形カーネル $\boldsymbol{X}\boldsymbol{X}^\top$ は「潜在座標の内積」であり、これは $\boldsymbol{y} = \boldsymbol{W}\boldsymbol{x}$ という線形写像と同じ構造を持つからです。
ここに面白い双対性があります。PPCAは観測ごとの共分散 $\boldsymbol{W}\boldsymbol{W}^\top + \sigma^2\boldsymbol{I}$（$D\times D$ 行列）を扱い、$\boldsymbol{W}$ を最尤推定します。一方、線形GP-LVMはデータ点間の共分散 $\boldsymbol{X}\boldsymbol{X}^\top + \beta^{-1}\boldsymbol{I}$（$N \times N$ 行列）を扱い、$\boldsymbol{X}$ を最尤推定します。同じPCA解を、「特徴方向の共分散」から見るか「サンプル間の共分散」から見るかの違いなのです。これは、PCAの主成分（特徴方向）と主座標（サンプルの埋め込み）の双対関係そのものです。
非線形カーネルにすると何が変わるか
この対応関係が美しいのは、「線形カーネル ⇔ PCA」を出発点として、カーネルを非線形に変えるだけで非線形次元削減に拡張できる点です。RBFカーネルにすれば、グラム行列 $\boldsymbol{X}\boldsymbol{X}^\top$ が非線形な $\alpha\exp(-\gamma\|\boldsymbol{x}_i-\boldsymbol{x}_j\|^2)$ に置き換わり、PCAでは表現できなかった曲がった多様体を捉えられるようになります。カーネルトリックと同じ精神です。
下の図は、手書き数字データ（0〜4）に対して、線形PCAとGP-LVMで2次元埋め込みを行い比較したものです。
左の線形PCAでは異なる数字のクラスタが重なり気味なのに対し、右のGP-LVMでは各数字のクラスタがより明瞭に分離しています。非線形カーネルが、数字の見た目の複雑な変動を低次元の潜在空間できれいにほどいているのです。線形と非線形の差が、可視化の質に直結することが一目でわかります。
GP-LVMの基本形では、潜在 $\boldsymbol{X}$ を点推定（最尤・MAP）しました。しかしベイズの立場からすれば、$\boldsymbol{X}$ もまた不確実な変数であり、本来は積分消去すべきものです。それを可能にしたのが、次に見るベイズGP-LVMです。
ベイズGP-LVM
潜在変数を積分消去する
基本のGP-LVMには弱点があります。潜在変数 $\boldsymbol{X}$ を点推定するため、潜在次元数 $Q$ を事前に決め打ちする必要があり、過学習のリスクもあります。Titsias & Lawrence (2010) は、$\boldsymbol{X}$ にも事前分布 $p(\boldsymbol{X}) = \prod_n \mathcal{N}(\boldsymbol{x}_n \mid \boldsymbol{0}, \boldsymbol{I})$ を置き、$\boldsymbol{X}$ を変分推論で積分消去するベイズGP-LVMを提案しました。
理想的には、$\boldsymbol{X}$ を周辺化した周辺尤度
$$ p(\boldsymbol{Y}) = \int p(\boldsymbol{Y} \mid \boldsymbol{X})\, p(\boldsymbol{X})\, d\boldsymbol{X} $$
を最大化したいところです。しかし $\boldsymbol{X}$ はカーネル行列 $\boldsymbol{K}$ の中に非線形に入り込んでいるため、この積分は解析的に解けません。ここで変分推論の出番です。$\boldsymbol{X}$ の事後分布を扱いやすいガウス分布 $q(\boldsymbol{X}) = \prod_n \mathcal{N}(\boldsymbol{x}_n \mid \boldsymbol{m}_n, \boldsymbol{S}_n)$ で近似し、対数周辺尤度の下界（ELBO, Evidence Lower BOund）を最大化します。
$$ \log p(\boldsymbol{Y}) \geq \mathbb{E}_{q(\boldsymbol{X})}\big[\log p(\boldsymbol{Y} \mid \boldsymbol{X})\big] – \mathrm{KL}\big(q(\boldsymbol{X}) \,\|\, p(\boldsymbol{X})\big) $$
右辺第1項はモデルがデータをどれだけよく説明するか、第2項は近似事後分布が事前分布からどれだけ離れているかのペナルティです。ただし第1項の期待値も、$\boldsymbol{X}$ がカーネルの中に入るため直接は計算できません。Titsias & Lawrence は、補助変数（inducing points） を導入したスパース近似と組み合わせることで、この期待値を解析的に評価できるようにしました（詳細はスパースガウス過程の議論に通じます）。
ARDによる自動次元選択
ベイズGP-LVMのもう一つの利点が、自動関連度決定（ARD, Automatic Relevance Determination） です。RBFカーネルを、潜在次元ごとに独立な逆長さスケール $w_q$ を持つARDカーネルに拡張します。
$$ k(\boldsymbol{x}_i, \boldsymbol{x}_j) = \alpha \exp\!\left(-\frac{1}{2}\sum_{q=1}^{Q} w_q\, (x_{iq} – x_{jq})^2\right) $$
各次元 $q$ の重み $w_q$ を学習すると、その次元がデータの説明にどれだけ寄与しているかが自動的にわかります。不要な次元の $w_q$ は0に近づき、その次元は実質的に無視されます。つまり、$Q$ を大きめに設定しておけば、本当に必要な潜在次元数をモデルが自動で決めてくれるのです。
この図は、本来2次元の構造を持つデータに対し、潜在次元を8まで許したときのARD重みです。第1・第2次元の重みだけが大きく残り（赤）、残りの次元の重みはほぼ0に落ちています（灰）。閾値（破線）を超える次元だけが有効と判断され、データの本質的な次元数「2」が自動的に発見されています。手で $Q$ を選ぶ必要がなくなるのは、実用上とても大きな利点です。
ベイズGP-LVMによって、次元数の自動決定と過学習の抑制が手に入りました。ここまで理論を積み上げてきたので、実際にスクラッチ実装して動かし、これらの主張を自分の目で確かめましょう。
Pythonでの実装
GP-LVMをスクラッチ実装する
まず、対数尤度を計算し、L-BFGSで潜在変数 $\boldsymbol{X}$ とカーネルハイパラを同時に最適化するGP-LVMをゼロから実装します。前節で導いた対数尤度 $\mathcal{L} = -\frac{D}{2}\log|\boldsymbol{K}| – \frac{1}{2}\mathrm{tr}(\boldsymbol{K}^{-1}\boldsymbol{Y}\boldsymbol{Y}^\top)$（定数項は省略）をそのままコードに落とします。
import numpy as np
from scipy.optimize import minimize
from scipy.spatial.distance import cdist
from sklearn.decomposition import PCA
def rbf_K(X, alpha, gamma, beta_inv=None):
"""RBFカーネルのグラム行列 K_ij = alpha*exp(-gamma*||xi-xj||^2)"""
d2 = cdist(X, X, "sqeuclidean")
K = alpha * np.exp(-gamma * d2)
if beta_inv is not None:
K = K + beta_inv * np.eye(len(X)) # 観測ノイズを対角に加える
return K
def neg_log_lik(params, Y, N, Q):
"""GP-LVM の負の対数尤度（X とハイパラを一緒に最適化）
params = [X.flatten(), log_alpha, log_gamma, log_beta_inv]"""
X = params[:N * Q].reshape(N, Q)
log_alpha, log_gamma, log_beta_inv = params[N * Q:]
alpha, gamma, beta_inv = np.exp(log_alpha), np.exp(log_gamma), np.exp(log_beta_inv)
D = Y.shape[1]
K = rbf_K(X, alpha, gamma, beta_inv)
# コレスキー分解で logdet と逆行列を安定に計算
L = np.linalg.cholesky(K + 1e-6 * np.eye(N))
Kinv = np.linalg.solve(L.T, np.linalg.solve(L, np.eye(N)))
logdetK = 2.0 * np.sum(np.log(np.diag(L)))
# L = -D/2 logdet - 1/2 tr(Kinv Y Y^T) の符号を反転
return 0.5 * D * logdetK + 0.5 * np.trace(Kinv @ Y @ Y.T)
ハイパラは正値制約があるので、log_alpha のように対数で持って最適化し、np.exp で正値に戻しています。カーネル行列の逆行列と行列式は、数値的に安定なコレスキー分解 np.linalg.cholesky 経由で計算しているのがポイントです。1e-6 * np.eye(N) は数値安定化のためのジッターです。
続いて、PCAで初期化してL-BFGSで最適化する学習関数を書きます。
def fit_gplvm(Y, Q=2, n_iter=200, seed=0):
"""PCA初期化 + L-BFGS で GP-LVM を学習"""
N, D = Y.shape
Yc = Y - Y.mean(0) # 中心化
X0 = PCA(n_components=Q, random_state=seed).fit_transform(Yc)
X0 = X0 / X0.std(0) # 潜在のスケールを正規化
init = np.concatenate([X0.ravel(),
[np.log(1.0), np.log(0.5), np.log(0.1)]])
history = []
def cb(p): # 収束ログ（対数尤度）を記録
history.append(-neg_log_lik(p, Yc, N, Q))
cb(init)
res = minimize(neg_log_lik, init, args=(Yc, N, Q),
method="L-BFGS-B", callback=cb,
options={"maxiter": n_iter})
X = res.x[:N * Q].reshape(N, Q)
return X, np.array(history), res
潜在配置をPCAで初期化し、L-BFGSの勾配は数値微分（SciPyの既定）に任せています。導出した解析勾配を渡せばさらに高速ですが、ここでは理解しやすさを優先しました。callback で各ステップの対数尤度を記録し、収束の様子を後で可視化できるようにしています。
尤度が単調増加することを確認する
まず、合成データ（1次元の潜在から5次元の非線形観測を作ったもの）で、最適化中に尤度が増加することを確かめます。
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng(3)
N = 50
t = np.linspace(-2.5, 2.5, N)
# 1次元潜在 t から非線形に5次元の観測を生成
Y = np.c_[np.sin(t), np.cos(t), t, t**2 / 3, np.sin(2 * t)]
Y += 0.03 * rng.standard_normal(Y.shape)
X, hist, res = fit_gplvm(Y, Q=2, n_iter=120, seed=3)
plt.figure(figsize=(9, 5))
plt.plot(np.arange(len(hist)), hist, marker="o", ms=3, lw=1.6)
plt.xlabel("L-BFGS イテレーション")
plt.ylabel("対数周辺尤度")
plt.title("GP-LVM の尤度最適化")
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()
print(f"最終の負の対数尤度: {res.fun:.2f}")
このコードを実行すると、対数尤度がイテレーションとともに単調に増加し、やがて平坦になって収束する様子が確認できます。前節の理論図（最適化の収束）と同じ振る舞いです。潜在配置とハイパラの同時最適化が、データを最もよく説明する低次元表現へと着実に向かっていることを意味します。非凸問題にもかかわらず、PCA初期化のおかげで安定して収束しています。
digitsデータを2次元に埋め込む
次に、本命の高次元データ可視化です。手書き数字（8×8=64次元）の0〜4を、GP-LVMで2次元に埋め込みます。
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_digits
digits = load_digits()
mask = digits.target < 5 # 0-4 だけ使う
Xd, yd = digits.data[mask], digits.target[mask]
rng = np.random.default_rng(42)
idx = rng.permutation(len(Xd))[:160] # 計算量のため160点に絞る
Xd, yd = Xd[idx] / 16.0, yd[idx] # 画素を [0,1] に正規化
X, hist, res = fit_gplvm(Xd, Q=2, n_iter=150, seed=42)
plt.figure(figsize=(8, 6.5))
sc = plt.scatter(X[:, 0], X[:, 1], c=yd, cmap="tab10",
s=45, edgecolors="k", linewidths=0.4)
plt.colorbar(sc, ticks=range(5), label="数字ラベル")
plt.xlabel("潜在次元 $x_1$"); plt.ylabel("潜在次元 $x_2$")
plt.title("GP-LVM による手書き数字の2次元埋め込み")
plt.tight_layout()
plt.show()
実行すると、64次元の画像が2次元の潜在空間に埋め込まれ、同じ数字どうしが近くに、異なる数字どうしが離れて配置されます。ラベル情報を一切使っていない（教師なし）にもかかわらず、数字ごとのクラスタが自然に形成されているのが見て取れます。GP-LVMが、画像のピクセル空間の複雑な構造を、本質的な少数の自由度に圧縮できている証拠です。
PCAと比較する
最後に、同じデータで線形PCAとGP-LVMの埋め込みを比較し、非線形化の効果を定量的に体感します。
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
digits = load_digits()
mask = digits.target < 5
Xd, yd = digits.data[mask], digits.target[mask]
rng = np.random.default_rng(42)
idx = rng.permutation(len(Xd))[:160]
Xd, yd = Xd[idx] / 16.0, yd[idx]
Xp = PCA(n_components=2, random_state=42).fit_transform(Xd - Xd.mean(0))
Xg, _, _ = fit_gplvm(Xd, Q=2, n_iter=150, seed=42)
fig, axes = plt.subplots(1, 2, figsize=(13, 5.5))
for ax, Xe, title in zip(axes, [Xp, Xg], ["線形PCA", "GP-LVM"]):
sc = ax.scatter(Xe[:, 0], Xe[:, 1], c=yd, cmap="tab10",
s=42, edgecolors="k", linewidths=0.4)
ax.set_title(title); ax.set_xlabel("次元1"); ax.set_ylabel("次元2")
ax.grid(alpha=0.3)
fig.colorbar(sc, ax=axes, ticks=range(5), shrink=0.8, label="数字ラベル")
plt.show()
このコードの出力（本文の比較図に対応）を見ると、PCAでは隣り合う数字のクラスタが重なりがちなのに対し、GP-LVMではクラスタの分離が改善しています。理論で述べた「線形カーネル ⇔ PCA、非線形カーネル ⇔ 非線形次元削減」という関係が、可視化の質の差として現れているのです。同じ教師なし設定でも、写像の柔軟性が表現力を大きく左右することがわかります。
応用
GP-LVMは確率モデルなので、可視化以外にも多彩な応用が開けます。代表的なものを見ておきましょう。
欠損値の補完
GP-LVMは観測 $\boldsymbol{Y}$ の一部が欠損していても、潜在多様体上で欠損次元を予測できます。各観測次元はGPなので、潜在 $\boldsymbol{X}$ を学習したあと、欠損していた次元の値をGP回帰の予測平均で埋められるのです。
この図は、ある2次元観測データの一部（白丸）で $y_2$ が欠損していたとき、GP-LVMが学習した潜在多様体を使って $y_2$ を予測した結果（赤いバツ）です。補完された値が真値（白丸の位置）にほぼ重なっており、潜在多様体が欠損情報をうまく補っていることがわかります。これは推薦システムやセンサデータの欠測補完に直結する能力です。
新しいデータの生成
潜在空間の任意の点 $\boldsymbol{x}_*$ を選び、各次元のGPで $f_d(\boldsymbol{x}_*)$ を予測すれば、新しい観測データを生成できます。たとえば手書き数字を学習したGP-LVMなら、潜在空間を連続的に動かすことで、ある数字から別の数字へ滑らかに変形する画像列を作れます。これは生成モデルとしての顔も持つことを意味し、後の深層生成モデル（VAEなど）の先駆けとも言えます。
動作データの多様体学習
GP-LVMが特に有名になったのは、人間の動作データ（モーションキャプチャ）のモデル化です。歩行や走行といった動作は、数十関節の高次元データですが、本質的には「歩行サイクルの位相」のような少数のパラメータで支配されています。GP-LVMはこの低次元の動作多様体を発見できます。
左がS字状の多様体として観測された高次元の動作軌跡、右がGP-LVMで展開した2次元の潜在表現です。連続的な動作パラメータが、低次元の滑らかな軌跡として復元されています。コンピュータグラフィックスのアニメーション生成や、ロボットの動作計画でこうした表現が活用されています。
これらの応用に共通するのは、「観測を生成する確率過程をモデル化している」という点です。単なる座標変換にとどまらず、不確実性・補完・生成までを統一的に扱えるのが、確率モデルとしてのGP-LVMの強みです。
まとめ
本記事では、ガウス過程潜在変数モデル（GP-LVM）について、確率的PCAから出発して理論の導出・実装・応用まで解説しました。
- モチベーション: PCAは線形写像なので曲がった多様体をほどけない。非線形次元削減を確率モデルとして定式化したい。
- PPCAの復習: 潜在 $\boldsymbol{x}_n$ から線形写像 $\boldsymbol{W}\boldsymbol{x}_n + \boldsymbol{\epsilon}$ で観測を生成し、$\boldsymbol{x}$ を周辺化すると $\boldsymbol{y}_n \sim \mathcal{N}(\boldsymbol{\mu}, \boldsymbol{W}\boldsymbol{W}^\top + \sigma^2\boldsymbol{I})$。最尤解はPCA。
- GP-LVMの定式化: PPCAの線形写像をガウス過程に置換。各観測次元 $\boldsymbol{y}_{:,d} \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{K})$ で、$\boldsymbol{K}$ は潜在 $\boldsymbol{X}$ から決まるカーネル行列。
- 尤度と最適化: 対数尤度 $\mathcal{L} = -\frac{D}{2}\log|\boldsymbol{K}| – \frac{1}{2}\mathrm{tr}(\boldsymbol{K}^{-1}\boldsymbol{Y}\boldsymbol{Y}^\top)$ を $\boldsymbol{X}$ について勾配法（L-BFGS）で最大化。PCA初期化が定石。
- PCAとの関係: 線形カーネルのGP-LVMはPPCA（PCA）に一致する。「特徴の共分散」と「サンプルの共分散」の双対関係であり、カーネルを非線形にするだけで非線形次元削減に拡張できる。
- ベイズGP-LVM: 潜在 $\boldsymbol{X}$ にも事前分布を置き変分推論で積分消去（Titsias & Lawrence 2010）。ARDカーネルにより潜在次元数を自動決定できる。
- 応用: 高次元データ可視化、欠損補完、データ生成、動作データの多様体学習。
数値実験では、GP-LVMの尤度がL-BFGSで単調に収束すること、手書き数字データの2次元埋め込みが線形PCAより明瞭にクラスタを分離することを確認しました。GP-LVMは、ガウス過程という確率的・非パラメトリックな枠組みを次元削減に持ち込んだエレガントな手法であり、その後のベイズ多様体学習や深層生成モデルへとつながる重要な礎石です。
次のステップとして、以下の記事も参考にしてください。