T-QARD Harbor

ボルツマンマシンの評価方法の検討

概要

「量子アニーリングによる画像生成の評価モデル」では，手書き文字画像の画像生成をD-Waveマシン，MCMC(マルコフ連鎖モンテカルロ法)を用いて行い，その画像を学習済みの分類器で評価していました．

ボルツマンマシンの評価は一般的にカルバック・ライブラー(KL)情報量や対数尤度を用いるのが一般的です．
確かに，対数尤度などでは生成画像そのものを評価することは難しいですが，ボルツマンマシンの評価はやはり尤度で評価するのが自然だと考えられます．

そこで，本記事では論文の評価指標と対数尤度による評価とを比較していきます．

文献情報

問題

ボルツマンマシンを用いて手書き数字画像を生成し，ボルツマンマシンを評価します．
ボルツマンマシンの評価として，上の論文では学習済みの分類モデルを用いて行っていましたが，本記事ではそれに加え，対数尤度による評価も行っていきます．

ボルツマンマシンの説明は「量子アニーリングによる画像生成の評価モデル」を参照してください．

実験1

まずは，論文で使用されているデータ，パラメータを用いて再現実験を行っていきます．

実験の流れとしては以下のとおりです．

  1. ボルツマンマシンを学習する
  2. 分類モデルを学習する
  3. 学習済みの両者を用いて，論文の評価指標の値を計算する．
  4. 対数尤度を計算する．

はじめに，必要なライブラリをインポートしておきましょう．

# 必要なライブラリのインポート
import os
import math
import itertools
import more_itertools
import random
import numpy as np
import tensorflow as tf
import matplotlib.pyplot as plt
from pathlib import Path
# プログレスバー用
from tqdm import tqdm
# PyTorch GPU使用のため(CPUでも実行可能)
import torch
import torch.nn.functional as F
from torch.utils.data import TensorDataset, DataLoader
# シード値の固定
def seed_everything(seed: int):
    tf.random.set_seed(seed)
    random.seed(seed)
    os.environ['PYTHONHASHSEED'] = str(seed)
    np.random.seed(seed)
    torch.manual_seed(seed)
    torch.cuda.manual_seed(seed)
    torch.backends.cudnn.deterministic = True
    torch.backends.cudnn.benchmark = False
    torch.use_deterministic_algorithms = True
    
SEED = 42
seed_everything(SEED)

def seed_worker(worker_id):
    worker_seed = torch.initial_seed() % 2**32
    np.random.seed(worker_seed)
    random.seed(worker_seed)

g = torch.Generator()
g.manual_seed(SEED)
<torch._C.Generator at 0x7f055cc33930>

データ準備

では，データのダウンロードから始めます．

# データダウンロード
! gdown https://drive.google.com/uc?id=18wXPH_lsKcDRohSI57iU4FA8Sf6aB5wk
Downloading...
From: https://drive.google.com/uc?id=18wXPH_lsKcDRohSI57iU4FA8Sf6aB5wk
To: /content/datasetdigit2.csv

  0% 0.00/1.19M [00:00<?, ?B/s]
100% 1.19M/1.19M [00:00<00:00, 99.0MB/s]

データは以下のような形をしています．

データの形 ： (\(N\), \(D\))

\(N\)：データ数
\(D\)：特徴量次元数 53次元 (手書き数字画像部分 48次元(8×6) + 数字の種類を表すOne-Hot表現 5次元)

手書き数字 5-9 を使用
例えば，8を表すOne-Hot表現は\((0, 0, 0, 1, 0)\)

では，データをロードして最初の5つの画像を表示してみます．

# データのロード
data = np.loadtxt('datasetdigit2.csv')
N, D = data.shape
# One-Hotの次元
N_SELECT = 5
# 手書き数字部分の次元
IMG_DIM = D - N_SELECT
# 手書き数字の画像サイズ
HEIGHT, WIDTH = 8, 6

# 最初の5つのデータを表示
plt.figure(figsize=(12, 5))
for i in range(5):
    num = i + N_SELECT
    plt.subplot(1, N_SELECT, i+1)
    plt.imshow(data[i, :IMG_DIM].reshape((HEIGHT, WIDTH)))

ボルツマンマシンを定義

ボルツマンマシンの概要については「量子アニーリングによる画像生成の評価モデル」を参照してください．

上のリンク先の記事には対数尤度の式を陽に示していなかったので，以下に示します．

対数尤度\(\log L(\theta)\)は以下のように変形できます．

$$
\begin{aligned}
\log L(\theta) &= \log \prod_{n = 1}^N p(\mathbf{x_n} | \theta) \\
&= \sum_{n = 1}^N \log p(\mathbf{x_n} | \theta) \\
&= \sum_{n = 1}^N \log \left( \frac{1}{Z(\theta)} \exp(-\Phi(\mathbf{x_n}, \theta)) \right) \\
&= \sum_{n = 1}^N \left( -\Phi(\mathbf{x_n}, \theta) – \log Z(\theta) \right) \\
&= \sum_{n = 1}^N -\Phi(\mathbf{x_n}, \theta) – N \log Z(\theta)
\end{aligned}
$$

ここで，

$$
\begin{aligned}
p(\mathbf{x} | \theta) &= \frac{1}{Z(\theta)} \exp (- \Phi(\mathbf{x}, \theta)) \\
\Phi(\mathbf{x}, \theta) &= – \sum_{i = 1}^{M} b_i x_i – \sum_{\{i, j\} \in E} w_{ij} x_i x_j \\
Z(\theta) &= \sum_{\mathbf{x}} \exp (-\Phi(\mathbf{x}, \theta))
\end{aligned}
$$

  • \(M\): ノード数
  • \(E\): エッジ集合

です．

以降，対数尤度という場合，便宜上両辺を\(N\)で割った\(\frac{1}{N}\log L(\theta)\)を指すこととします．

では，ボルツマンマシンのクラスを定義します．

主なメソッドは以下のとおりです．

  • train : 学習用メソッド
  • generate : 指定した数字を生成するメソッド
  • calc_likelihood : 対数尤度を計算するメソッド
class BolzmannMachine(torch.nn.Module):
    """
    n_feats: 画像部分の特徴量次元
    n_select: 使用する数字の種類
    momentum: モーメンタム法の係数
    weight_decay: l2 正則化の係数
    burn_in: バーンインタイム
    sampling_interval: サンプリング間隔
    eta: 学習率
    """
    # 論文の値がデフォルト値
    def __init__(self, 
                 img_dim, 
                 n_select,
                 output_dir,
                 eta=2.5e-2,
                 momentum=0.5, 
                 weight_decay=1e-5, 
                 burn_in=200,
                 sampling_interval=10,
                 batch_size=256):
        super().__init__()
        self.img_dim = img_dim
        self.n_select = n_select
        self.dim = img_dim + n_select
        self.momentum = momentum
        self.weight_decay = weight_decay
        self.burn_in = burn_in
        self.sampling_interval = sampling_interval
        self.eta = eta
        self.batch_size = batch_size

        self.output_dir = output_dir
        
        # CPUかGPUを使用するか判定
        # GPUが使えればGPU
        self.device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
        # 下三角 + 対角成分が0, それ以外は1の行列を定義
        self.mask = 1. - torch.tril(torch.ones(self.dim, self.dim, device=self.device))
        # 重み行列
        # 重みは下三角部分と対角部分の成分は持たないようにする
        w = torch.randn(self.dim, self.dim, device=self.device) * self.mask
        # バイアス
        b = torch.randn(self.dim, device=self.device)

        # この処理でクラスの属性としている．
        self.register_buffer('w', w)
        self.register_buffer('b', b)

        # 前ステップのvを保持するための変数
        self.v_w = 0
        self.v_b = 0
    
    # 訓練メソッド
    @torch.no_grad()
    def train(self, num_epochs, data):
        """
        num_epochs: エポック数
        x: 手書き文字部分の特徴量 (N, D) - N: データ数, D: データの次元
        y: 数字の種類を示すOne-Hot ベクトルの集合 (N, 5) - N: データ数
           今回の場合は5-9の数字を使用するため5次元 
        """
        x, y = data[:, :self.img_dim], data[:, self.img_dim:]
        x = torch.tensor(x, dtype=torch.float)
        y = torch.tensor(y, dtype=torch.float)
        # データ準備
        ds = TensorDataset(x, y)
        dl = DataLoader(
            ds,
            batch_size=self.batch_size,
            shuffle=True,
            worker_init_fn=seed_worker,
            generator=g,
        )
        # 学習
        for e in tqdm(range(num_epochs), total=num_epochs, desc='Epoch : '):
            for x, y in dl:
                x, y = x.to(self.device), y.to(self.device)
                self.update(x, y)
            # 10エポック毎にバイアスと重みを保存
            if (e + 1) % 10 == 0:
                torch.save(self.state_dict(), self.output_dir / f'ckpt_{e+1:04d}.pth')
    
    # 生成用メソッド
    @torch.no_grad()
    def generate(self, num):
        """
        num: 生成したい数字(5-9)
        """
        # One-Hot表現に変換
        y = F.one_hot(torch.LongTensor([num - self.n_select]).view(1, -1), num_classes=self.n_select)
        # ギブスサンプリングを行い, 返す
        return self.sampling(y)[-1, :-self.n_select].cpu().numpy().reshape(-1)

    @torch.no_grad()
    def generate_multi(self, num_list):
        y = F.one_hot(torch.LongTensor(num_list) - 5, num_classes=5)
        sample = self.sampling(y).cpu().numpy()[:, :-self.n_select]
        return sample, y.cpu().numpy()

    # パラメータ更新用メソッド
    def update(self, x, y):
        """
        x: 手書き文字部分の特徴量
        y: 文字の種類を表すOne-Hot特徴量
        """
        # ギブスサンプリングを行う
        samples = self.sampling(y)
        # 勾配を計算
        grad_b, grad_w = self.grad_b(x, y, samples), self.grad_w(x, y, samples)
        self.v_b = self.momentum * self.v_b + (1 - self.momentum) * grad_b
        self.v_w = self.momentum * self.v_w + (1 - self.momentum) * grad_w
        # 変化量を計算
        # 第1項目: 学習率 * モーメンタム法における速度
        # 第3項目: L2正則化項
        delta_b = self.eta * self.v_b - self.weight_decay * self.b
        delta_w = self.eta * self.v_w - self.weight_decay * self.w
        # パラメータ更新
        self.b += delta_b
        self.w += delta_w
        self.w *= self.mask
        return grad_b, grad_w

    # バイアス ｂ の勾配を求めるメソッド
    def grad_b(self, x, y, samples):
        """
        x: 手書き文字部分の特徴量
        y: 文字の種類を表すOne-Hot特徴量
        samples: サンプリングしたもの
        """
        N = x.size(0)
        # 手書き文字の部分とOne-Hotを結合: cat((N, C'), (N, 5)) => (N, C)
        _x = torch.cat([x, y], dim=-1)
        # データ方向に平均を取る: mean((N, C), dim=0) => (C,)
        e_data = _x.mean(dim=0)
        # サンプルについても同様に求める: mean((N, C), dim=0) => (C,)
        e_model = samples.mean(dim=0)
        # バイアスの勾配を返す
        return e_data - e_model
    
    # 重み w の勾配を求めるメソッド
    def grad_w(self, x, y, samples):
        """
        x: 手書き文字部分の特徴量
        y: 文字の種類を表すOne-Hot特徴量
        samples: サンプリングしたもの
        """
        # データ数
        N = x.size(0)
        # 手書き文字の部分とOne-Hotを結合: cat((N, C'), (N, 5)) => (N, C)
        _x = torch.cat([x, y], dim=-1)
        # 行列積を求めて，データ方向に平均を求める: mean((N, C, 1) x (N, 1, C), dim=0) => (C, C)
        e_data = (_x[:, :, None] @ _x[:, None, :]).mean(dim=0)
        # サンプルについても同様に求める: mean((N, C, 1) x (N, 1, C), dim=0) => (C, C)
        e_model = (samples[:, :, None] @ samples[:, None, :]).mean(dim=0)
        # 重みの勾配を返す
        return e_data - e_model
    
    # ギブスサンプリングを行うメソッド
    def sampling(self, y):
        """
        y: 文字の種類を表すOne-Hotベクトル
        """
        # サンプルを更新するインナー関数
        def _sample(x, b, w, y):
            # バイアス + 重み x 特徴量: (1, C, 1) + (1, C, C) x (N, C, 1) => (N, C, 1)
            _x = b[None, :, None] + w[None, :, :] @ x[:, :, None]
            # xの要素が1の確率 p(x = 1 | \theta)を求める
            p = torch.sigmoid(_x.squeeze(-1))
            rand = torch.rand(_x.size(0), self.dim, device=self.device)
            # U(0, 1) < p(x = 1 | \theta)の時1， それ以外を0
            x = (rand < p).float()
            x[:, -self.n_select:] = y
            return x
        # 最初のサンプルとしてランダムに初期化
        x = (torch.rand(y.size(0), self.dim, device=self.device) > 0.5).float()
        x[:, -self.n_select:] = y
        w = self.w + self.w.T
        # バーンインだけ更新
        for _ in range(self.burn_in):
            x = _sample(x, self.b, w, y)
        # サンプリング間隔だけ更新し，それをサンプルとする．
        for _ in range(self.sampling_interval):
            x = _sample(x, self.b, w, y)
        return x
    
    # 尤度計算用メソッド
    @torch.no_grad()
    def calc_likelihood(self, data, batch_size=2**14):
        """
        data: 計算に使用するデータ
        batch_size: (default: 2^14) 尤度計算のバッチサイズ
        """
        if isinstance(data, np.ndarray):
            data = torch.tensor(data, dtype=torch.float, device=self.device)
        w = self.w + self.w.T

        # 第1項目
        # (1, 1, C) x (N, C, 1) + (N, 1, C) x (1, C, C) x (N, C, 1) => (N, 1, 1) => `mean` => (1,)
        first = (self.b[None, None, :] @ data[:, :, None] + data[:, None, :] @ w[None, :, :] @ data[:, :, None]).mean()
        # 第2項目
        Z = 0
        for num in [5, 6, 7, 8, 9]:
            z_x = itertools.product([0, 1], repeat=self.img_dim)
            chunked_z_x = more_itertools.chunked(z_x, batch_size)
            for d_x in chunked_z_x:
                d = torch.zeros(len(d_x), self.dim, device=self.device)
                d[:, :-self.n_select] = torch.tensor(d_x, device=self.device, dtype=torch.float)
                d[:, num-10] = 1
                # (1, 1, C) x (N, C, 1) + (N, 1, C) x (1, C, C) x (N, C, 1) => (N, 1, 1)
                minus_phi = (self.b[None, None, :] @ d[:, :, None] + d[:, None, :] @ w @ d[:, :, None])
                Z += torch.exp(minus_phi).sum()
        second = torch.log(Z)
        # 尤度 = 第1項目 - 第2項目
        likelihood = first - second
        return float(likelihood.cpu().numpy())

ボルツマンマシンの学習

実際にボルツマンマシンを学習していきます．

# 各種パラメータを設定
N_EPOCHS = 4000
exp1_dir = Path('./exp1')
exp1_dir.mkdir(parents=True, exist_ok=True)
# ボルツマンマシン初期化
bm = BolzmannMachine(img_dim=IMG_DIM, n_select=N_SELECT, output_dir=exp1_dir)
# 学習開始
bm.train(N_EPOCHS, data)
Epoch : 100%|██████████| 4000/4000 [07:41<00:00,  8.66it/s]

次に，100エポックごとの生成結果を見ていきます．

# 100エポック毎の生成結果をプロット
n_gen = N_EPOCHS // 100
plt.figure(figsize=(8, 2 * n_gen))
for n in range(n_gen):
    bm.load_state_dict(torch.load(exp1_dir / f'ckpt_{(n + 1) * 100:04d}.pth'))
    for i in range(5):
        num = i + N_SELECT
        gen = bm.generate(num)
        plt.subplot(n_gen, N_SELECT, N_SELECT * n + i + 1)
        plt.imshow(gen.reshape((HEIGHT, WIDTH)), aspect='auto')
  • 序盤は学習が進むに連れて生成画像も数字に近づいていそうですが，途中で崩壊しているため，安定性に乏しいといえます．

分類モデルの学習

次に，分類モデルで評価を行うためはじめに分類モデルを学習させます．

# 分類モデル用にライブラリをインポート
from sklearn.model_selection  import train_test_split

from keras.callbacks import ModelCheckpoint
from keras.models import Sequential, load_model
from keras.layers import Dense, Dropout, InputLayer
# 再度データを用意
X = data[:, :IMG_DIM]
y = data[:, IMG_DIM:]
X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.1, shuffle=True, random_state=SEED)

X_train = X_train.astype(np.float32)
X_valid = X_valid.astype(np.float32)

分類モデルのニューラルネットワークを定義します．

def build_model(img_dim, n_select):
    model = Sequential()
    model.add(InputLayer(input_shape=(img_dim,)))
    model.add(Dense(512, activation='relu'))
    model.add(Dropout(0.4))
    model.add(Dense(512, activation='relu'))
    model.add(Dropout(0.4))
    model.add(Dense(n_select, activation='softmax'))
    return model

実際に分類モデルを学習します．

# 学習エポック数
N_EPOCHS_CLS = 20
# バッチサイズ
BATCH_SIZE = 32

# モデル初期化
model = build_model(IMG_DIM, N_SELECT)
model.compile(loss='categorical_crossentropy', optimizer='adam', metrics=['accuracy'])
ckpt_callback = ModelCheckpoint(
    filepath=exp1_dir / 'best.h5',
    monitor='val_accuracy',
    mode='max',
    save_best_only=True
)
# 学習
history = model.fit(
    X_train,
    y_train,
    batch_size=BATCH_SIZE, 
    epochs=N_EPOCHS_CLS, 
    verbose=0, 
    validation_data=(X_valid, y_valid),
    callbacks=[ckpt_callback]
)

LossとAccuracyをエポックごとの推移で見ていきます．

plt.plot(range(1, N_EPOCHS_CLS+1), history.history['loss'],  marker='.', label='train')
plt.plot(range(1, N_EPOCHS_CLS+1), history.history['val_loss'], marker='.', label='valid')
plt.legend(loc='best', fontsize=10)
plt.grid()
plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.show()
plt.plot(range(1, N_EPOCHS_CLS+1), history.history['accuracy'], marker='.', label='train')
plt.plot(range(1, N_EPOCHS_CLS+1), history.history['val_accuracy'], marker='.', label='valid')
plt.legend(loc='best', fontsize=10)
plt.grid()
plt.xlabel('Epoch')
plt.ylabel('Accuracy')
plt.show()
  • LossとAccuracyの推移を見ても訓練(train)，検証(valid)間で差は大きくなく問題なく学習できているといえます．

分類モデルによる評価

では，実際に評価していきます．
論文の評価指標\(R_a\)は以下のように表されるのでした．

$$
R_a = \frac{D_m}{D_t}
$$

  • \(D_m\): 生成で指定した数字と分類結果が一致した数
  • \(D_t\): 生成した画像の総数

論文の評価指標で10エポックずつ，モデルを評価してみましょう．

  • 5 – 9までの数字を100枚ずつ生成し，分類モデルでAccuracyを計算します．
n_gen = N_EPOCHS // 10
epoch_list = list()
score_list = list()
bm = BolzmannMachine(img_dim=IMG_DIM, n_select=N_SELECT, output_dir=exp1_dir)
model = load_model(exp1_dir / 'best.h5')
for n in range(n_gen):
    bm.load_state_dict(torch.load(exp1_dir / f'ckpt_{(n + 1) * 10:04d}.pth'))
    num_list = torch.LongTensor([5] * 100 + [6] * 100 + [7] * 100 + [8] * 100 + [9] * 100)
    gen, y = bm.generate_multi(num_list)
    epoch_list.append((n + 1) * 10)
    score = model.evaluate(gen, y, verbose=0)
    score_list.append(score)
_, scores = tuple(zip(*score_list))
plt.plot(epoch_list, scores, marker='.', label='train')
plt.grid()
plt.xlabel('Epoch')
plt.ylabel('R_a')
plt.show()
  • 生成画像でみていたように途中で学習が崩壊し，振動している様子が見て取れます．
  • 学習を続ければ評価指標の値が上昇し続けていくわけではなさそうです．
  • 最高値は論文より高くなっていますが，論文の結果は安定して評価指標が推移していたため，そういった点では再現を行うことができませんでした．

エポックごとの評価指標の値がわかったため，一番評価指標の値が良いモデルの生成画像を見てみましょう．

best_epoch = epoch_list[np.argmax(scores)]
print(f'Best Epoch : {best_epoch}')
Best Epoch : 1690
plt.figure(figsize=(8, 6))
bm = BolzmannMachine(img_dim=IMG_DIM, n_select=N_SELECT, output_dir=exp1_dir)
bm.load_state_dict(torch.load(exp1_dir / f'ckpt_{best_epoch:04d}.pth'))
for i in range(3): 
    for j in range(N_SELECT):
        num = j + N_SELECT
        gen = bm.generate(num)
        plt.subplot(3, N_SELECT, 5 * i + j + 1)
        plt.imshow(gen.reshape((HEIGHT, WIDTH)), aspect='auto')

対数尤度計算

次に，対数尤度を求めていくのですが，ここで問題があります．

対数尤度を求めるためにはボルツマン分布の分配関数(\(Z(\theta)\))を正確に評価する必要があります．
ここで，画像の次元数は48次元なので，一回の尤度計算に\(2^{48}\)個の場合分けを計算する必要が出てきます．
これは，たとえGPUで計算したとしても現実時間で終わらせることが難しいです．

そこで，本記事では次元数を小さくしたデータを人工的に作成し，そのデータを用いることで尤度計算を行います．
そのため，新しいデータに対してもう一度学習及び評価を行っていきます．

実験2

上述したように，論文データの対数尤度を計算するには計算時間の問題がありました．
そこで，実験2として論文データよりも小さいデータに対して同じ実験を行い，評価指標\(R_a\)と対数尤度の計算を行っていきます．

データ準備

5-9の数字について，20次元(5×4)のデータを作成していきます．

d5 = np.array([[0, 1, 1, 0],
               [1, 0, 0, 0],
               [0, 1, 1, 0],
               [0, 0, 0, 1],
               [0, 1, 1, 0]])
edge_flag5 = [1, 1, 1, 1, 1, 1]
plt.imshow(d5)
d6 = np.array([[0, 1, 1, 0],
               [1, 0, 0, 0],
               [0, 1, 1, 0],
               [1, 0, 0, 1],
               [0, 1, 1, 0]])
edge_flag6 = [1, 1, 1, 1, 1, 1]
plt.imshow(d6)

d7 = np.array([[0, 1, 1, 0],
               [1, 0, 0, 1],
               [0, 0, 0, 1],
               [0, 0, 0, 1],
               [0, 0, 0, 1]])
edge_flag7 = [1, 1, 0, 1, 0, 1]
plt.imshow(d7)
 
d8 = np.array([[0, 1, 1, 0],
               [1, 0, 0, 1],
               [0, 1, 1, 0],
               [1, 0, 0, 1],
               [0, 1, 1, 0]])
edge_flag8 = [1, 1, 1, 1, 1, 1]
plt.imshow(d8)
d9 = np.array([[0, 1, 1, 0],
               [1, 0, 0, 1],
               [0, 1, 1, 0],
               [0, 0, 0, 1],
               [0, 1, 1, 0]])
edge_flag9 = [1, 1, 1, 1, 1, 1]
plt.imshow(d9)

5枚では足りないので，不自然にならないように画素を増やして枚数を増やします．

_data = list()
num_list = [d5, d6, d7, d8, d9]
edge_flag_list = [edge_flag5, edge_flag6, edge_flag7, edge_flag8, edge_flag9]
edge_idx_list = [(0, 0), (0, 3), (2, 0), (2, 3), (4, 0), (4, 3)]
labels = np.eye(5)
for i in range(5):
  num = num_list[i]
  edge_flag = edge_flag_list[i]
  for j, flag in enumerate(edge_flag):
    d = num.copy()
    if flag:
      d[edge_idx_list[j]] = 1
    d = np.concatenate([d.reshape(-1), labels[i]], axis=-1)
    _data.append(d)
  d = np.concatenate([num.reshape(-1), labels[i]], axis=-1)
  _data.append(d)
data = np.stack(_data, axis=0)
data.shape
(35, 25)
# 各種定数を定義しておく
IMG_DIM = 20
N_SELECT = 5
HEIGHT = 5
WIDTH = 4

実際に使用するデータについてプロットしてみます．

plt.figure(figsize=(8, 12))
for i in range(data.shape[0]):
  plt.subplot(7, 5, i+1)
  plt.imshow(data[i, :IMG_DIM].reshape((HEIGHT, WIDTH)))

ボルツマンマシンの学習

すでに実験1でボルツマンマシンは定義したので，それを上で作成したデータで学習します．

N_EPOCHS = 4000
exp2_dir = Path('./exp2')
exp2_dir.mkdir(parents=True, exist_ok=True)
bm = BolzmannMachine(img_dim=IMG_DIM, n_select=N_SELECT, output_dir=exp2_dir)
bm.train(N_EPOCHS, data)
Epoch : 100%|██████████| 4000/4000 [02:01<00:00, 32.91it/s]

同様に100エポックごとの生成画像を見てみましょう．

# 100エポック毎の生成結果をプロット
n_gen = N_EPOCHS // 100
plt.figure(figsize=(7, 2 * n_gen))
for n in range(n_gen):
    bm.load_state_dict(torch.load(exp2_dir / f'ckpt_{(n+1) * 100:04d}.pth'))
    for i in range(5):
        num = i + N_SELECT
        gen = bm.generate(num)
        plt.subplot(n_gen, N_SELECT, N_SELECT * n + i + 1)
        plt.imshow(gen.reshape((HEIGHT, WIDTH)), aspect='auto')
  • 5, 6, 8, 9は非常に似ているため，ところどころ混ざってしまっているものがあります．
  • 7については比較的問題なく学習できていそうです．

分類モデルの学習

分類モデルについても，上のデータを用いて同様の学習をさせます．

# 再度データを用意
X = data[:, :IMG_DIM]
y = data[:, IMG_DIM:]
X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.1, shuffle=True, random_state=SEED)

X_train = X_train.astype(np.float32)
X_valid = X_valid.astype(np.float32)
N_EPOCHS_CLS = 20
BATCH_SIZE = 32

model = build_model(IMG_DIM, N_SELECT)
model.compile(loss='categorical_crossentropy', optimizer='adam', metrics=['accuracy'])
ckpt_callback = ModelCheckpoint(
    filepath=exp2_dir / 'best.h5',
    monitor='val_accuracy',
    mode='max',
    save_best_only=True
)
history = model.fit(
    X_train,
    y_train,
    batch_size=BATCH_SIZE, 
    epochs=N_EPOCHS_CLS, 
    verbose=0, 
    validation_data=(X_valid, y_valid),
    callbacks=[ckpt_callback]
)

LossとAccuracyのエポックごとの推移を見てみます．

plt.plot(range(1, N_EPOCHS_CLS+1), history.history['loss'],  marker='.', label='train')
plt.plot(range(1, N_EPOCHS_CLS+1), history.history['val_loss'], marker='.', label='valid')
plt.legend(loc='best', fontsize=10)
plt.grid()
plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.show()
plt.plot(range(1, N_EPOCHS_CLS+1), history.history['accuracy'], marker='.', label='train')
plt.plot(range(1, N_EPOCHS_CLS+1), history.history['val_accuracy'], marker='.', label='valid')
plt.legend(loc='best', fontsize=10)
plt.grid()
plt.xlabel('Epoch')
plt.ylabel('Accuracy')
plt.show()
  • そもそもの学習データ数が少ないため，スコアの推移が極端になっています．
  • ただ終盤のエポックでは問題なく分類できていそうです

分類モデルによる評価

では10エポックごとに論文の評価指標の推移を見ていきます．

n_gen = N_EPOCHS // 10
epoch_list = list()
score_list = list()
bm = BolzmannMachine(img_dim=IMG_DIM, n_select=N_SELECT, output_dir=exp2_dir)
model = load_model(exp2_dir / 'best.h5')
for n in range(n_gen):
    bm.load_state_dict(torch.load(exp2_dir / f'ckpt_{(n+1) * 10:04d}.pth'))
    num_list = torch.LongTensor([5] * 100 + [6] * 100 + [7] * 100 + [8] * 100 + [9] * 100)
    gen, y = bm.generate_multi(num_list)
    epoch_list.append((n+1) * 10)
    score = model.evaluate(gen, y, verbose=0)
    score_list.append(score)
_, scores = tuple(zip(*score_list))
plt.figure(figsize=(8, 4))
plt.plot(epoch_list, scores, marker='.')
plt.grid()
plt.xlabel('Epoch')
plt.ylabel('R_a')
plt.show()
  • 実験1の論文データのときと比較すると値の推移は非常に安定しています．

一番評価指標の値が良いモデルの生成画像を見てみましょう．

best_epoch = epoch_list[np.argmax(scores)]
print(f'Best Epoch : {best_epoch}')
Best Epoch : 3920
plt.figure(figsize=(6, 7))
bm = BolzmannMachine(img_dim=IMG_DIM, n_select=N_SELECT, output_dir=exp2_dir)
bm.load_state_dict(torch.load(exp2_dir / f'ckpt_{best_epoch:04d}.pth'))
for i in range(3):
    for j in range(N_SELECT):
        num = j + N_SELECT
        gen = bm.generate(num)
        plt.subplot(3, N_SELECT, 5 * i + j + 1)
        plt.imshow(gen.reshape((HEIGHT, WIDTH)), aspect='auto')
  • 問題なく生成できていそうです．
  • 端のピクセルが多くなるのはデータ拡張の影響だと考えられます．

尤度計算

では，実験1ではできなかった尤度計算に移ります．

尤度計算の実装はすでに定義したボルツマンマシンのクラス内で行ったので，20エポックずつ尤度計算をして推移を確認してみましょう．

N = N_EPOCHS // 20
epoch_list = list()
likelihood_list = list()
bm = BolzmannMachine(img_dim=IMG_DIM, n_select=N_SELECT, output_dir=exp2_dir)
for n in tqdm(range(N)):
    epoch_list.append((n + 1) * 20)
    bm.load_state_dict(torch.load(exp2_dir / f'ckpt_{(n+1) * 20:04d}.pth'))
    likelihood = bm.calc_likelihood(data, batch_size=2**18)
    likelihood_list.append(likelihood)
100%|██████████| 200/200 [18:39<00:00,  5.60s/it]
plt.plot(epoch_list, likelihood_list, marker='.')
plt.grid()
plt.xlabel('Epoch')
plt.ylabel('Likelihood')
plt.tight_layout()

学習が進むにつれて尤度は上昇していることがわかりますが，評価指標\(R_a\)と比べると上昇は速いことがわかります．これは，評価指標\(R_a\)の計算にはボルツマンマシンの生成結果を使用しているため，生成のランダム性により，少ないエポック数のモデルでは画像が崩れやすく，\(R_a\)の値は低くなりやすい一方で，尤度計算には実際のデータを使用して評価しているため，そこにランダム性はなくそこそこの品質の画像を生成できるだけの重みとバイアスがあれば尤度は上がりやすいため，このような違いが出たと考えられます．

そういった点から，分類モデルの性能に依存しているものの，生成モデルとしての性能を確認するという点では，生成画像そのものを使用している論文の評価指標の方が優れていると言えます．

また，尤度計算では実験1でわかったようにある程度の次元のデータになると，尤度自体の計算が困難な場合があり，計算量の観点からも評価指標\(R_a\)の方が優れています．

結論

本記事では論文の評価指標と対数尤度を比較しましたが，論文の評価指標\(R_a\)の方が生成画像そのものを評価に組み込め，計算時間も少ないため，生成モデルとしての性能を測る上では優れていると言えます．
一方で学習は対数尤度ベースで行われているため，学習が適切に行われているかを確認する上では尤度を見ることは妥当であるといえます．
しかしながら，尤度は非常に計算量が大きいため，そういった確認にも簡単に計算を行うことの出来る論文の評価指標\(R_a\)を使用することは可能だと考えられます．

あとがき

論文のデータでは学習が安定しないという問題が起きましたが，数字画像の生成自体は行うことができました．
学習の安定性の問題は，論文のパラメータの値が著者実装に依存していることから起きていると思うので，本実装用にパラメータを調整することで改善すると考えられます．

また，ボルツマンマシンの評価で使われる対数尤度については自作のデータセットを用いることで論文の評価指標と比較することができました．
対数尤度の計算量についてはどうしても大きくなってしまうため，そういった点では次元が大きくなってもそこまで計算時間の大きくならない論文の評価指標の方が様々な面で優れていると感じました．

今回は5-9の数字を使用しましたが，今後は0-9のようなクラス数の多いものや，データの次元数が大きいデータセット，また，論文で行われていたようにD-Waveを使用した実験に挑戦したいと思います．

本記事の担当者

清水怜央

Post Views: 237

© 2018-2023 T-QARD: Center for Quantum Application Research and Development, Graduate School of Information Sciences, Tohoku University