GPyTorch ではじめる深層ガウス過程入門

ガウス過程（GP: Gaussian Process）は、関数そのものに確率分布を定義するノンパラメトリックなモデルです。このモデルの最大の強みは、単なる予測値だけでなく、その不確実性（信頼区間）を定量的に示せる点にあります。これにより、予測の信頼性が求められる場面で、非常に有用なアプローチとなります。

しかし、従来のガウス過程の実装には、大規模データに対する計算コスト（通常 $O(n^3)$）が課題になります。そこで注目されるのが、PyTorch上に構築されたスケーラブルなGP推論ライブラリ GPyTorch です。GPyTorchは、BBMM（Blackbox Matrix-Matrix multiplication）などの手法を用いることで、この計算上のボトルネックを大幅に解消し、柔軟でカスタム可能なモデル構築を可能にします。

今回は、まずガウス過程モデルの解説から始めます。それから、GPyTorchの概要を解説し、具体的な実装コードを通して深く掘り下げます。この強力なフレームワークを使いこなし、複雑なデータに対しても、不確実性を伴う信頼性の高い予測を実現する方法を探求しましょう。

ガウス過程とは

ガウス過程は、関数そのものに確率分布を定義するノンパラメトリックなモデルです。これは、機械学習において広く使われるアプローチであり、単に入力点 $x$ に対する出力 $y$ を予測するだけでなく、関数全体の集合に対して確率的な振る舞いを定義します。

**図1. ガウス過程回帰の例（https://arxiv.org/pdf/1805.07152 より引用）**

確率過程としての定義

ガウス過程の定義は非常に簡潔です。任意の有限個の入力点 $X = {x_1, \dots, x_n}$ を選んだとき、それらの点における関数値 $f(X) = {f(x_1), \dots, f(x_n)}$ が、必ず多変量ガウス分布 (Multivariate Gaussian Distribution) に従う、という性質を持つ確率過程として定義されます。

この関数分布は、次の2つの関数によって特徴づけられます。

平均関数 $\mu(x)$：関数の平均的な形状（期待値）を定義します。
共分散関数（カーネル関数） $k(x, x’)$：任意の2つの入力点 $x$ と $x’$ の間の「類似性」や「相関」を定義します。

したがって、ガウス過程は以下のように記述されます。

$$ f(x) \sim \mathcal{GP}(m(x), k(x, x’)) $$

特にカーネル関数は、予測される関数の滑らかさや周期性といった基本的な性質を決定する、モデルの根幹とも言える部分です。

特徴と応用分野

ガウス過程はノンパラメトリックモデルに分類されます。これは、モデルの複雑さ（実質的なパラメータ数）がデータの量に比例して柔軟に変化する（増大する）ため、データの構造に合わせた精度の高い表現が可能です。

ガウス過程の主要な応用分野は多岐にわたります。

回帰（ガウス過程回帰, GPR: Gaussian Process Regression）
分類（GPC: Gaussian Process Classification） （非ガウス尤度を用いる）
ベイズ最適化 (Bayesian Optimization) や不確実性推定 (Uncertainty Estimation)

利点と克服すべき課題

利点

ガウス過程の最大の魅力は、予測値と同時に、その予測に対する不確実性（信頼区間）を定量的に提供できる点です。予測の信頼性が求められる医療や金融、あるいはベイズ最適化のような探索的決定が必要な分野において、非常に強力な情報となります。また、少ない訓練データ量でも、カーネルによって定義された事前分布の恩恵を受け、高い汎化性能を発揮しやすい傾向があります。

課題

一方で、ガウス過程には従来から大きな計算上の課題があります。訓練データ数 $n$ が増加すると、予測分布の計算に必要なカーネル行列の逆行列計算や対数行列式の計算に、計算量が $O(n^3)$ 必要となります。この$O(n^3)$ のボトルネックと、カーネル行列の格納に $O(n^2)$ のメモリが必要となるため、ガウス過程はこれまで数千点以下のデータセットに限定されてきました。このスケーラビリティの問題を解決するために、様々な近似手法や、GPUを効率的に利用する手法が研究されてきました。

GPyTorch とは

GPyTorch は、PyTorchのエコシステム上に構築された、柔軟性とスケーラビリティを両立させたガウス過程ライブラリです。従来のガウス過程実装が抱えていた、大規模データに対する計算コスト $O(n^3)$ の課題を克服し、最新のGPUハードウェアを最大限に活用することを目指しています。

PyTorchとの密接な連携と設計思想

GPyTorchの設計は、従来のガウス過程パッケージとは一線を画しています。このライブラリは、完成されたガウス過程モデルを提供するのではなく、ユーザーが必要なコンポーネントを組み合わせてカスタムモデルを効率的に構築するためのツールを提供します。PyTorchで深層学習のモデルを構築するアプローチと類似しており、モデル設計に高い柔軟性をもたらします。

GPyTorchのモデルは、PyTorchの torch.nn.Module を直接拡張しています。そのため、ガウス過程のハイパーパラメータ（カーネルの長さスケールやノイズ分散など）は、訓練可能な torch.nn.Parameter として扱われます。これらのハイパーパラメータの最適化（通常、タイプII最尤推定, Type-II MLEを用います）には、Adamなどの標準的なPyTorchオプティマイザ（torch.optim）を利用できます。

ガウス過程モデルを構成する主要コンポーネント

GPyTorchにおいて、厳密な推論（Exact Inference）を行うガウス過程モデルを構築する際には、主に以下の5つのオブジェクトを定義する必要があります。

コンポーネント	クラス/モジュール	役割
ガウス過程モデル	gpytorch.models.ExactGP	推論処理の中核を担い、学習データに条件付けられた潜在ガウス過程の事後分布を返します。
尤度 (Likelihood)	gpytorch.likelihoods.GaussianLikelihood	最も一般的なガウス過程回帰の尤度です。観測ノイズが均一である（誤差項の分散が常に一定）と仮定します。
平均 (Mean)	gpytorch.means.ConstantMean	ガウス過程の事前平均を定義します。多くの場合、ConstantMean（定数平均）が推奨されます。
カーネル (Kernel)	gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel	ガウス過程の事前共分散（共分散関数）を定義します。RBFカーネルに出力スケールを適用したものが推奨される出発点です。
分布	gpytorch.distributions.MultivariateNormal	多変量ガウス分布を表現するために用いられます。

大規模データへの対応とスケーラビリティ技術

GPyTorchは、GPUを最大限に活用し、大規模データセットでもガウス過程推論を可能にする複数のスケーラビリティ技術を組み込んでいます。

BBMM (Blackbox Matrix-Matrix inference): GPyTorchの中心的な推論エンジンであり、行列-行列乗算を活用することで、厳密なガウス過程推論の計算複雑度を$O(n^3)$から$O(n^2)$に削減します。これは、Cholesky分解よりもGPUアクセラレーションに適しています。
GPU/Multi-GPU活用: ユーザーがデータとモデルをGPUに配置するだけで、GPU上での計算が自動的に有効化されます。さらに、MultiDeviceKernelを使用することで、大規模データセットに対し複数のGPUにカーネル行列を分散させた推論も可能です。
SKI/KISS-GP: Structured Kernel Interpolation (SKI)、またはKISS-GPと呼ばれる手法は、カーネルをGridInterpolationKernelでラップするだけで、10万点以上のデータセットにも対応可能な近似手法です。この方法は漸近的にほぼ線形時間でスケールするという利点があります。
LOVE (Lanczos Variance Estimates): 予測分散を高速に計算するための機能です。特にKISS-GPと組み合わせることで、gpytorch.settings.fast_pred_varコンテキスト内で、予測平均と分散の両方を定数時間で計算できる場合があります。
KeOps統合: 外部ライブラリであるKeOpsとの統合により、GPU上でカーネル行列乗算を効率的に実行できます。これにより、数10万点のデータポイントに対する厳密なガウス過程訓練を、単一のGPUで短時間で実行することが可能になります。

GPyTorchの使い方

このセクションでは、GPyTorchを使った実践的なモデル構築の方法を、具体的なコードを交えながらステップバイステップで解説します。まずは基本的なガウス過程回帰モデルから始め、次に応用例として深層カーネル学習（DKL）へと進みます。

セットアップ

はじめに、GPyTorchと関連ライブラリをインストールします。PyTorchが必須なので、環境に合わせて先にインストールしておきましょう。

# 使用したいバージョンのPyTorchがあれば先にインストールします
pip install torch==2.6.0 torchvision==0.21.0 torchaudio==2.6.0 --index-url https://download.pytorch.org/whl/cu124

# GPyTorch本体をインストールします
pip install gpytorch

# プロット用にmatplotlibもインストールします
pip install matplotlib

シンプルなガウス過程回帰 (Exact GP Regression)

ここでは、GPyTorchの最も基本的な機能を用いて、厳密なガウス過程回帰モデルを構築します。sin関数から生成したデータにモデルを適合させ、予測の平均と不確実性（信頼区間）を可視化するまでの一連の流れを見ていきましょう。

ステップ1: データの準備

まず、モデルの訓練に使用するデータを生成します。今回は、sin関数にガウスノイズを加えたシンプルなデータセットを用意します。

import math
import torch
import gpytorch
from matplotlib import pyplot as plt

# 訓練データを作成します
# torch.linspaceで0から1の範囲で100点のデータを生成します
train_x = torch.linspace(0, 1, 100)
# 真の関数（sin関数）に、標準偏差sqrt(0.04)=0.2のガウスノイズを加えます
train_y = torch.sin(train_x * (2 * math.pi)) + torch.randn(train_x.size()) * math.sqrt(0.04)

ステップ2: ガウス過程モデルの定義

次に、ガウス過程モデルを定義します。GPyTorchでは、gpytorch.models.ExactGPを継承してカスタムモデルを構築するのが一般的です。__init__でモデルの構成要素（平均関数、共分散関数）を定義し、forwardで入力xに対するガウス過程の事前分布を返します。

# gpytorch.models.ExactGPを継承して、独自のGPモデルを作成します
class ExactGPModel(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood):
        super(ExactGPModel, self).__init__(train_x, train_y, likelihood)
        
        # 平均関数を定義します。ここでは定数平均（ConstantMean）を使用します。
        self.mean_module = gpytorch.means.ConstantMean()
        
        # 共分散関数（カーネル）を定義します。
        # RBFカーネルは滑らかな関数を仮定する一般的な選択肢です。
        # ScaleKernelは、カーネル全体に出力スケールを適用します。
        self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())

    # forwardメソッドで、GPの事前分布を定義します
    def forward(self, x):
        # 平均と共分散を計算
        mean_x = self.mean_module(x)
        covar_x = self.covar_module(x)
        # 多変量正規分布としてガウス過程を返します
        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)

# 観測ノイズが正規分布に従うと仮定する、最も一般的な尤度（GaussianLikelihood）を定義します
likelihood = gpytorch.likelihoods.GaussianLikelihood()
# モデルをインスタンス化します
model = ExactGPModel(train_x, train_y, likelihood)

ステップ3: モデルの訓練

GPyTorchのモデルはPyTorchのnn.Moduleなので、訓練ループもPyTorchの標準的なスタイルで行います。オプティマイザ（例: Adam）を定義し、損失関数として周辺対数尤度（Marginal Log-Likelihood）を最大化（つまり、負の周辺対数尤度を最小化）することで、カーネルのハイパーパラメータなどを学習します。

# 訓練モードに設定
model.train()
likelihood.train()

# Adamオプティマイザを使用し、モデルの全パラメータ（GPのハイパーパラメータ）を渡します
optimizer = torch.optim.Adam(model.parameters(), lr=0.1)

# 損失関数として、厳密なGP向けの周辺対数尤度（ExactMarginalLogLikelihood）を使用します
mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)

training_iter = 50
for i in range(training_iter):
    optimizer.zero_grad()
    # モデルから出力（多変量正規分布）を取得
    output = model(train_x)
    # 負の周辺対数尤度を損失として計算
    loss = -mll(output, train_y)
    loss.backward()
    
    # パラメータの学習状況をモニタリング
    print('Iter %d/%d - Loss: %.3f   lengthscale: %.3f   noise: %.3f' % (
        i + 1, training_iter, loss.item(),
        model.covar_module.base_kernel.lengthscale.item(),
        model.likelihood.noise.item()
    ))
    optimizer.step()

ステップ4: 予測と結果の可視化

訓練が完了したら、モデルを評価モードに切り替えて、新しいデータ点に対する予測を行います。likelihood(model(test_x))を呼び出すことで、予測分布（平均と分散を含む）が簡単に得られます。ガウス過程の最大の特長である「不確実性」を信頼区間としてプロットしてみましょう

# 評価モードに設定
model.eval()
likelihood.eval()

# `torch.no_grad()`と`gpytorch.settings.fast_pred_var()`コンテキスト内で予測を実行
with torch.no_grad(), gpytorch.settings.fast_pred_var():
    test_x = torch.linspace(0, 1, 51)
    # likelihood(model(test_x)) で、観測ノイズも考慮した予測分布を取得します
    observed_pred = likelihood(model(test_x))

# 結果をプロットします
with torch.no_grad():
    plt.style.use("ggplot")
    f, ax = plt.subplots(1, 1, figsize=(6, 5))

    # 95%信頼区間の上下限を取得
    lower, upper = observed_pred.confidence_region()
    # 訓練データを黒い点でプロット
    ax.plot(train_x.numpy(), train_y.numpy(), 'k*')
    # 予測平均を青線でプロット
    ax.plot(test_x.numpy(), observed_pred.mean.numpy(), 'b')
    # 信頼区間を影付きで描画
    ax.fill_between(test_x.numpy(), lower.numpy(), upper.numpy(), alpha=0.5)
    ax.set_ylim([-3, 3])
    ax.legend(['Observed Data', 'Mean', 'Confidence'])
    plt.title("Simple Gaussian Process Regression")
    plt.xlabel("$x$")
    plt.ylabel("$y = \\sin(x)$")
    plt.grid()
    plt.show()

以下のプロット結果では、全領域について信頼区間の大きさがほぼ一定ですが、一般的に訓練データが密な領域では信頼区間が狭く、データがない領域では広がる傾向になります。

Deep Kernel Learning (DKL) の実装

次に、ニューラルネットワークを特徴抽出器としてガウス過程回帰に組み込む深層カーネル学習（DKL）を実装します。これにより、複雑な非線形パターンを持つデータに対しても、ガウス過程の表現力を大幅に向上させることができます。

基本的な流れは、ニューラルネットワークで入力データを低次元の特徴空間に射影し、その特徴空間上でガウス過程を適用するというものです。GPyTorchでは、このエンドツーエンドのモデルをシームレスに構築し、ニューラルネットワークの重みとガウス過程のハイパーパラメータを同時に最適化できます。

ステップ1: 特徴抽出器とDKLモデルの定義

まず、入力データを変換するための小規模のニューラルネットワーク（特徴抽出器）をtorch.nn.Sequentialで定義します。その後、DKL-GPモデルを定義します。このモデルは内部に特徴抽出器を持ち、forwardメソッドの中で、入力データをニューラルネットワークに通した結果をガウス過程に渡す、という処理を行います。

# DKLでは、入力データを低次元の特徴空間に射影するNNを定義します
# この例では、1次元の入力を2次元の特徴量に変換します
feature_extractor = torch.nn.Sequential(
    torch.nn.Linear(1, 100),
    torch.nn.ReLU(),
    torch.nn.Linear(100, 50),
    torch.nn.ReLU(),
    torch.nn.Linear(50, 2),
)

# DKL-GPモデルの定義
class DKLGPModel(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood):
        super(DKLGPModel, self).__init__(train_x, train_y, likelihood)
        self.mean_module = gpytorch.means.ConstantMean()
        
        # 特徴抽出器をモデルに組み込みます
        self.feature_extractor = feature_extractor
        
        # 共分散モジュール: 特徴空間の次元（この場合は2次元）に合わせてカーネルを定義します
        # SKI/KISS-GP (GridInterpolationKernel) を使うことで、多次元空間でも計算を高速化します
        self.covar_module = gpytorch.kernels.GridInterpolationKernel(
            gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel(ard_num_dims=2)),
            num_dims=2, grid_size=100
        )

    def forward(self, x):
        # 1. 入力xを特徴抽出器に通して射影します
        projected_x = self.feature_extractor(x)
        
        # 2. 射影された特徴量をGPに渡します
        mean_x = self.mean_module(projected_x)
        covar_x = self.covar_module(projected_x)
        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)

likelihood = gpytorch.likelihoods.GaussianLikelihood()
model = DKLGPModel(train_x.unsqueeze(-1), train_y, likelihood)

ステップ2: モデルの訓練と予測

訓練と予測のコードは、先ほどのシンプルなガウス過程回帰とほとんど同じです。大きな違いは、optimizerがニューラルネットワークのパラメータとガウス過程のハイパーパラメータの両方を同時に学習対象とすることです。これにより、データに最もフィットする特徴表現そのものを学習することができます。

# 訓練と予測のプロセスをまとめて実行します

# 訓練モードに設定
model.train()
likelihood.train()

# オプティマイザにはNNとGPの全パラメータが渡されます
optimizer = torch.optim.Adam(model.parameters(), lr=0.01)
mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)

# 訓練ループ
training_iter = 50
for i in range(training_iter):
    optimizer.zero_grad()
    # DKLモデルは(N, D_in)の形式で入力を受け取るため、unsqueeze(-1)で次元を追加
    output = model(train_x.unsqueeze(-1))
    loss = -mll(output, train_y)
    loss.backward()
    print(f'Iter {i + 1}/{training_iter} - Loss: {loss.item():.3f}')
    optimizer.step()

# 評価モードに設定
model.eval()
likelihood.eval()

# テストデータで予測とプロット
test_x = torch.linspace(0, 1, 101)
with torch.no_grad(), gpytorch.settings.fast_pred_var():
    observed_pred = likelihood(model(test_x.unsqueeze(-1)))

with torch.no_grad():
    f, ax = plt.subplots(1, 1, figsize=(6, 5))
    lower, upper = observed_pred.confidence_region()
    ax.plot(train_x.numpy(), train_y.numpy(), 'k*')
    ax.plot(test_x.numpy(), observed_pred.mean.numpy(), 'b')
    ax.fill_between(test_x.numpy(), lower.numpy(), upper.numpy(), alpha=0.5)
    ax.set_ylim([-3, 3])
    ax.legend(['Observed Data', 'Mean', 'Confidence'])
    plt.show()

おわりに

本記事では、ガウス過程が持つ不確実性の定量化という強力な利点と、その課題であったスケーラビリティを克服するライブラリ、GPyTorchの概要と基本機能について解説してきました。

GPyTorchは、PyTorchのモジュール構造を継承することで、従来のガウス過程パッケージにはなかった柔軟なモデル構築を可能にしています。そして、BBMM（Blackbox Matrix-Matrix multiplication）や KISS-GP（Structured Kernel Interpolation）といったスケーラビリティ技術、およびGPUアクセラレーションを活用することにより、ガウス過程を数万から数10万点規模のデータセットに適用することを可能にしました。

GPyTorchの持つ柔軟なエコシステムは、今後さらに高度な応用を可能にします。

不確実性の定量化: ガウス過程の特性上、予測と同時に 予測の不確実性（分散）を自然に提供します。これにより、予測の信頼性が実務上求められる、意思決定支援の分野で、極めて信頼性の高い情報を提供できます。
高度なカーネル設計: ユーザーは gpytorch.kernels.Kernel クラスを継承し、自作カーネルを容易に実装できます。
複合カーネル構造: 複数のカーネルを AdditiveKernel（加法構造）や ProductKernel（積構造）として組み合わせることが可能です。これにより、モデルの表現力を高めたり、予測の解釈性を向上させたりすることができます。
マルチタスク学習: 複数の関連タスク間の相関関係を捉える MultitaskKernel を使用したマルチタスクガウス過程モデルの構築をサポートしています。
フルベイズ推論: Pyroライブラリと統合することで、ハイパーパラメータに対してNUTSサンプリングを用いた完全なベイズ推論（Fully Bayesian Inference）を実行でき、ハイパーパラメータの不確実性も考慮したモデリングが可能です。

GPyTorchは、単なる予測モデルとしてだけでなく、不確実性情報が不可欠なベイズ最適化や、ニューラルネットワークと融合させた深層カーネル学習（DKL）の研究開発基盤としても活用されています。ぜひこの強力なライブラリを通じて、ガウス過程の可能性を探求してください。

More Information

arXiv:1511.02222, Andrew Gordon Wilson et al., 「Deep Kernel Learning」, https://arxiv.org/abs/1511.02222

arXiv:1611.00336, Andrew Gordon Wilson et al., 「Stochastic Variational Deep Kernel Learning」, https://arxiv.org/abs/1611.00336

arXiv:1803.06058, Geoff Pleiss et al., 「Constant-Time Predictive Distributions for Gaussian Processes」, https://arxiv.org/abs/1803.06058

arXiv:1809.11165, Jacob R. Gardner et al., 「GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration」, https://arxiv.org/abs/1809.11165

arXiv:2006.11267, Geoff Pleiss et al., 「Fast Matrix Square Roots with Applications to Gaussian Processes and Bayesian Optimization」, https://arxiv.org/abs/2006.11267

arXiv:2007.10731, Arthur Leroy et al., 「MAGMA: Inference and Prediction with Multi-Task Gaussian Processes」, https://arxiv.org/abs/2007.10731

codemajinのえんとろぴぃ

Blog

GPyTorch ではじめる深層ガウス過程入門

ガウス過程とは

確率過程としての定義

特徴と応用分野

利点と克服すべき課題

利点

課題

GPyTorch とは

PyTorchとの密接な連携と設計思想

ガウス過程モデルを構成する主要コンポーネント

大規模データへの対応とスケーラビリティ技術

GPyTorchの使い方

セットアップ

シンプルなガウス過程回帰 (Exact GP Regression)

ステップ1: データの準備

ステップ2: ガウス過程モデルの定義

ステップ3: モデルの訓練

ステップ4: 予測と結果の可視化

Deep Kernel Learning (DKL) の実装

ステップ1: 特徴抽出器とDKLモデルの定義

ステップ2: モデルの訓練と予測

おわりに

Blog

GPyTorch ではじめる深層ガウス過程入門

ガウス過程とは

確率過程としての定義

特徴と応用分野

利点と克服すべき課題

利点

課題

GPyTorch とは

PyTorchとの密接な連携と設計思想

ガウス過程モデルを構成する主要コンポーネント

大規模データへの対応とスケーラビリティ技術

GPyTorchの使い方

セットアップ

シンプルなガウス過程回帰 (Exact GP Regression)

ステップ1: データの準備

ステップ2: ガウス過程モデルの定義

ステップ3: モデルの訓練

ステップ4: 予測と結果の可視化

Deep Kernel Learning (DKL) の実装

ステップ1: 特徴抽出器とDKLモデルの定義

ステップ2: モデルの訓練と予測

おわりに

関連記事

LoRAによる大規模言語モデルのファインチューニング

UnslothではじめるLLMのFine-tuning

LiteASR: 低ランク近似による効率的な自動音声認識の実現