ベイズ最適化のお気持ち

TLDL;

ベイズ最適化をどう数理的とコードをつなげるかをNumPyで作成をしたコードを掲示をします。

用語説明

• ベイズ最適化では、ベイズ推論をしているNoise等のハイパラを最適化をする事を指します
• ベイズ推論とは、不確実性を含んだ予測モデルです

具体的な例としてはアンケートの最適化を考えてみましょう。
アンケート項目をどうするかを考える際に、愚直に集めた解答から次の項目を考えてみると、例えば適当にアンケートを”っぽい”と解答をした人がいたとすると、この結果が明らかに悪影響が出ると思います。これをどうにかするのが、ベイズ推論の具体例です。

問題設定

では、問題設定を考えてみましょう！
noiseは$\sigma \sim \mathcal{N}(0, 0.01)$として、設定をしています。

$$
y = \sin(2\pi{x}) + \sigma \cdot \sqrt{0.01}
$$
これでNoiseを含んだデータになっています。
今回はこのNoiseを推定をして、どれくらいの範囲にデータがあるか等を見ていくのが目的です！

まずは観測データを生成をしてみましょう。

np.random.seed(42) # 再現性確保の為に、Seed値固定
train_x = np.linspace(0,1,100) # 0から1の間で100点選択。
# 真の関数：sin(2*pi*x) に少しのノイズを加える
# np.random.randnは標準正規分布(平均0, 分散1)の値から選択をするので、$\sqrt{0.01}$を掛ける事で、分散を0.1の分布を作成できます
train_y = np.sin(train_x * (2 * np.pi)) + np.random.randn(len(train_x)) * np.sqrt(0.01) 
true_y = np.sin(train_x * (2 * np.pi)) 
# 次の点を探すための探索空間（候補点200個）
test_x = np.linspace(0, 1, 200)

では可視化をしてみましょう！

# ベイズ最適化のMotivationはこのノイズの含んだデータでも予測をしたい。という事です。
plt.title('ノイズがどう作用をしているかを理解する。')
plt.plot(train_x, train_y, label="実験データのPlot")
plt.plot(train_x, true_y, label="本来の理想のデータ")
plt.legend()
plt.show()

出力結果が以下です。
以下の出力結果を見るとわかる通り、やはり、Noiseが入っているので、本来の関数からはギザギザでずれています。
このギザギザを推定をしてみようがモチベーションです。

ベイズ最適化

カーネル関数

ベイズと切っても切り離せないのが、カーネル関数です。
これは非常に便利ですが、カーネルトリック等を理解をすると良いとは思いますが、恐らくこの記事では書ききれないので、後日記載をした記事のリンクを挿入をしておきます。
今回は、基本的なRBFカーネルを使います。
$$
RBF(x_1, x_2) = \sigma^2 exp(\frac{|x_1 - x_2|^2}{2l^2})
$$
上記で表せます。
コードでは以下と記載をします。

def rbf_kernel(x1, x2, output_scale, lengthsale):
    sq_dist = (x1[:, None], x2[None, :]) 
    return output_scale * np.exp(sq_dist **2 / 2 * lengthscale **2)

獲得関数

次に獲得関数について説明をします。
獲得関数のイメージとしては、$x_1$という点を選択をした際に、情報量が増えるように$x_2$を選択をする。その選択をする記述です。
今回はUCBを使用をしてみます。
$\mu$は平均、$\sigma$は分散をあらわし、$\beta$は、こちらもハイパラです。
$$
UCB(\x) = \mu + \beta \sigma
$$

コードでは

def upper_confidence_bound(mu, std, beta=2.0)
    return mu + beta * std

周辺尤度最大化

まず周辺尤度の式を記載をします。導出はLLMに聞いてみる事をお勧めします。
線形代数、統計学と色々な知識が必要で楽しいです。$\theta$はハイパラメータ(Noise, Lengthscale, outputscale)の事を指し示します。
$$
log_p(y | X, \theta) = -\frac{1}{2}y^T\Sigma^{-1}y - \frac{1}{2}log|\Sigma| - \frac{N}{2}log(2 \pi)
$$

で表せます。実装は以下です。

def negative_marginal_log_likelihood(paramas, x, y):
    lengthscale, outputscale, noise = params
    N = len(y)
    K_noise = rbf_kernel(x, y, outputscale, lengthscale) + noise * np.eye(N) # eyeは単位行列を作成します。
    try: 
        L = np.linalg.cholesky(K_noise) # 正定値行列であるかの保証を確認をしていないので、念の為tryでエラーを回避します。
    except np.linalg.LinAlgError:
        return 1e10 # 失敗をした際に大き目の値で更新をやり直させます。
    alpha = np.linalg.solve(L.T, np.linalg.solve(L, y))
    logdet = 2 * np.sum(np.log(np.diag)) # 行列式計算
    return 0.5 * y @ alpha + 0.5 * logdet + 0.5 * N * np.log(2 * np.pi)

事後分散・平均の算出

ベイズ推論では、事後分散、事後平均を求める事で、予測の不確実性を明らかにします。
まず以下で事後分散・平均が算出をされます。
事後平均

$$
\mu_{* } = K_{X_{*}X}(K_{XX} + \sigma_{n}^2 I)^{-1}y
$$

mu_s = K_sX @ np.linalg.inv(K_XX + noise * np.eye(N)) @ train_y

事後分散

$$
\Sigma_{* } = K_{x_{* }X}(K_{XX} + \sigma_{n}^2I)^{-1}K_{XX}
$$

sigma_s = K_ss - K_sX @ np.linalg.inv(K_XX + noise * np.eye(N)) @ K_sX.T

予測全体のコードを見ます。

def predict(test_x, train_x, train_y, lengthscale, output_scale, noise):
    N_train = len(train_x)
    K_XX = rbf_kernel(train_x, train_x, outputscale, lengthscale)
    K_sX = rbf_kernel(test_x, train_x, outputscale, lengthscale)
    K_ss = rbf_kernel(test_x, test_x, outputscale, lengthscale)
    L = np.linalg.solve(L.T, np.linalg.solve(L, train_y))
    mu_s = K_sX @ alpha
    v = np.linalg.solve(L, K_sX.T)
    sigma_s = K_ss - v.T @ v
    pred_var = np.sqrt(np.diag(sigma_s), 1e-8)
    return mu_s, np.sqrt(pred_var)

では最適化と予測、訓練を実施をします。
コードはAppendexにて掲示をしますが、ここでは先に訓練後の図を掲示をします。
以下を見ればある程度分かりますが、上では訓練をした後のガウス過程モデルで予測をした図です。
このように、青色で不確実性を予測をできており、観測点に関しても、今回は運よく、すべて不確実性で予測をされている部分に入っています。
上記の通り簡単にですが、説明をしました。
理論的な部分は恐らく今後勉強の為に書いていく予定ですが、Appendexにてコードをすべて掲示をしているので、暇があれば動かしてみてください。

Appendex

今回使用をしたすべてのコード

# %% [markdown]
# ## import

# %%
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = "IPAGothic"

# %% [markdown]
# ## Generate data

# %%
# ==========================================
# 1. 初期データの設定（あえて少ない6点からスタート）
# ==========================================
np.random.seed(42)
train_x = np.linspace(0,1,10)
# 真の関数：sin(2*pi*x) に少しのノイズを加える
train_y = np.sin(train_x * (2 * np.pi)) + np.random.randn(len(train_x)) * np.sqrt(0.01)
true_y = np.sin(train_x * (2 * np.pi)) 
# 次の点を探すための探索空間（候補点200個）
test_x = np.linspace(0, 1, 200)

# %%
# ベイズ最適化のMotivationはこのノイズの含んだデータでも予測をしたい。という事です。
plt.title('ノイズがどう作用をしているかを理解する。')
plt.plot(train_x, train_y, label="実験データのPlot")
plt.plot(train_x, true_y, label="本来の理想のデータ")
plt.legend()
plt.show()

# %% [markdown]
# ## Gaussian Regression

# %%
# ==========================================
# 2. ガウス過程の基本関数の定義
# ==========================================
def rbf_kernel(x1, x2, lengthscale, output_scale):
    """RBF（高斯）カーネルの計算"""
    sq_dist = (x1[:, None] - x2[None, :]) ** 2 # (N, 1) - (1, N)で計算をしている。これって正しいの？->for ループを使わないで処理をしようとするという動機で合致をする。
    return output_scale * np.exp(-sq_dist / (2 * lengthscale ** 2))


def negative_marginal_log_likelihood(params, x, y):
    """paramsは生の値 (lengthscale, output_scale, noise)"""
    lengthscale, output_scale, noise = params
    N = len(y)
    # 共分散行列の左上
    K_noise = rbf_kernel(x, x, lengthscale, output_scale) + noise * np.eye(N) # $K_{XX} + noise * I$
    try:
        L = np.linalg.cholesky(K_noise)          # 逆行列でなくCholeskyを使う
    except np.linalg.LinAlgError:
        return 1e10
    # L^{-1} * yを算出
    alpha = np.linalg.solve(L.T, np.linalg.solve(L, y))
    
    logdet = 2 * np.sum(np.log(np.diag(L)))
    return 0.5 * y @ alpha + 0.5 * logdet + 0.5 * N * np.log(2 * np.pi)


def predict(test_x, train_x, train_y, lengthscale, output_scale, noise):
    # 行数を解答、ここで行数はBatch数をあらわしている。ただ、行列にし、うまく計算をする事で、本来\Order{N}のものが\Order{N^{1/2}}等になりうれしい。
    N_train = len(train_x) 
    K_XX = rbf_kernel(train_x, train_x, lengthscale, output_scale) + noise * np.eye(N_train)
    K_sX = rbf_kernel(test_x, train_x, lengthscale, output_scale)
    K_ss = rbf_kernel(test_x, test_x, lengthscale, output_scale)
    L = np.linalg.cholesky(K_XX)
    alpha = np.linalg.solve(L.T, np.linalg.solve(L, train_y))
    mu_s = K_sX @ alpha
    v = np.linalg.solve(L, K_sX.T)
    sigma_s = K_ss - v.T @ v
    pred_var = np.maximum(np.diag(sigma_s), 1e-8)
    return mu_s, np.sqrt(pred_var)

# %% [markdown]
# ## Aquisition Function

# %%
def upper_confidence_bound(mu, std, beta=2.0):
    return mu + beta * std

# %% [markdown]
# ## Marginal Log Likelihood

# %%
def numpy_minimize(train_x, train_y, init_params, epochs=10, lr=0.05, h=1e-5):
    """log空間で最適化。正値制約は exp() により自動的に満たされる。"""
    log_theta = np.log(np.array(init_params, dtype=float)) # hyperparameters
    for epoch in range(epochs):
        raw = np.exp(log_theta)
        current_loss = negative_marginal_log_likelihood(raw, train_x, train_y)
        gradients = np.zeros_like(log_theta)
        for i in range(len(log_theta)):
            pert = log_theta.copy(); pert[i] += h
            perturbed_loss = negative_marginal_log_likelihood(np.exp(pert), train_x, train_y)
            gradients[i] = (perturbed_loss - current_loss) / h
        log_theta -= lr * gradients
        if (epoch + 1) % 10 == 0 or epoch == 0:
            print(f"Epoch {epoch+1:3d}/{epochs} - Loss: {current_loss:.3f} | "
                  f"lengthscale: {raw[0]:.3f} | output_scale: {raw[1]:.3f} | noise: {raw[2]:.4f}")
    return np.exp(log_theta)

# %%
def numpy_minimize_with_bounds(train_x, train_y, init_params, epochs=10, lr=0.05, h=1e-5):
    log_theta = np.log(np.array(init_params), dtype=float)
    gradients = np.zeros_like(log_theta)
    
    # --- 1. 実数空間での Bounds を定義 ---
    # (例: lengthscale, output_scale, noise の探索範囲をそれぞれ [1e-3, 10.0] に制限)
    lower_bounds = np.array([0.3, 0.9, 1e-4])
    upper_bounds = np.array([0.4, 1.0, 1.0])
    # raw
    # --- 2. 対数空間用の Bounds に変換 ---
    log_lower = np.log(lower_bounds)
    log_upper = np.log(upper_bounds)
    for epoch in range(epochs):
        current_loss = negative_marginal_log_likelihood(np.exp(log_theta), train_x, train_y)        
        for i in range(len(log_theta)):
            pert = log_theta.copy()
            pert[i] += h 
            perturbed_loss = negative_marginal_log_likelihood(np.exp(pert), train_x, train_y)
            gradients[i] = (perturbed_loss - current_loss) / h
            
        # パラメータの更新
        log_theta -= lr * gradients
        # --- 3. 【ここを追加】Bounds の範囲内に丸め込む（投影） ---
        log_theta = np.clip(log_theta, log_lower, log_upper)
        # if (epoch + 1) % 30 == 0 or epoch == 0:
        #     print(f"Epoch {epoch+1:3d}/{epochs} - Loss: {current_loss:.3f} | "
        #           f"lengthscale: {raw[0]:.3f} | output_scale: {raw[1]:.3f} | noise: {raw[2]:.4f}")
    return np.exp(log_theta)

# %% [markdown]
# ## training & prediction & aqisition func calc

# %%
# ==========================================
# 4. メイン処理：モデル訓練 ➔ 予測 ➔ 獲得関数で次の一手を選択
# ==========================================

# ── Step 1: モデルの訓練 (MLLの最小化) ──
initial_params = [0.2, 1.0, 0.05]
bounds = [(1e-3, None), (1e-3, None), (1e-6, None)]
opt_params = numpy_minimize_with_bounds(train_x, train_y, initial_params, epochs=10, lr=0.01)
opt_lengthscale, opt_output_scale, opt_noise = opt_params

# ── Step 2: 予測の計算 ──
pred_mean, pred_std = predict(test_x, train_x, train_y, opt_lengthscale, opt_output_scale, opt_noise)

# ── Step 3: 獲得関数の計算と「次の点」の選定 ──
# 探索空間全体（test_x）に対してUCBスコアを計算
ucb_scores = upper_confidence_bound(pred_mean, pred_std, beta=2.5)
# print('UCB Score')
# print(ucb_scores)
# print('-'*40)

# UCBスコアが「最も高いインデックス」を取得 ＝ 次に実験すべき場所
next_index = np.argmax(ucb_scores)
next_x = test_x[next_index]
next_y_score = ucb_scores[next_index]

print(f"候補となる次の実験点 x: {next_x:.4f} (UCBスコア: {next_y_score:.3f})")

# UCBも分布になるはず。
# plt.plot(ucb_scores)
# plt.show()

# %%
# ==========================================
# 5. 可視化（予測分布と獲得関数を上下に並べる）
# ==========================================
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(7, 6), sharex=True)

# 上段：ガウス過程の予測結果
ax1.plot(train_x, train_y, 'k*', markersize=10, label='Observed Data')
ax1.plot(test_x, np.sin(test_x * (2 * np.pi)), 'r--', alpha=0.5, label='True Function')
ax1.plot(test_x, pred_mean, 'b-', label='GP Mean')
ax1.fill_between(test_x, pred_mean - 2 * pred_std, pred_mean + 2 * pred_std, alpha=0.2, color='blue', label='Confidence ($2\sigma$)')
ax1.set_title("1. Gaussian Process Regression")
ax1.set_ylabel("y")
ax1.legend(loc="lower left")
ax1.grid(True)

# 下段：獲得関数（UCB）の形状と、次の一手

ax2.plot(test_x, ucb_scores, 'g-', label='UCB Score')
ax2.axvline(x=next_x, color='red', linestyle=':', linewidth=2, label=f'Next Sampling Point (x={next_x:.2f})')
ax2.plot(next_x, next_y_score, 'ro')  # 最大値の点に赤丸
ax2.set_title("2. Acquisition Function (UCB)")
ax2.set_xlabel("x")
ax2.set_ylabel("Acquisition Value")
ax2.legend(loc="lower left")
ax2.grid(True)

plt.tight_layout()
plt.show()

# %%
!pip list | grep numpy

# %%