botorch入門2 モデル作成2

前回に引き続きモデル作成について記述します。

https://qiita.com/tanaka_benkyo/items/de21d421ebad17bb0b78

今回は多目的モデルについてまとめていきます。

致命的なミスなどあればご指摘いただければと思います。

## 5. MultiOutput(SingleTaskGP)

SingleTaskGP は本来1次元出力のモデルですが、train_Y に (N, M) 形状の多次元出力を与えることで、1つのモデル内で多目的出力（M目的）を同時に学習することも可能です。この場合、BoTorch は出力次元を別々の目的として処理します。
ただし出力間の相関は明示的にモデリングされないため、独立な多目的最適化と解釈されます。

次のような関数からデータ点をサンプリングしてモデルの作成を行います。
負の相関を持つような2つの目的変数を考えます。


```python
# # トイデータの定義（2変数入力・目的の大小が反転）
# def blackbox1(x):
#     return np.sin(x[:, 0:1]*2) + np.cos(x[:, 0:1]/2 - 2) - x[:, 0:1]*0.1 + 2*np.exp(-np.abs(x[:, 1:2]-5)) + np.random.rand(len(x), 1)

# def blackbox2(x):
#     return -0.9*np.sin(x[:, 0:1]*2) - 1.1*np.cos(x[:, 0:1]/2 - 2) - x[:, 0:1]*0.1 + 2*np.exp(-np.abs(x[:, 1:2]-5)) + np.random.rand(len(x), 1)


def blackbox1(x):
    return np.sin(x[:, 0:1]*2) + np.cos(x[:, 0:1]/2 - 2) - x[:, 0:1]*0.1 + np.random.rand(len(x), 1)

def blackbox2(x):
    return -0.8*np.sin(x[:, 0:1]*2) - 1.2*np.cos(x[:, 0:1]/2 - 2) - x[:, 0:1]*0.15 + np.random.rand(len(x), 1)

x = np.arange(0,10,0.1).reshape(-1,1)
y1 = blackbox1(x)
y2 = blackbox2(x)
plt.plot(y1)
plt.plot(y2)
```


![output_41_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/3808cf7d-27c6-4321-a4be-740a9d925219.png)


```python
train_x = torch.rand(10,1, dtype=torch.double)*10
train_y = torch.cat([blackbox1(train_x), blackbox2(train_x)], dim=1)
```


```python
bounds = torch.stack([
    torch.min(train_x, dim=0)[0],
    torch.max(train_x, dim=0)[0]
])

# 入力データを正規化
train_x_scaled = normalize(train_x, bounds=bounds)

# モデルの作成
model = SingleTaskGP(
    train_x_scaled,
    train_y,
    outcome_transform=Standardize(
        m=train_y.size(-1)
    )
)

# モデルの最尤推定を行う
mll = ExactMarginalLogLikelihood(
    model.likelihood,
    model
)
fit_gpytorch_mll(mll)
```







```python
x = torch.arange(0,10,0.1).reshape(-1,1)
x_norm = normalize(x, bounds)

pred = model.posterior(x_norm)
pred_mean = pred.mean.detach().numpy()
pred_std = pred.variance.sqrt().detach().numpy()

plt.plot(x.ravel(), pred_mean[:,0], color='orange')
plt.fill_between(x.ravel(), pred_mean[:,0]+pred_std[:,0], pred_mean[:,0]-pred_std[:,0], color='orange', alpha=.4)

plt.plot(x.ravel(), pred_mean[:,1], color='green')
plt.fill_between(x.ravel(), pred_mean[:,1]+pred_std[:,1], pred_mean[:,1]-pred_std[:,1], color='green', alpha=.4)

plt.scatter(train_x, train_y[:,0], color='red')
plt.scatter(train_x, train_y[:,1], color='purple')
plt.xlabel('x', size=20)
plt.ylabel('y', size=20)
```


![output_44_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/20a23c06-f3c2-42dd-b775-63434f725dc2.png)

    


## 6. ModelListGP

**概要**:
`ModelListGP` は複数の `SingleTaskGP` モデルなどを一括管理できるラッパーであり、出力が独立している場合に有効です。

**主な用途**:

* 多目的ベイズ最適化（それぞれの目的ごとに独立モデルを構築）
* 出力間の依存を考慮しないマルチ出力問題

**特徴**:

* 各出力に異なるモデル・カーネルを設定可能
* 多目的獲得関数（例：EHI、NEHI）と組み合わせやすい


次のような関数からデータ点をサンプリングしてモデルの作成を行います。
負の相関を持つような2つの目的変数を考えます。


```python
from botorch.models import ModelListGP
from botorch.models.transforms import Standardize
from gpytorch.mlls.sum_marginal_log_likelihood import SumMarginalLogLikelihood
```


```python
train_x = torch.rand(10,1, dtype=torch.double)*10
train_y1 = blackbox1(train_x)
train_y2 = blackbox2(train_x)
```


```python
bounds = torch.stack([
    torch.min(train_x, dim=0)[0],
    torch.max(train_x, dim=0)[0]
])

# 入力データを正規化
train_x_scaled = normalize(train_x, bounds=bounds)

# それぞれの目的に対する GP モデルを作成
model1 = SingleTaskGP(
    train_x_scaled,
    train_y1,
    outcome_transform=Standardize(m=1)
)
model2 = SingleTaskGP(
    train_x_scaled,
    train_y2,
    outcome_transform=Standardize(m=1)
)

# ModelListGP にまとめる（出力は独立と仮定）
model = ModelListGP(model1, model2)

# 複数モデルのロス関数は SumMarginalLogLikelihood を使用
mll = SumMarginalLogLikelihood(model.likelihood, model)
fit_gpytorch_mll(mll)
```





```python
x = torch.arange(0,10,0.1).reshape(-1,1)
x_norm = normalize(x, bounds)

pred = model.posterior(x_norm)
pred_mean = pred.mean.detach().numpy()
pred_std = pred.variance.sqrt().detach().numpy()

plt.plot(x.ravel(), pred_mean[:,0], color='orange')
plt.fill_between(x.ravel(), pred_mean[:,0]+pred_std[:,0], pred_mean[:,0]-pred_std[:,0], color='orange', alpha=.4)

plt.plot(x.ravel(), pred_mean[:,1], color='green')
plt.fill_between(x.ravel(), pred_mean[:,1]+pred_std[:,1], pred_mean[:,1]-pred_std[:,1], color='green', alpha=.4)

plt.scatter(train_x, train_y[:,0], color='red')
plt.scatter(train_x, train_y[:,1], color='purple')
plt.xlabel('x', size=20)
plt.ylabel('y', size=20)
```

![output_49_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/4f687184-3cf2-4547-be4e-66164703a65c.png)

    


## 7. MultitaskGP

**概要**:
`MultitaskGP` は、複数の出力（タスク）を同時にモデリングし、タスク間の共通構造を活用するモデルです。

**主な用途**:

* 条件付き実験（温度、圧力などのタスク間比較）
* 複数目的の相関があるケース

**特徴**:

* 内部的には `MultitaskKernel` を用いて出力間の共通パターンを捉える
* タスク間で情報共有できるため、少数データにも有効

次のような関数からデータ点をサンプリングしてモデルの作成を行います。
負の相関を持つような2つの目的変数を考えます。データ数に偏りがある場合を想定します。


```python
from botorch.models.multitask import MultiTaskGP
```


```python
train_x1 = torch.rand(20,1, dtype=torch.double)*10
train_y1 = blackbox1(train_x1)

train_x2 = torch.rand(5,1, dtype=torch.double)*10
train_y2 = blackbox2(train_x2)


train_x = torch.cat([
    torch.cat([train_x1, torch.zeros_like(train_x1)], dim=1),
    torch.cat([train_x2, torch.ones_like(train_x2)], dim=1),
])
train_y = torch.cat([train_y1, train_y2])
```


```python
bounds = torch.stack([
    torch.min(train_x, dim=0)[0],
    torch.max(train_x, dim=0)[0]
])

# 入力データを正規化
train_x_scaled = normalize(train_x, bounds=bounds)

model = MultiTaskGP(
    train_x_scaled,
    train_y,
    task_feature=-1,
    outcome_transform=Standardize(
        m=train_y.size(-1)
    )
)

mll = ExactMarginalLogLikelihood(model.likelihood, model)
fit_gpytorch_mll(mll)
```

   


```python
x = torch.arange(0,10,0.1).reshape(-1,1)
x_norm = normalize(x, bounds)[:,0:1]

x_norm = torch.cat([
    torch.cat([x_norm, torch.zeros_like(x_norm)], dim=1),
    torch.cat([x_norm, torch.ones_like(x_norm)], dim=1),
])

pred = model.posterior(x_norm)
pred_mean = pred.mean.detach().numpy().reshape(2,-1).T
pred_std = pred.variance.sqrt().detach().numpy().reshape(2,-1).T

plt.plot(x.ravel(), pred_mean[:,0], color='orange')
plt.fill_between(x.ravel(), pred_mean[:,0]+pred_std[:,0], pred_mean[:,0]-pred_std[:,0], color='orange', alpha=.4)

plt.plot(x.ravel(), pred_mean[:,1], color='green')
plt.fill_between(x.ravel(), pred_mean[:,1]+pred_std[:,1], pred_mean[:,1]-pred_std[:,1], color='green', alpha=.4)

plt.scatter(train_x1[:,0], train_y1, color='red')
plt.scatter(train_x2[:,0], train_y2, color='purple')
plt.xlabel('x', size=20)
plt.ylabel('y', size=20)
```


![output_54_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/0294572c-1b22-4eec-9095-ec66b1d28179.png)

    


## 8. KroneckerMultitaskGP

**概要**:
`KroneckerMultitaskGP` は、出力数が多いときでも計算効率を維持できるように設計されたモデルで、Kronecker 積を用いてカーネルを表現します。

**主な用途**:

* 多数の出力を持つマルチタスク問題
* タスク間構造を効率的に学習したい場合

**特徴**:

* モデルサイズの増加に対してスケーラブル
* 出力間のカーネル構造を明示的に分離・統合可能

次のような関数からデータ点をサンプリングしてモデルの作成を行います。
負の相関を持つような2つの目的変数を考えます。同じ条件で2つの目的変数が得られているt場合を考えます。


```python
from botorch.models.multitask import KroneckerMultiTaskGP
```


```python
train_x = torch.rand(10,1, dtype=torch.double)*10
train_y1 = blackbox1(train_x)
train_y2 = blackbox2(train_x)

train_y = torch.cat([train_y1, train_y2], dim=1)
```


```python
bounds = torch.stack([
    torch.min(train_x, dim=0)[0],
    torch.max(train_x, dim=0)[0]
])

# 入力データを正規化
train_x_scaled = normalize(train_x, bounds=bounds)

model = KroneckerMultiTaskGP(
    train_x_scaled,
    train_y,
    outcome_transform=Standardize(
        m=train_y.size(-1)
    )
)

mll = ExactMarginalLogLikelihood(model.likelihood, model)
fit_gpytorch_mll(mll)
```





```python
x = torch.arange(0,10,0.1).reshape(-1,1)
x_norm = normalize(x, bounds)

pred = model.posterior(x_norm)
pred_mean = pred.mean.detach().numpy()
pred_std = pred.variance.sqrt().detach().numpy()

plt.plot(x.ravel(), pred_mean[:,0], color='orange')
plt.fill_between(x.ravel(), pred_mean[:,0]+pred_std[:,0], pred_mean[:,0]-pred_std[:,0], color='orange', alpha=.4)

plt.plot(x.ravel(), pred_mean[:,1], color='green')
plt.fill_between(x.ravel(), pred_mean[:,1]+pred_std[:,1], pred_mean[:,1]-pred_std[:,1], color='green', alpha=.4)

plt.scatter(train_x, train_y[:,0], color='red')
plt.scatter(train_x, train_y[:,1], color='purple')
plt.xlabel('x', size=20)
plt.ylabel('y', size=20)
```


![output_59_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/4e2c42c1-25fe-4690-a7cf-e35fdf6d4049.png)



## 9. MultiFidelityGP

**概要**:
`MultiFidelityGP` は、異なる信頼性・コストの情報源（例：低精度なシミュレーション vs 高精度な実験）を同時に学習し活用するモデルです。

**主な用途**:

* マルチフィデリティベイズ最適化
* 精度・コストトレードオフのある設計問題

**特徴**:

* Fidelity index を入力に加えることで階層情報を表現
* 高コストな高精度データの節約が可能

次のような関数からデータ点をサンプリングしてモデルの作成を行います。
2つの目的変数間にバイアスがあるような場合を考えます。データ数に偏りがあり、多く得られる方がノイズが大きいものとします。


```python
from botorch.models.gp_regression_fidelity import SingleTaskMultiFidelityGP
```


```python
def blackbox_high(x):
    return np.sin(x[:, 0:1]*2) + np.cos(x[:, 0:1]/2 - 2) - x[:, 0:1]*0.1 + 0.5*np.random.rand(len(x), 1)

def blackbox_low(x):
    return np.sin(x[:, 0:1]*2) +np.cos(x[:, 0:1]/2 - 2) - x[:, 0:1]*0.15 + +1.5 + 2*np.random.rand(len(x), 1)

x = np.arange(0,10,0.1).reshape(-1,1)
y1 = blackbox_high(x)
y2 = blackbox_low(x)
plt.plot(y1)
plt.plot(y2)
```




![output_63_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/c82b51e8-c691-4dcc-a449-f3c2992e6a47.png)



```python
train_x1 = torch.rand(50,1, dtype=torch.double)*10
train_y1 = blackbox_low(train_x1)

train_x2 = torch.rand(5,1, dtype=torch.double)*10
train_y2 = blackbox_high(train_x2)


train_x = torch.cat([
    torch.cat([train_x1, torch.zeros_like(train_x1)], dim=1),
    torch.cat([train_x2, torch.ones_like(train_x2)], dim=1),
])
train_y = torch.cat([train_y1, train_y2])
```


```python
bounds = torch.stack([
    torch.min(train_x, dim=0)[0],
    torch.max(train_x, dim=0)[0]
])

# 入力データを正規化
train_x_scaled = normalize(train_x, bounds=bounds)

model = SingleTaskMultiFidelityGP(
    train_x_scaled,
    train_y,
    outcome_transform=Standardize(m=1),
    data_fidelities=[1]
)

mll = ExactMarginalLogLikelihood(
    model.likelihood,
    model
)
fit_gpytorch_mll(mll)
```





```python
x = torch.arange(0,10,0.1).reshape(-1,1)
x_norm = normalize(x, bounds)[:,0:1]

x_norm = torch.cat([
    torch.cat([x_norm, torch.zeros_like(x_norm)], dim=1),
    torch.cat([x_norm, torch.ones_like(x_norm)], dim=1),
])

pred = model.posterior(x_norm)
pred_mean = pred.mean.detach().numpy().reshape(2,-1).T
pred_std = pred.variance.sqrt().detach().numpy().reshape(2,-1).T

plt.plot(x.ravel(), pred_mean[:,0], color='orange')
plt.fill_between(x.ravel(), pred_mean[:,0]+pred_std[:,0], pred_mean[:,0]-pred_std[:,0], color='orange', alpha=.4)

plt.plot(x.ravel(), pred_mean[:,1], color='green')
plt.fill_between(x.ravel(), pred_mean[:,1]+pred_std[:,1], pred_mean[:,1]-pred_std[:,1], color='green', alpha=.4)

plt.scatter(train_x1[:,0], train_y1, color='red')
plt.scatter(train_x2[:,0], train_y2, color='purple')
plt.xlabel('x', size=20)
plt.ylabel('y', size=20)
```


![output_66_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/79a568f4-5675-45d7-ae6d-e4c9cb7f486d.png)

    


## 10. SaasFullyBayesianMultiTaskGP

**概要**:
`SaasFullyBayesianMultiTaskGP` は、マルチタスク構造と完全ベイズ推論を組み合わせた高性能モデルで、BoTorch で高度なベイズ最適化を行うために設計されています。

**主な用途**:

* タスク間構造を持つ高次元問題の不確実性推論

**特徴**:

* `Saas` モデルの利点（変数選択・事後分布）にタスク構造のモデリングを加える
* 学習は `fit_fully_bayesian_model_nuts` を通じて実行

次のような関数からデータ点をサンプリングしてモデルの作成を行います。
20次元のデータを用意しますが、目的変数は4つの説明変数のみに依存しているものを想定します。


```python
from botorch.models.fully_bayesian_multitask import SaasFullyBayesianMultiTaskGP
```


```python
f_hd1 = lambda y: 3*torch.sin(y[:,[4]]*6) + 2*torch.sin(y[:,[9]]*2+torch.pi/6) - torch.sin(y[:,[14]]*2+torch.pi/1.5) + 2.5*torch.sin(y[:,[19]]*4) + 0.2*torch.randn(len(y),1)
f_hd2 = lambda y: -2*torch.sin(y[:,[3]]*6) - 2.5*torch.sin(y[:,[8]]*2+torch.pi/6) + 0.5*torch.sin(y[:,[14]]*2+torch.pi/1.5) - 3*torch.sin(y[:,[19]]*4) + 0.3*torch.randn(len(y),1)

train_x1 = torch.rand(20,20, dtype=torch.double)
train_y1 = f_hd1(train_x1)

train_x2 = torch.rand(5,20, dtype=torch.double)
train_y2 = f_hd1(train_x2)

train_x = torch.cat([
    torch.cat([train_x1, torch.zeros_like(train_x1)], dim=1),
    torch.cat([train_x2, torch.ones_like(train_x2)], dim=1),
])
train_y = torch.cat([train_y1, train_y2])
```


```python
model = SaasFullyBayesianMultiTaskGP(
    train_x,
    train_y,
    task_feature=-1,
    outcome_transform=Standardize(m=train_Y.size(-1))
)

fit_fully_bayesian_model_nuts(
    model,
    warmup_steps=128,
    num_samples=64,
    thinning=16,
    disable_progbar=False
)
```

    Sample: 100%|█████████████████████████████████████████| 192/192 [04:50,  1.51s/it, step size=3.09e-01, acc. prob=0.860]
    


```python
plot_lengthscale_importance(model)
```

![output_71_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/acf47ad7-2620-46d7-b27b-1130108ce8fc.png)

    


## 11. HigherOrderGP

**概要**:
`HigherOrderGP` は、BoTorch の拡張で導入された高次構造（例: 時系列画像やテンソルなど）に対応したモデルです。観測出力がベクトルやテンソルである場合に有効です。

**主な用途**:

* 時系列ベクトル、2次元画像のような構造化出力の回帰
* GPによる動的システムの応答予測など

**特徴**:

* 出力がテンソル構造（例：画像、時系列）であってもそのまま扱える
* FlattenedStandardize によりテンソル標準化を適用
* 潜在空間における構造保持が可能

次のような関数からデータ点をサンプリングしてモデルの作成を行います。
一つの説明変数に対して、出力は画像である場合を考えます。説明変数が大きくなほど、画像の模様が大きな値を取るものとします。


```python
from botorch.models import HigherOrderGP
from botorch.models.higher_order_gp import FlattenedStandardize
from botorch.optim.fit import fit_gpytorch_mll_torch
from functools import partial
from torch.optim import Adam
from linear_operator.settings import _fast_solves
from botorch.models.transforms import Normalize, Standardize
```


```python
def make_img(t):
    x1, x2 = np.meshgrid(np.arange(10), np.arange(10))
    return np.sin(t/8) * (2 ** np.cos(x1*1.5) + 3**np.sin(x2*1) + np.exp(x1 / 10)) + np.random.rand(10,10)

t = 10
plt.imshow(make_img(t))
```



![output_74_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/cd7b986c-6278-498f-ac9d-ac0048224173.png)

    



```python
train_x=torch.randint(high=100, size=(10,1)).double()
train_y=torch.tensor([make_img(t.detach().numpy()) for t in train_x])

bounds = torch.stack([
    torch.min(train_x, dim=0)[0],
    torch.max(train_x, dim=0)[0]
])

# 入力データを正規化
train_x_scaled = normalize(train_x, bounds=bounds)

model = HigherOrderGP(
            train_x_scaled,
            train_y,
            outcome_transform=FlattenedStandardize(train_y.shape[1:]),
            input_transform=Normalize(train_x_scaled.shape[-1]),
            latent_init="gp",
        )

mll = ExactMarginalLogLikelihood(model.likelihood, model)

with _fast_solves(True):
    fit_gpytorch_mll_torch(
        mll, step_limit=200, optimizer=partial(Adam, lr=0.01)
    )
```


```python
train_x
```




    tensor([[81.],
            [63.],
            [39.],
            [86.],
            [36.],
            [10.],
            [44.],
            [23.],
            [90.],
            [81.]], dtype=torch.float64)



画像の予測結果と実際の画像の比較と、画素値のyyプロットを確認します。


```python
T = torch.randint(high=100, size=(10,1)).double()
T_scaled = normalize(T, bounds=bounds)

true_y=[make_img(t.detach().numpy()) for t in T]
pred_y = [model.posterior(torch.tensor([[t]])).mean.detach().numpy()[0] for t in T_scaled]
```


```python
for i in range(len(T)):
    fig, ax = plt.subplots(1,3,figsize=(15,5))
    ax[0].imshow(true_y[i])
    ax[1].imshow(pred_y[i])
    ax[2].scatter(true_y[i].ravel(), pred_y[i].ravel())
    ax[2].set_title(f"t={T[i]}")
    plt.show()
```


 ![output_78_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/59021082-365a-43f0-b72e-feaa6c14d4d7.png)
![output_78_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/f2f59e12-c43f-4dc0-b664-ce26a2e5a327.png)
![output_78_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/07f1cdda-057d-491a-aec3-443e4617ecd6.png)
![output_78_3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/4c21e695-ada1-42e1-8b9b-174aed6f1d98.png)
![output_78_4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/9e3ced9b-e445-4968-8589-3d45313545f2.png)
![output_78_5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/67366bbd-d399-4068-a5b5-d1d2da75e941.png)
![output_78_6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/05770bfa-08a5-42f0-9416-db15758a8501.png)
![output_78_7.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/abed0be7-be0e-4cb5-b5d1-1e379a6ea89a.png)
![output_78_8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/3637b9fc-0a87-46a9-8586-9295e9421d59.png)
![output_78_9.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/9c9b7787-5969-4b4c-adec-bfa181bf1d45.png)


今回は以上となります。
次回は深層ガウス過程などについて説明します。

https://qiita.com/tanaka_benkyo/items/7dbbaa508e1af87ac030

ご意見あればコメントお願いします