制約付きガウス過程回帰 — 物理法則や単調性をGPに「正しく」入れる方法

ガウス過程回帰(GPR)は、少ないデータから「予測値+不確かさ」を返してくれる便利な道具です。でも、素のGPを使っているとこんな不満が出てきませんか。

  • 流体の速度場を回帰したら、質量保存($\nabla\cdot\bm{v}=0$)を破る速度場が出てきた
  • 吸着量や効用みたいに「増えることはあっても減らない」と分かっている量なのに、データが少ないと予測が波打って一部で減少してしまう
  • 濃度や確率みたいに「絶対に負にならない」量なのに、予測の信頼区間がマイナスまで伸びる

私たちは「物理法則」や「単調性」という事前知識を持っているのに、素のGPはそれを知りません。だったら、その知識をGPに教え込めないか?——これが制約付きガウス過程回帰の動機です。

そして本記事で一番伝えたいのは、こういう問いへの答えです。

なぜ、線形の制約(微分や発散ゼロ)は「カーネルに織り込むだけ」で、サンプルも予測も全部、全域で厳密に満たすようにできるのか? 一方で、なぜ単調性や正値性のような不等式制約は「近似」でしか入れられないのか?

この違いの正体は、GPが持つ「線形作用素で閉じている」という性質にあります。そこを演算レベルで腹落ちさせるのが本記事のゴールです。制約をうまく入れると、データが少なくても予測が事前知識に支えられて安定し、外挿でも暴れにくくなります。応用は広く、非圧縮流れ・磁場・電場の回帰、計算機実験のサロゲートモデル、物理を尊重するベイズ最適化などで効きます。

本記事の内容

  • 制約を「線形(等式)」と「不等式・形状」の2系統に整理する
  • GPが線形作用素で閉じている、という核心の性質
  • 微分(傾き)観測・発散ゼロ場を厳密に入れる方法とPython実装
  • 単調性を近似で入れる方法(仮想微分点+probit尤度)とPython実装

前提知識

この記事はGPRの基礎の上に立ちます。以下を先に読むと、本記事の話がスッと入ってきます。

制約は大きく2系統に分かれる

最初に地図を持っておきましょう。GPに入れられる制約は、扱いやすさで2つに大別できます。

GPに入れられる制約の2系統マップ:線形等式制約は厳密に全域で成立、不等式・形状制約は近似で点ごとに促す

  • ① 線形(等式)制約:微分値の観測、勾配の観測、発散ゼロ・回転ゼロの場、線形PDE、境界条件など。これらはカーネルに織り込めて、サンプルも予測も全域で厳密に満たす。しかもガウス性が保たれるので、普通のGPと同じ閉形式で解けます。
  • ② 不等式・形状制約:単調性 $f’\ge0$、凸性 $f”\ge0$、正値・上下限 $a\le f\le b$ など。これらはガウス性を壊すため、仮想点と確率的な尤度を使って「点ごとに、近似的に」促すしかありません。推論には EP(期待値伝播)や MCMC が要ります。

この記事の構成も、この地図に沿います。まずは「なぜ①は厳密にできるのか」を生む、GPの一番大事な性質から始めます。図の左側(厳密に守れる方)を理解するための土台です。

核心:GPは線形作用素で「閉じている」

制約付きGPのすべての出発点が、この一文です。

ガウス過程に線形作用素を作用させても、結果はまたガウス過程になる。

ここで線形作用素 $\mathcal{L}$ とは、微分 $\partial/\partial x$、積分、定数倍、和——「線形性 $\mathcal{L}(af+bg)=a\mathcal{L}f+b\mathcal{L}g$ を満たす操作」のことです。

ガウス過程は線形作用素で閉じている:gにL(微分・積分・線形結合)を作用させてもGPのまま

なぜ閉じるのか。直感的には、ガウス分布が線形変換で閉じている(ガウスを定数倍して足してもガウス)ことの、無限次元版だと思えば腑に落ちます。GPは「どんな有限個の点を取り出しても多変量ガウスになる」確率過程でした。微分も「近い2点の差の極限」という線形操作なので、ガウスの線形結合の極限、つまりやっぱりガウスなのです。

大事なのは、変換後のGPの共分散が閉形式で書けることです。$g\sim\mathcal{GP}(0,k_g)$ に作用素 $\mathcal{L}$ をかけると、

$$ \begin{equation} \mathrm{Cov}[\mathcal{L}g(x),\,\mathcal{L}g(x’)] = \mathcal{L}_x\,\mathcal{L}_{x’}\,k_g(x,x’) \end{equation} $$

つまり「カーネルを $x$ と $x’$ の両方で同じ作用素にかける」だけ。微分の制約なら、カーネルを微分するだけで新しいカーネルが手に入るのです。

この性質さえ認めれば、線形制約GPは全部その応用です。まず一番素朴な「微分(傾き)を観測に混ぜる」ところから手を動かしてみましょう。

線形制約①:微分(傾き)を観測する

普通のGPは関数の $f(x_i)$ を観測します。でも、線形作用素で閉じているなら、傾き $f'(x_i)$ も同じ枠組みで観測値として混ぜられます。導関数 $f’$ もGPで、$f$ との間に共分散があるからです。

RBFカーネル $k(x,x’)=\sigma^2\exp\!\big(-\frac{(x-x’)^2}{2\ell^2}\big)$ で、必要な共分散を式(1)に従って計算します。

$$ \frac{\partial k}{\partial x’} = \frac{x-x’}{\ell^2}\,k,\qquad \frac{\partial^2 k}{\partial x\,\partial x’} = \Big(\frac{1}{\ell^2}-\frac{(x-x’)^2}{\ell^4}\Big)k $$

1つ目は「値 $f(x)$ と傾き $f'(x’)$ の共分散」、2つ目は「傾きどうし $f'(x),f'(x’)$ の共分散」です。これらをブロックに並べて同時共分散行列を作り、値の観測 $\bm{y}$ と傾きの観測 $\bm{y}’$ をまとめて条件付ければよいだけです。

import numpy as np

def rbf(X, Y, l, s=1.0):
    d = X[:, None] - Y[None, :]
    return s**2 * np.exp(-d**2 / (2*l**2))

def dk_dy(X, Y, l, s=1.0):   # ∂k/∂y = (x-y)/l^2 * k  (値と傾きの共分散)
    d = X[:, None] - Y[None, :]
    return (d / l**2) * rbf(X, Y, l, s)

def dk_dx(X, Y, l, s=1.0):   # ∂k/∂x = -(x-y)/l^2 * k
    return -dk_dy(X, Y, l, s)

def d2k(X, Y, l, s=1.0):     # ∂^2k/∂x∂y  (傾きどうしの共分散)
    d = X[:, None] - Y[None, :]
    return (1/l**2 - d**2/l**4) * rbf(X, Y, l, s)

l, s, noise = 1.2, 1.0, 1e-4
f  = lambda x: np.sin(x)
df = lambda x: np.cos(x)
Xtr = np.array([-3.0, 0.0, 3.0])           # 観測点(たった3点)
Xs  = np.linspace(-5, 5, 200)

def predict(use_deriv):
    if not use_deriv:                       # 値だけ観測(普通のGP)
        K  = rbf(Xtr, Xtr, l, s) + noise*np.eye(len(Xtr))
        Ks = rbf(Xs, Xtr, l, s)
        mean = Ks @ np.linalg.solve(K, f(Xtr))
        cov  = rbf(Xs, Xs, l, s) - Ks @ np.linalg.solve(K, Ks.T)
    else:                                   # 値+傾きを観測
        Kff, Kfd = rbf(Xtr, Xtr, l, s), dk_dy(Xtr, Xtr, l, s)
        Kdf, Kdd = dk_dx(Xtr, Xtr, l, s), d2k(Xtr, Xtr, l, s)
        K = np.block([[Kff, Kfd], [Kdf, Kdd]]) + noise*np.eye(2*len(Xtr))
        y = np.concatenate([f(Xtr), df(Xtr)])        # 値と傾きを並べる
        Ks = np.hstack([rbf(Xs, Xtr, l, s), dk_dy(Xs, Xtr, l, s)])
        mean = Ks @ np.linalg.solve(K, y)
        cov  = rbf(Xs, Xs, l, s) - Ks @ np.linalg.solve(K, Ks.T)
    sd = np.sqrt(np.clip(np.diag(cov), 0, None))
    return mean, sd

_, sd0 = predict(False)
_, sd1 = predict(True)
print(f"平均±2σ幅: 値だけ {sd0.mean():.3f} → 値+傾き {sd1.mean():.3f}")
# -> 値だけ 0.550 → 値+傾き 0.277

微分観測ありなしの比較:傾きも観測すると事後分布が締まり真の関数によく一致する

左(値だけ)と右(値+傾き)を比べてください。同じ3点を観測していても、各点で傾き(橙の線分)を一緒に教えると、事後分布が目に見えて締まります。実際、平均的な不確かさ(±2σ幅)は 0.550 から 0.277 へと半分近くに縮みました。傾きという1次の情報が、点と点の間の振る舞いを強く拘束しているのです。

ここでのポイントは、傾きの観測が「制約」として厳密に効いていること。$f’$ もGPの一員だから、値の観測とまったく同じ条件付けで扱えました。では、もっと物理的な制約——「発散がゼロ」のような場の法則——も同じ精神で入れられるでしょうか。ここからが本記事の山場です。

線形制約②:発散ゼロの場を「構造的に」入れる

非圧縮の流れや磁場は、$\nabla\cdot\bm{f}=0$(発散ゼロ)という法則を満たします。2次元なら

$$ \frac{\partial f_1}{\partial x_1}+\frac{\partial f_2}{\partial x_2}=0 $$

です。これをGPに入れたい。ナイーブには「2つの成分 $f_1,f_2$ を別々のGPで回帰する」ところですが、それだと両者が無関係に動くので、発散ゼロはまず満たされません。

ここで Jidling らの線形制約GP(NeurIPS 2017)の発想が効きます。アイデアはこうです。目標の場 $\bm{f}$ を、ある下地のスカラーGP $g$ を線形作用素 $G$ で変換したものとして定義します。

$$ \begin{equation} \bm{f}(\bm{x}) = G_{\bm{x}}\,g(\bm{x}) \end{equation} $$

制約を $F_{\bm{x}}\bm{f}=0$(発散ゼロなら $F_{\bm{x}}=[\partial/\partial x_1,\ \partial/\partial x_2]$)と書くと、私たちが望むのは「どんな $g$ を選んでも $F_{\bm{x}}\bm{f}=0$ になる」ことです。式(2)を代入すると $F_{\bm{x}}\bm{f}=F_{\bm{x}}G_{\bm{x}}\,g$。これが任意の $g$ でゼロになる条件は、

$$ \begin{equation} F_{\bm{x}}\,G_{\bm{x}} = 0 \end{equation} $$

作用素として $F G=0$ であること。言い換えると、$G$ の列が「制約作用素 $F$ の零空間(ヌルスペース)を張る」ことです。これさえ満たせば、$\bm{f}=Gg$ と書ける場はもれなく発散ゼロになります。

2次元の発散ゼロでは、次の $G$ が条件(3)を満たします。

$$ G_{\bm{x}}=\begin{bmatrix}-\,\partial/\partial x_2\\[2pt]\partial/\partial x_1\end{bmatrix} \quad\Rightarrow\quad F_{\bm{x}}G_{\bm{x}}=-\frac{\partial^2}{\partial x_1\partial x_2}+\frac{\partial^2}{\partial x_2\partial x_1}=0 $$

偏微分の順序交換でぴったり打ち消し合ってゼロ。これは物理でおなじみ、「スカラーポテンシャル $g$ の回転(カール)を取った場は発散ゼロ」という事実そのものです($\bm{f}=\nabla\times(g\,\hat{\bm{z}})$ の2次元版)。

下地GPを $g\sim\mathcal{GP}(0,k_g)$ とすれば、式(1)より目標場 $\bm{f}$ の共分散($2\times2$ 行列)は

$$ \begin{equation} K_{\bm{f}}(\bm{x},\bm{x}’)=G_{\bm{x}}\,k_g(\bm{x},\bm{x}’)\,G_{\bm{x}’}^{\!\top} \end{equation} $$

RBFの $k_g$ を入れて各成分を書き下すと、$\bm{r}=\bm{x}-\bm{x}’$ として

$$ K_{11}=\Big(\tfrac{1}{\ell^2}-\tfrac{r_2^2}{\ell^4}\Big)k_g,\quad K_{22}=\Big(\tfrac{1}{\ell^2}-\tfrac{r_1^2}{\ell^4}\Big)k_g,\quad K_{12}=K_{21}=\tfrac{r_1 r_2}{\ell^4}\,k_g $$

になります。あとはこの $K_{\bm{f}}$ を使って普通にGP回帰するだけ。カーネルを差し替えただけで、出てくる予測場は構造的に発散ゼロです。実装して確かめましょう。真の場として回転場 $\bm{f}=(-x_2,\,x_1)$(発散ゼロ)を14点だけ観測し、(a) 制約GP と (b) 成分独立の素朴GP を比べます。

import numpy as np

def kf_block(xi, xj, l, s=1.0):       # 発散ゼロカーネルの 2x2 ブロック
    r = xi - xj
    k = s**2 * np.exp(-(r @ r) / (2*l**2))
    K = np.empty((2, 2))
    K[0,0] = (1/l**2 - r[1]**2/l**4) * k
    K[1,1] = (1/l**2 - r[0]**2/l**4) * k
    K[0,1] = K[1,0] = (r[0]*r[1]/l**4) * k
    return K

def build_K(A, B, l, s=1.0):
    M = np.zeros((2*len(A), 2*len(B)))
    for i, xi in enumerate(A):
        for j, xj in enumerate(B):
            M[2*i:2*i+2, 2*j:2*j+2] = kf_block(xi, xj, l, s)
    return M

rng = np.random.default_rng(3)
l, s, noise = 1.3, 1.2, 1e-3
true_f = lambda P: np.stack([-P[:,1], P[:,0]], axis=1)   # 発散ゼロの回転場
Xtr = rng.uniform(-2, 2, size=(14, 2))
Ytr = true_f(Xtr) + 0.03*rng.standard_normal((14, 2))

K = build_K(Xtr, Xtr, l, s) + noise*np.eye(2*len(Xtr))
alpha = np.linalg.solve(K, Ytr.ravel())
mean_at = lambda P: (build_K(P, Xtr, l, s) @ alpha).reshape(-1, 2)   # 予測平均場

# 予測場の発散を中心差分で測る
def divergence(fn, P, h=1e-3):
    e1, e2 = np.array([[h,0]]), np.array([[0,h]])
    d1 = (fn(P+e1)[:,0] - fn(P-e1)[:,0]) / (2*h)
    d2 = (fn(P+e2)[:,1] - fn(P-e2)[:,1]) / (2*h)
    return d1 + d2

gx = np.linspace(-2.2, 2.2, 14)
grid = np.stack(np.meshgrid(gx, gx), -1).reshape(-1, 2)
print(f"制約GPの最大|発散| = {np.abs(divergence(mean_at, grid)).max():.2e}")
# -> 制約GPの最大|発散| = 1.63e-07  (素朴GPは 4.31e-01)

予測場の発散を格子上で測ると、制約GPは $1.6\times10^{-7}$(中心差分の打ち切り誤差レベル=事実上ゼロ)。一方、成分を独立に回帰した素朴GPは $0.43$ にもなります。実に270万倍の差です。図で見ると一目瞭然です。

真の非圧縮(発散ゼロ)場と14点の観測をベクトル矢印で表した図

これが真の場(灰)と14点の観測(赤)。この少ないデータから、2つのGPがどんな場を復元するかを見ます。

素朴GPの予測場:発散が大きく非物理的な湧き出し・吸い込みが各所に現れる

成分独立の素朴GP。色は発散の大きさで、赤は「湧き出し」、青は「吸い込み」を表します。観測の隙間で非物理的な発散があちこちに出て、回転場の構造も崩れています。

発散ゼロ制約GPの予測場:発散はほぼゼロ(白一色)で滑らかな回転場を正しく復元

発散ゼロ制約GP。色がほぼ真っ白=全域で発散がゼロ。同じ14点から、ちゃんと非圧縮な回転場を復元できています。事前知識(質量保存)を構造に埋め込んだ効果は絶大です。

素朴GPと制約GPの最大発散を対数軸の棒グラフで比較:制約GPは数値誤差レベルまで消える

最大発散を対数軸で並べた図。制約GPの発散は数値誤差の床まで落ちています。なぜここまで効くのか、理屈を一段だけ深掘りしておきましょう。

なぜ「全域で厳密」なのか

予測平均は $\bm{m}(\bm{x}_*)=\sum_j K_{\bm{f}}(\bm{x}_*,\bm{x}_j)\,\bm{\alpha}_j$ という、カーネル列の線形結合です。各列 $K_{\bm{f}}(\bm{x}_*,\bm{x}_j)=G_{\bm{x}_*}k_g\,G_{\bm{x}_j}^\top$ は、$\bm{x}_*$ について $G_{\bm{x}_*}$ の像にいます。だから発散を取ると $F_{\bm{x}_*}\bm{m}=\sum_j (F_{\bm{x}_*}G_{\bm{x}_*})k_g\,G^\top\bm{\alpha}_j=0$。条件(3) $FG=0$ がそのまま効いて、観測点だけでなく任意の $\bm{x}_*$ で、平均もサンプルも厳密にゼロ発散になります。「データで近づける」のではなく「そもそも制約を破る関数が確率1で存在しない」——これが線形制約の強さです。

ここがバニラGPとの決定的な違いです。素朴なGPは「観測点では制約っぽいデータに合わせる」ことしかできず、観測の隙間や外挿領域では制約が崩れます(さきほどの発散の図がまさにそれ)。対して線形制約GPは、仮説空間そのものを「制約を満たす関数の集合」に絞り込んでいる。サンプリングしても、事後平均を取っても、どこを切り取っても制約の外には出られません。データ量とも無関係で、観測が0点でも事前分布の段階ですでに発散ゼロです。

同じ枠組みで作れる、ほかの線形制約

この $\bm{f}=Gg,\ FG=0$ という設計は、発散ゼロ専用ではありません。$F$(守りたい線形制約)を決めて、その零空間を張る $G$ を見つければ、いろいろな物理が同じレシピで入ります。

  • 回転ゼロ(curl-free)の場:静電場や保存力場のように $\nabla\times\bm{f}=0$ を満たす場は、スカラーポテンシャルの勾配 $\bm{f}=\nabla g$ として書けます。つまり $G_{\bm{x}}=[\partial/\partial x_1,\ \partial/\partial x_2]^\top$。発散ゼロのときの「回転を取る」を「勾配を取る」に替えただけで、同じ構造です。$F=[\partial/\partial x_2,\ -\partial/\partial x_1]$(2次元の回転)に対して $FG=0$ が偏微分の順序交換で成り立ちます。
  • 線形PDEの解:$\mathcal{D}f=0$ となる線形微分作用素 $\mathcal{D}$(熱方程式・波動方程式・マクスウェル方程式など)も、$F=\mathcal{D}$ と置けば同じ枠組みに乗ります。$G$ を解析的に求めるのが難しい場合は、$F$ の係数から $G$ の係数を零空間の計算で機械的に求める手続きが Jidling らの論文で与えられています。
  • 境界条件:「領域の縁で値がゼロ」のような線形の境界条件も、その条件を満たす基底でカーネルを構成すれば厳密に組み込めます。

共通する発想は一つ。「制約を満たす関数だけで張られた空間」を先に用意し、その中でGPを定義する。データに頼らず構造で守る、という思想です。

ここまでが図の左側(厳密に守れる線形制約)。では右側、単調性や正値性のような不等式制約は、なぜ同じようにいかないのでしょうか。

不等式・形状制約:なぜ近似が必要なのか

「発散ゼロ」は $F\bm{f}=0$ という等式で、しかも $F$ が線形でした。だから零空間という線形代数の道具が使え、ガウス性も保たれた。

ところが「単調性 $f'(x)\ge0$」や「正値性 $f(x)\ge0$」は不等式です。不等式は「ある半空間に制限する」操作で、これは線形変換ではありません。ガウス分布を半分に切り落とすと、もう正規分布ではなく切断正規分布になります。事後分布がガウスでなくなる以上、これまでのような閉形式の条件付けは使えず、近似推論(EP・変分・MCMC)が必要になります。これが「②は近似」の理由です。

もう少し直感的に言うと、線形等式制約は「平面(部分空間)に乗れ」という拘束でした。ガウス分布を部分空間に射影してもガウスのまま——だから式一発で解けた。一方、不等式は「この半分の側にいろ」という拘束です。ベルのような正規分布の山を真ん中でスパッと切ると、断面はもうベル型ではない。平均も分散もずれ、解析的な公式が閉じない。さらに困るのは、制約点を増やすほど「複数の半空間の共通部分」に確率を押し込めることになり、その体積(規格化定数)を解析的に計算できなくなる点です。だからサンプリングや近似が要る。この「等式=部分空間、不等式=半空間」という幾何の違いが、厳密さの差を生む根っこです。

まず、制約なしGPで何が困るかを見ます。単調増加が分かっている関数(ロジスティック曲線)を、ノイズの乗った6点だけから回帰します。

制約なしGP:少データだと事後平均が非単調に波打ち、両端で減少してしまう

赤い帯が「事後平均が減少している区間」です。データが少なくノイズもあるせいで、増加するはずの関数なのに、あちこちで減少しています。両端の外挿では大きく下がってしまい、単調性という事前知識が完全に無視されています。

これを直す定番が、Riihimäki と Vehtari の仮想微分観測(AISTATS 2010)です。考え方はシンプルで、「ここで傾きは正であってほしい」という仮想点を置き、そこで $f'(x_v)>0$ という情報を尤度として与えます。ただし「正」という不等式をそのまま扱えないので、probit尤度で「傾きが正である確率」を表します。

$$ p\big(\text{単調}\mid f'(x_v)\big)=\Phi\!\Big(\frac{f'(x_v)}{\nu}\Big) $$

$\Phi$ は標準正規分布の累積分布関数で、$\nu$ は緩さを決める小さな定数です。

probit尤度の概念:微分が正なら強く支持、負なら罰せられる滑らかな尤度

この図のように、$f’>0$ なら尤度が1に近づき(強く支持)、$f’<0$ なら0に近づく(罰せられる)。微分の符号にだけ反応する、滑らかな「単調性の審判」です。この尤度はガウスでないので、厳密な事後分布は手に入らず、EP などで近似します。

ここでは仕組みを腹落ちさせるため、仮想点で「正の傾き」を緩く観測するという、線形観測による近似版を実装します(厳密な符号制約ではなく、正の目標傾きをゆるい誤差で与えることで単調性を促す)。微分観測の枠組みをそのまま再利用できます。

import numpy as np
# rbf, dk_dy, dk_dx, d2k は前の節と同じ定義を使う

true_f = lambda x: 1.0/(1.0 + np.exp(-2.5*x))      # 単調増加
l, s, noise = 1.1, 1.0, 0.05**2
Xtr = np.array([-2.5, -1.2, -0.2, 0.3, 1.4, 2.6])
ytr = true_f(Xtr) + np.array([0.06,-0.10,0.10,-0.09,0.07,-0.05])
Xs  = np.linspace(-3.2, 3.2, 200)

# 制約なし
K = rbf(Xtr, Xtr, l, s) + noise*np.eye(len(Xtr))
m0 = rbf(Xs, Xtr, l, s) @ np.linalg.solve(K, ytr)

# 仮想微分点で「傾き>0」を促す(正の微分を緩く観測)
Xv, dslope, dnoise = np.linspace(-3, 3, 21), 0.15, 0.08**2
Kff, Kfd = rbf(Xtr, Xtr, l, s), dk_dy(Xtr, Xv, l, s)
Kdf, Kdd = dk_dx(Xv, Xtr, l, s), d2k(Xv, Xv, l, s)
Kj = np.block([[Kff + noise*np.eye(len(Xtr)), Kfd],
               [Kdf, Kdd + dnoise*np.eye(len(Xv))]])
yj = np.concatenate([ytr, np.full(len(Xv), dslope)])   # 仮想点で傾き=正
Ks = np.hstack([rbf(Xs, Xtr, l, s), dk_dy(Xs, Xv, l, s)])
m1 = Ks @ np.linalg.solve(Kj, yj)

print("制約なしは単調か:", bool(np.all(np.diff(m0) >= -1e-6)))   # -> False
print("単調化後は単調か:", bool(np.all(np.diff(m1) >= -1e-6)))   # -> True

仮想微分点で単調化したGP:事後平均が全域で単調増加になり真の関数に沿う

仮想微分点(緑の三角)で「傾きは正」と教えた結果、事後平均は全域で単調増加になりました(np.diff がすべて非負)。両端の外挿も素直に増加しています。事前知識が予測を救った好例です。

ただし注意してほしいのは、これは「点ごとに緩く促した」結果で、線形制約のような全域・厳密の保証ではないこと。仮想点の間ではわずかに崩れうるし、厳密にやるなら probit 尤度+EP が要ります。この「厳密に守れる線形制約 vs 近似で促す不等式制約」という非対称性こそ、制約付きGPの肝です。

さらに踏み込むには(より厳密な不等式制約)

不等式を「全域で厳密に」守らせたい場合の代表が、Maatouk と Bay(2017)の有限次元近似です。GPを区分線形などの基底関数で展開し、その係数を切断ガウス分布から引きます。係数に課す不等式(係数が単調増なら関数も単調増、など)が、基底の性質を通じて領域全体の不等式制約に翻訳される仕組みです。これなら全条件付きサンプルが領域全域で制約を満たします。代償として、推論は切断ガウスからのサンプリング(MCMC等)になります。

正値性 $f>0$ や上下限のような「範囲」の制約には、もっと手軽なwarping(変数変換)もよく使われます。たとえば「必ず正」なら、観測対象そのものではなく $g=\log f$ をGPでモデル化し、予測を $f=\exp(g)$ と戻すだけ。指数関数の出力は必ず正なので、予測も信頼区間も負になりません(対数正規過程)。$[a,b]$ に収めたいならロジット変換を噛ませます。手軽な反面、warpingは出力の分布形を歪める(対称な信頼区間が非対称になる)ので、不確かさの解釈には少し注意が要ります。「厳密さ・全域性・実装の手軽さ」のどれを取るかで、手法を選ぶことになります。

制約の型 代表的な入れ方 守られ方
微分・勾配の観測 カーネルの微分を観測に追加 厳密
発散ゼロ・回転ゼロ・線形PDE $\bm{f}=G g,\ FG=0$ で構造的に 厳密・全域
境界条件(線形) 境界で値0などをカーネルに反映 厳密
単調性・凸性 仮想微分点+probit尤度+EP 近似・点ごと
正値・上下限 warping(対数GP等)/有限次元+切断ガウス 近似 or 全域(手法による)

まとめ

「GPに制約は入れられるのか?」への答えは、入れられる。ただし制約が線形(等式)か不等式かで、効かせ方と厳密さがまるで違う、でした。

  • GPは線形作用素で閉じている。だから微分・積分・線形結合を通してもGPのままで、共分散は $\mathcal{L}_x\mathcal{L}_{x’}k$ と閉形式で出る。
  • 線形(等式)制約は、$\bm{f}=Gg$ と「制約を満たす関数だけ」で場を張り直す($FG=0$)ことで、サンプルも予測も全域で厳密に満たせる。発散ゼロGPはカーネルを差し替えるだけで実現でき、予測場の発散は数値誤差レベルまで消えた。
  • 不等式・形状制約(単調・凸・正値)は半空間への制限でガウス性を壊すため、仮想点+probit尤度などで近似的に促す。全域厳密にしたいなら有限次元+切断ガウスへ。

物理法則を尊重する回帰は、データが少ない領域での外挿や、安全性が問われるサロゲートモデルで効きます。「事前知識をモデルの構造に埋め込む」というこの発想は、物理情報ニューラルネット(PINN)とも地続きです。

次のステップとして、以下もどうぞ。

参考文献:Jidling, Wahlström, Wills, Schön, “Linearly constrained Gaussian processes”, NeurIPS 2017 / Riihimäki, Vehtari, “Gaussian processes with monotonicity information”, AISTATS 2010 / Swiler, Gulian, Frankel, Safta, Jakeman, “A Survey of Constrained Gaussian Process Regression”, 2020 / Maatouk, Bay, “Gaussian process emulators for computer experiments with inequality constraints”, Math. Geosci. 2017。

ガウス過程回帰の理論と実装をわかりやすく解説
ガウス過程回帰の基礎。事前分布・事後分布・カーネルから予測の不確かさまでを導出と実装で解説します。
画像なし
ベイズ最適化の理論と実装 — ガウス過程で効率的に最適解を探索する
ガウス過程を使って少ない試行で最適解を探すベイズ最適化を、獲得関数の設計から実装まで解説します。