グラフ機械学習と強化学習について

主にグラフ機械学習や強化学習手法を記載します。

モンテカルロ木探索についてのまとめ (データ構造とアルゴリズム Advent Calendar 2020)

 この記事は「データ構造とアルゴリズム Advent Calendar 202020日目の記事です。

19日目は@takilogさん; グラフ上の合流に関する問題とアルゴリズム ,
21日目は@tmaeharaさんです。

概要

モンテカルロ木探索 (MCTS) は、木探索にモンテカルロ (ランダム) 要素を加味した評価関数を用いることで、効率よく探索を行うことのできるアルゴリズムです1。特に囲碁など2プレイヤーの完全情報ゼロサムゲームで使用されており、プロ囲碁棋士を破ったAlphaGO, AlphaZeroや、Atariゲーム、将棋、チェスなどでも圧倒的スコアを出したMuZero2もこの方法を用いてます。別のタスクでは、バンディット問題やモデルベース強化学習においての探索にも使われたりします。ポーカーのような不完全情報ゲームでもMCTSを拡張した方法もあります3

データ構造としては木構造になるので、初心者の方は実装するのが少し難しいかもしれませんが、幅優先探索深さ優先探索を理解しているならば組み込むのは難しくありません。また、幅優先探索に評価関数を組み込んだ最良優先探索の考え方も用います。有名なアルゴリズムA*アルゴリズムなどです。なお、一番難しい箇所は最適な評価関数を設計するところです。効率よく並列計算をしようとすると難しいです4

様々な人がMCTSを解説していると思いますが、私も記事を書いてみようと思います。ボードゲームに限らず、探索が必要な場面では使える技術ですので、覚えておいて損はないと思います(計算コストは非常に大きいですが・・・)。

木探索と強化学習

木探索は、下図のように、望む最終結果を得るために最適なノードを選択しながらルートを探索していきます。通常の幅優先探索深さ優先探索では、組み合わせ数が多すぎて探索しきれません。

beam_search

引用元 DL2AI: ビームサーチの例

評価関数値に基づきノードを選択し、次のノードをさらに選択するという方法の幅優先探索です。貪欲法では、評価関数の中で最も良い値を選び続けていきます。この図は、最終的な文字列を探索していることを表しており、ビームサーチと呼ばれる方法で、ビーム幅2(評価値TOP2のノードを選択する)で探索をしています。

例えば、ユーザーが望む文章を生成したいとき、木探索とRNNによる評価値を組み合わせることで文章生成ができます。

評価関数は通常、最終結果をもとに更新されることが多いです。基本的にはこの枠組みで問題を解いていくことになります。

バンディット問題

MCTSを適用する例題としてバンディット問題を考えてみます。これは以前、強化学習入門:多腕バンディット問題という記事を書いています。要点は、K個のスロットマシンにあたりのある台が1つあり、多数の試行を通してどのスロットの腕を選んだらよいかを探索する問題です。こちらも木探索を使ってスロットのアームを選択しており、最終結果に基づき評価関数値を更新するという方法になっています。

モデルベース強化学習

ここで、強化学習機械学習手法の1つであり、「環境」からの「状態」「報酬」を受け取り、それをもとにエージェント(モデル)が「行動」を行う方策を学習していきます。モデルベースでは状態遷移が既知、すなわち上図での枝が分かっている状態、マルコフ決定過程(MDP)で表現されるものです。しっかりした理論に基づいて体系化されていますが、ここでは省略します。評価関数・報酬の設計に関わるところで強化学習の考え方を使っていくことになります。

アルゴリズム

MCTSは4つのステップからなります。Wikipediaから図を引用します。

mcts

Monte Carlo tree search - Wikipedia

  • Selection: ルートから始めて、子ノードを選択していきます。選択基準には複数あります。
  • Expansion: 選択したノードから、次のノードをどれにするか選択します。シミュレーション後の結果をもとに展開します。
  • Simulation: 一様分布などランダムにノードを選択していきリーフノードまで計算を行います。このSimulation何度も繰り返します。rolloutと呼ばれます。
  • Backup: シミュレーションで得られた最終結果をもとに評価値を更新していきます。この図では(勝った回数) / (試行回数)となっており、最終的に勝つことができたノードほど選択されやすくなります。

この操作を繰り返していくことで、木を拡大していきます。十分なシミュレーション回数を行うと、よりよい探索を行うことができるようになります。

評価関数

探索と活用(Exploration and Exploitation)のバランスがうまく取れた評価関数を使っていきます。バンディット問題やベイズ最適化といった場面でも出てくる話です。「次の試行(実験)点をどこから選べばよいのか」というプラニングの問題では、この探索と活用というジレンマに悩まされることになります。機械学習モデルは、基本的に学習データ範囲内(内装領域)でしかモデル精度を保証できません。学習データ外に素晴らしい候補点があった場合は、探索を行ってその点を見つけていくことになります。そのための方法が複数提案されています。

UCB

代表的なものはUpper Confidence Bound (UCB)です。

$$ \begin{align} a_i = \arg \max_i (\frac{w_i}{n_i} + c \sqrt{\frac{logN}{n_i}}) \end{align} $$

 w_iが勝った回数、 n_iがノードを選択した回数、$N_i$が全シミュレーション回数、$c$探索のパラメータで、この値が大きいほどランダムなノードが選ばれやすくなります。通常は$\sqrt2$です。第1項は活用を表します。あるノードを選んだときにどれだけ勝っているかの確率で、この手だけを選んでいけば、ある種の貪欲法になります。第2項は探索を表します。選ばれた回数が多いほど、この値は小さくなるので、少ない手が選ばれやすくなります。

Crazy Stone

Crazy Stoneという囲碁のプログラムがありますが、この選択基準は次の式になります。

$$ \mu_i = \exp \left(-2.4 \frac{\mu_0 - \mu_i}{\sqrt{2 ( σ_{0}^{2} + σ_{i}^{2} ) }} \right) + \varepsilon_i \ $$

この式はガウス分布を仮定して得られるものの近似です(2.4が正規分布と一致するように選ばれている)。これは、n腕バンディットで用いられるボルツマン分布と類似したものです。εは0にならないようにする定数で、下式で定義されています。

$$ \begin{align} \varepsilon_i = \frac{0.1 + 2^{-1}+a_i}{N} \end{align} $$

ここでa_iはあたりなら1、その他は0という囲碁のルールに基づいた値で、ドメイン知識で決定されています。また、バリアンスは下式で計算します。

\begin{align} σ^{2} = \frac{Σ_{2} - S \mu^{2} + 4P^{2}}{S+1} \end{align}

ここで、Pが盤面の点数。Σ_{2} はノードのsum of squaredです。Sはシミュレーション回数です。不確定要素を考慮するためにこの式が得られています。かなり工夫がされています。

AlphaZero

arxiv.org

AlpahZeroは評価関数をニューラルネットワークで決定しています。ゲームのルールのみ教えています。個人だとまず学習できないレベルのマシンスペックによって計算されています。

MCTS疑似コード

以下のコードを見るとAlphaZeroで何をやっているのか概要が分かるかと思います。基本的にはMCTSを実行していますが、どのノード(手)を選択するかでニューラルネットワークで決定されたスコアをもとに決定されているのが分かるかと思います。

def run_mcts(config: AlphaZeroConfig, game: Game, network: Network):
  root = Node(0)
  evaluate(root, game, network) // ニューラルネットワークで値を取得して方策を予測
  add_exploration_noise(config, root)

  for _ in range(config.num_simulations):
    node = root
    scratch_game = game.clone()
    search_path = [node]

    while node.expanded():
      action, node = select_child(config, node)  // UCBスコアに基づいてノードの選択
      scratch_game.apply(action)
      search_path.append(node)

    value = evaluate(node, scratch_game, network)
    backpropagate(search_path, value, scratch_game.to_play())
  return select_action(config, game, root), root

Selection: UCB score

\begin{align} score &= Q(s, a) + C(s)P(a|s) \frac{\sqrt{N(s)}}{1+N(s, a)} \\ C(s) &= \log \frac{1+N(s) + c_{base}}{c_{base}} + c_{init} \end{align}

Q(s,a): 状態 s で行為 a を行った際の平均報酬, P(a|s)は方策です。すなわち、状態 s で行為 a を選択する確率で、ニューラルネットワークの出力となっています。 Nは訪問回数です。実際の疑似コードでは下のようになります。

# The score for a node is based on its value, plus an exploration bonus based on
# the prior.
def ucb_score(config: AlphaZeroConfig, parent: Node, child: Node):
  pb_c = math.log((parent.visit_count + config.pb_c_base + 1) /
                  config.pb_c_base) + config.pb_c_init
  pb_c *= math.sqrt(parent.visit_count) / (child.visit_count + 1)

  prior_score = pb_c * child.prior
  value_score = child.value()
  return prior_score + value_score

Backupにて更新されていく値になっています。詳細はブログ記事や論文を参照ください。

↓ブログ記事

deepmind.com

MuZero

MuZeroもAlphaZeroと同様に評価関数をニューラルネットワークで決定しています。AlpahZeroとの違いはAtraiやゲームルールすら教えずゼロから学習している点が違います。

  • Selection

探索木のノードごとに内部状態sがあり、sからの行動aはエッジ(s, a)で保存されています。 $\{N(s, a), Q(s, a), P(s, a), R(s, a), S(s, a)\}$は、それぞれ訪問回数N, 報酬の平均Q, 方策P, 報酬Rと状態遷移 Sを表しています。完全に強化学習の枠組みで得られる値です。

f:id:udnp:20201219211529p:plain

c1, c2は、よく訪問されるノードのQ(s, a)値に関連するP(s, a)の影響をコントロールするために用いられる定数です。

  • Expansion

最終ステップ$l$で、葉ノードに到達したとき、報酬と状態が別の関数gにて計算され、同時にテーブルに保存されます。また、方策とQ値は予測関数fで計算され、新しいノードが探索木に追加されることになります。そして, $N(s^{l}, a)=0, Q(s^{l}, a)=0, P(s^{l}, a)=p^{l}$で初期化されます。

  • Backup

累積報酬に関しては、時間割引率$\gamma$を考慮した推定値です。 f:id:udnp:20201220102234p:plain

$(s^ {k -1}, a^ k)$に対するQ, Nに関しては、次のように更新しています。

f:id:udnp:20201220101513p:plain

その他の応用

MCTSを使うとき、ゲームでは勝ったか負けたかのゼロサムですが、最終スコアが連続値でも問題なくつかうことができます。例えば、最終スコアが10になるように手を選択したり、特定の文字となるように探索したりなどです。 前者の場合では最終結果がどれだけ10に近いかで評価関数を更新してきます。後者でも、文字列を何らかの関数でスカラーにすれば同じプロセスとなります。Chemoinfomaticsの分野では、分子構造を表す文字列生成 → ChemTSといった方法もあります。

最後に

MCTSの考え方は非常にシンプルなのですが、selection, expansionでどのようにノードを選択し拡大していくかという評価関数を作りこんでいくところが非常に難しいです。backupのところでも、単純に訪問回数に+1をするだけでなく行動価値の値を更新する際など、かなりの工夫が必要になっていきます。AlphaZeroやMuZeroに関する処理を実装することでさらに理解を深めていきたいです。

References


  1. Remi Coulom. Efficient selectivity and backup operators in monte-carlo tree search. In International conference on computers and games, pages 72–83. Springer, 2006.

  2. Mastering Atari, Go, Chess and Shogi by Planning with a Learned Model, https://arxiv.org/abs/1712.01815

  3. Combining Prediction of Human Decisions with ISMCTS in Imperfect Information Games, AAMAS ‘18: Proceedings of the 17th International Conference on Autonomous Agents and MultiAgent SystemsJuly pages 1874–1876, 2018.

  4. On Effective Parallelization of Monte Carlo Tree Search, arXiv:2006.08785

ECFPとNeural-fingerprintの比較

はじめに

ケモインフォマティクスやマテリアルズインフォマティクスでは良く使用されるExtended Connectivity Circular Fingerprints (ECFP)及び Neural Graph Fingerprint (NFP) を実装してみました。ZincデータセットのlogP-SAを予測することで、両手法を比較しています。ECFPとNFPは下論文で提案された方法です。

ECFP自体は類似する方法がRDKitやCDKに実装されおり、正確かつ高速に使うことができます。NFPは実装上の都合上、CPUオンリーでバッチに対応していません。私の実装したものは、あくまで理解のためですので、アルゴリズムは多少間違っているかもしれません。ECFPについてはQiitaにも載せています。

フィンガープリント

化学の分野で用いられるフィンガープリントは、分子構造を表現する方法の1つで2次元構造に基づいています。主に、分子の「類似性検索」や定量的構造活性・物性相関(QSAR, Quantitative Structure-Activity Relationship or QSPR, Quantitative Structure-Property Relationship)モデル構築に使用されます。分子構造はSMILESと呼ばれる文字列情報で表現されます。また、分子はグラフ構造でも表せ、原子情報をノード、結合情報をエッジとして表現することもできます。原子=ノード、結合=エッジとして考えて構いません。私がよく用いるフィンガープリントはECFPですが、グラフニューラルネットワークによるフィンガープリントは多数ありますが、本記事では、ECFPとNFPについて記載します。

f:id:udnp:20200408183011p:plain

図1はECFPとNFPの概念図です。左図はECFPの概略図で分子グラフのメッセージパッシングによりノード(原子の)情報が伝達され、最終的に固定長の0, 1のビットに変換されます。右図はNFPで情報伝達はニューラルネットワークで重み付けされ情報が伝達されます。

アルゴリズム

実際にアルゴリズムを見るのがわかりやすいと思います。

f:id:udnp:20200408182955p:plain

Circular fingerprints (ECFP)

上図のアルゴリズムを詳しく見ていきます。まず、分子の半径Rとfingerprintの長さSを決めます。この半径は隣接ノードへ伝達する回数で、この操作はメッセージパッシングと呼ばれます。半径2なら2つ隣りのノードまで、R=3なら3分子離れたノード情報を伝達していく操作になります。fingerprintの配列fを0で初期化しておきます。固定長Sは1024や2048 bitsとなることが多いです。分子中の原子ごとにAtom Propertiesであるノード特徴量をリストに追加していきます。例えば、

  • 原子番号
  • 隣接重原子(水素を除いた原子)の数
  • 原子価
  • 環構造であるかどうか

などです。他にも追加していくことができます。隣接ノード情報を注目しているノードへ集め、concatenationしていきます。これをハッシュに通し、fingerprintの長さSで余りを求めます。この余りのインデックス部分に1を書き込みます。これを全原子に対して行い(for each atom in molecule)、何度も繰り返す(for 1 to R)ことで情報をメッセージパッシングさせることができます。この操作はMorganアルゴリズムとも呼ばれています。グラフ同型問題で非常に高い識別精度をもつWeisfeiler-Lehman graph isomorphism testと操作は非常に似ています。ECFPやNFPでは、この特徴量をどう決めるかによって結果が大きく変わるので、正しく決めることが大事です。

Neural graph fingerprints (NFP)

NFPとECFPとの違いは、concatenate以降のところがニューラルネットワークに決定される所になります。ECFPでは単純に結合し、ハッシュを通して、その余りのインデックスを1としている決定的な方法でした。NFPでは学習によってノード情報をランク付けし、ソフトマックスで確率的にインデックスを決定し、フィンガープリントの値としています。そのため、0, 1のベクトルではなく連続値となります。論文名のConvolutional Networks on Graphともあるように、 これはグラフ畳み込みニューラルネットワーク(GCN)の考え方です。現在はより一般的なフレームワーク(DGL, DeepChem, GraphNetsやPyTorch geometricなど)があるので、そちらで実装した方が良いと思います。NFPはGCNとほとんど同じアルゴリズムですが、NFPでは最終層でsoftmaxを取り、sum poolingを行っている点が異なっています。

ECFPの実装

分子を扱う場合はやはり、RDKitがとても便利です。RDKitでノード特徴量を作っていきます。

RDKit版 ECFP

from __future__ import print_function
from rdkit import Chem
from rdkit.Chem import AllChem, Draw
import numpy as np

smiles = 'c1ccccn1'

mol = Chem.MolFromSmiles(smiles)
bit_morgan1 = {}
fp1 = AllChem.GetMorganFingerprintAsBitVect(mol, radius=2, nBits=12, bitInfo=bit_morgan1)
bit1 = list(fp1)
print(bit1)
>> [0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0]

RDKitを用いれば、このように簡単にSMILESからECFPを行うことができます。このフィンガープリントを説明変数として目的変数を回帰したものが、QSAR / QSPRと呼ばれるものになります。

参考までに分子を描画します。

def mols2grid_image(mols, molsPerRow=1):
    mols = [e if e is not None else Chem.RWMol() for e in mols]
    for mol in mols:
        AllChem.Compute2DCoords(mol)

    return Draw.MolsToGridImage(mols, molsPerRow=molsPerRow, subImgSize=(150,150))

mols2grid_image([mol], molsPerRow=2)

f:id:udnp:20200408183032p:plain

自作版 ECFP

自作するためには、SMILESから原子の情報とエッジの情報を表す特徴を抽出しなければなりません。複数分子に対応できるようにmolオブジェクトをリストで渡すようにしていますが、その後の処理は行っていません。

  • 隣接行列
def GetAdjacencyMatrix(smols: list, connected=True):
    max_length = max(mol.GetNumAtoms() for mol in mols)
    bond_labels = [Chem.rdchem.BondType.ZERO] + list(sorted(set(bond.GetBondType() for mol in mols for bond in mol.GetBonds())))
    bond_encoder_m = {l: i for i, l in enumerate(bond_labels)}

    A = np.zeros(shape=(max_length, max_length), dtype=np.int32)
    begin, end = [b.GetBeginAtomIdx() for b in mol.GetBonds()], [b.GetEndAtomIdx() for b in mol.GetBonds()]
    bond_type = [bond_encoder_m[b.GetBondType()] for b in mol.GetBonds()]
    A[begin, end] = bond_type
    A[end, begin] = bond_type
    return A 

GetAdjacencyMatrixで隣接行列を取得します。smilesからChem.MolFromSmilesによりmolオブジェクトにした後、結合種類によってラベル付けをし、分子のインデックス、結合情報が得られるので書き込んでいきます。なお、rdkit.Chem.rdmolops.GetAdjacencyMatrix()だと上記の操作が行えます。

m = Chem.MolFromSmiles('c1ccccc1')
Chem.rdmolops.GetAdjacencyMatrix(m)
  • 特徴量抽出

RDKitの便利な関数により、分子中の各原子ごとの情報を行列にしていきます。RDKitのECFPがどのような情報を用いているかわからなかったので、適当に選択しています。Module Hierarchyが参考になります。

def createNodeFeatures(mol) -> np.ndarray:
    """ノード特徴量を生成する

    :param mol: rdkit mol object
    :return: np.array
    """
    features = np.array([[
        *[a.GetDegree() == i for i in range(5)],
        *[a.GetExplicitValence() == i for i in range(9)],
        *[int(a.GetHybridization()) == i for i in range(1, 7)],
        *[a.GetImplicitValence() == i for i in range(9)],
        a.GetIsAromatic(),
        a.GetNoImplicit(),
        *[a.GetNumExplicitHs() == i for i in range(5)],
        *[a.GetNumImplicitHs() == i for i in range(5)],
        *[a.GetNumRadicalElectrons() == i for i in range(5)],
        a.IsInRing(),
        a.GetAtomicNum(),
        a.GetDegree(),
        a.GetExplicitValence(),
        a.GetImplicitValence(),
        a.GetFormalCharge(),
        a.GetTotalNumHs()] for a in mol.GetAtoms()], dtype=np.int32)

    return features

特徴抽出の結果は、次のようになります。rowがノード数、columnが特徴数です。

>> print(createNodeFeatures(mol))
[[0. 1. 0. ... 0. 3. 1.]
 [0. 0. 1. ... 0. 2. 2.]
 [0. 0. 0. ... 0. 0. 3.]
 ...
 [0. 0. 1. ... 0. 1. 2.]
 [0. 0. 0. ... 0. 0. 3.]
 [0. 1. 0. ... 0. 0. 1.]]
  • MyECFPクラス 上記の隣接行列と原子の特徴行列からECFPを実装します。
class ECFP:
    """Extended Connectivity Fingerpritnsを計算する

    :param mol: rdkit mol object
    :param radius: 隣接する原子情報をどこまで見るか. メッセージパッシングの回数に相当する
    :param nbits: 得られた部分構造の格納数するためのビット数
    :param n_feat: 特徴量行列 (ノード数 × 特徴量数 )
    """

    def __init__(self, mol, radius: int = 2, nbits: int = 2048, n_feat: np.array = None):
        self.mol = mol
        self.radius = radius
        self.nbits = nbits
        self.fps = np.zeros(shape=(self.nbits,), dtype=np.int32)

        if n_feat is None:
            n_feat = createNodeFeatures(mol)

        n_feat = np.array(n_feat, dtype=np.int32)
        n_atoms = n_feat.shape[0]
        self.adj = Chem.GetAdjacencyMatrix(mol)
        deg = np.array(np.sum(self.adj, axis=1), dtype=str)

        # concatenate node features.
        n_feat = np.array([''.join([str(f) for f in n_feat[atom]]) for atom in range(n_atoms)], dtype=str)
        self.n_feat = np.array([n_feat[i] + deg[i] for i in range(n_atoms)])

    def _concat_neighbor(self, atom: int, n_feat: list) -> list:
        """隣接情報を付与する。メッセージパッシングとaggregationに相当する。concatenationで次元を減らす。

        :param atom: 注目している原子のindex
        :param n_feat: 更新前の特徴量情報
        :return: ノード特徴量を返す
        """
        nei_id = np.nonzero(self.adj[atom])[0]
        new = [str(nid) + str(feat) + str(self.adj[atom][nei_id][nid]) for nid, feat in enumerate(n_feat[nei_id])]
        vec = ''.join([str(ind) for ind in new])
        return vec

    def calculate(self, nei_info: bool = False) -> np.ndarray:
        """ECFPを計算する

        :param nei_info: メッセージパッシング前の情報も加えるかどうか。無いほうが精度上がる。
        :return: result of fingerprint.
        """
        n_atoms = self.n_feat.shape[0]
        identifier = copy.deepcopy(self.n_feat)
        for _ in range(0, self.radius):
            for atom in range(n_atoms):
                if _ == 0:
                    v = self.n_feat[atom]
                else:
                    v = self.n_feat[atom] if not nei_info else self._concat_neighbor(atom, identifier)
                identifier[atom] = hash(v)
                index = int(identifier[atom]) % self.nbits
                self.fps[index] = 1

        self.identifier = identifier
        return self.fps

concat_neighbor()では、隣接する原子の情報を集めます。隣接のidと情報、結合次数をconcatenateしています。calculate()でECFPのアルゴリズム通り実装していきます。各原子ごとに隣接する特徴量を取得し、hash()を通し、これをその原子の識別子とします。これを固定長で割った余りが部分構造のindexになります。これを全ての原子に対して行います。メッセージパッシング操作を繰り返すほど、遠くの原子の情報を伝えることに相当します。

NFP実装

NFPに関してはDGL+PyTorchを用いて実装していきます。このライブラリは化学構造に特化していおり、JTNNやMPNNの学習モデルもあり、非常に使いやすいです。DGLはグラフニューラルネットワークを統一的なフレームワークで実装できるようになっています。

前処理

SMILESをグラフデータにする必要があります。簡単なのはmol_tobigraph()を用いることです。ノードやエッジの特徴量に関してはCanonicalAtomFeaturizer()`CanonicalBondFeaturizer()やを用いると自動的に特徴量を作ることができます。

mols_graph = [mol_to_bigraph(
    mol, 
    node_featurizer=CanonicalAtomFeaturizer(), 
    edge_featurizer=CanonicalBondFeaturizer()) for mol in mols]

ECFPと比較するために、上記createNodeFeatures()での特徴量を初期値とします。

モデル

import dgl.function as fn

gcn_msg = fn.copy_src(src='h', out='m')
gcn_reduce = fn.sum(msg='m', out='h')


class NFP(nn.Module):
    def __init__(self, input_dim, hidden_dim=64, depth=2, nbits=16):
        super(NFP, self).__init__()
        self.linear1 = nn.Linear(input_dim, hidden_dim)
        self.linear2 = nn.Linear(hidden_dim, nbits)
        self.softmax = nn.Softmax(dim=1)
        
        self.depth = depth
        self.nbits = nbits
        
        self.linear3 = nn.Linear(nbits, hidden_dim)        
        self.linear4 = nn.Linear(hidden_dim, 1)      
         
    def forward(self, g, n_feat):
        with g.local_scope():
            fps = torch.zeros([1, self.nbits])
            for _ in range(self.depth):
                g.ndata['h'] = n_feat
                g.update_all(gcn_message, gcn_reduce)
                h = g.ndata['h']
                
                r = F.relu(self.linear1(h))
                i = self.softmax(self.linear2(r))
                fps += torch.sum(i, dim=0)

            out = F.relu(self.linear3(fps))
            out = self.linear4(out).squeeze(0)
        return fps, out

モデルは非常にシンプルで、GCNを行った後に線形関数に通し、ソフトマックを通します。この状態だと(バッチ数×ノード数×特徴量)という次元なのでノードに対してsum()を取ってやることで次元を落としてやります。得られたフィンガープリントは、そのままMLPを用いて回帰に用いることができます。学習は、通常のニューラルネットの学習方法と同じです。なお、この実装ではミニバッチに対応していません。DGLではミニバッチを作るときは、グラフ同士の隣接行列を対角行列にし、大きなグラフとすることで対応しています。ミニバッチに対応させるのは少々面倒なので、実装はしていません。また、softmax後にsumを取る必要があるためcudaに対応できていません。

学習に関しては、通常のPyTorchの実装と同じですが、ここではearly stoppingは実装していません。

loss_func = nn.MSELoss(reduction="none")
optimizer = optim.Adam(model.parameters(), lr=0.001)
model.train()
epoch_losses = []
for epoch in range(10):
    epoch_loss = 0
    for ite, batch in enumerate(train_loader):
        _, bg, label, masks = batch
        n_feat = bg.ndata['h']
        fps, prediction = model(bg, n_feat)
        loss = (loss_func(prediction, label) * (masks != 0).float()).mean()
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()
        epoch_loss += loss.detach().item()
    epoch_loss /= (ite + 1)
    print(f'Epoch {epoch}, loss {epoch_loss:.4f}'
    epoch_losses.append(epoch_loss)

実験条件

zincデータ22万のうち教師データ1800, テストデータ200としています。logP(溶解性の指標)に対して

  1. NFP + MLP
  2. NFP + Random Forest
  3. 特徴量初期値のノードに対する和 + Random Forest
  4. ECFP + Random Forest
  5. RDkit Morgan Fingerprint + Random Forest

で比較を行っています。NFPの条件は、radius=2, nbits = 16です。MLPは2層の隠れ層のノード数は64で、エポック数は10です。Random Forestはscikit-learnのデフォルトです。ECFP, RDKit Morgan fingerprint はともにradius=2, nbits=2048です。

結果

テストデータに対する予測精度のスコアを載せます。

model r2 RMSE
1 0.8945 0.3577
2 0.8902 0.3725
3 0.7827 0.7371
4 0.7505 0.8463
5 0.7346 0.9002

NFPは学習させている分、予測精度が上がっています。分子グラフの特徴量をうまく埋め込むことができているようです。途中で取り出した16次元のベクトルをRandom Forestで回帰しても精度は同じくらい出ていますが、データ数を増やしてみるとMLPの方が精度が高めに出ることの方が多いです。ECFPの精度は自作の方が良かったのですが、データセットに依存します。また、radiusを増加させていくと精度が落ちていきます。

実際には、エッジ情報を用いた方が精度がもう少し向上します。MPNNやSchnetはDGLでAlchemyの学習済みモデルもあるので、転移学習させても面白いかもしれません。

最後に

GNNを用いることは多いのですが、実際のタスクに適用すると大抵の場合は求めたい精度に達しません。これはグラフのノード・エッジ特徴量の初期値に依存し、そもそもこれらの特徴量がある程度、目的変数に対して相関がないと、回帰・識別タスクに適用しても悲惨なことになることが多いです。やはり楽な方法はなく、特徴量エンジニアリング、第一原理計算やMD計算など原理原則に基づいたパラメータの探索と色々試行錯誤が必要となってきます。低分子に対するGNNの研究はもうやりつくされた感があります。ポリマー、タンパク質など巨大な分子はDFT計算でも計算が難しいので、GNNをどう適用するかが大事になってきそうです。とりあえずGNN+NASや分散学習を試していきたい。

ソースコードはこちら。

github.com

参考文献

  • ChemAxon ECFP
  • RDKit Fingerprints Document
  • 化合物でもDeep Learningがしたい

Chemical VAE

一年くらい前に書いた記事の実装を行ってみました。

udnp.hatenablog.com

コード:

github.com

正しいものを使用したいときは本家本元から参照してください。

簡単に解説

分子構造はSMILESという文字列で表すことができます。例えば、ヒドロキシ基がついたベンゼン環(フェノール)はOc1ccccc1で表すことができます。文字列ならば、生成モデルを使ってOc1ccccc1→潜在空間→Oc1ccccc1となるようなモデルを作ってしまえば、多様な分子を作り出すことができるというわけです。目的の活性や物性をもつような潜在変数を探してやって、そこからデコードしてやれば目的の分子を作り出すことができるという考えです。

実装するためのポイント

前処理

まずはSMILESを機械が読み込めるような形にしていく必要があります。

化合物DBであるZINCからSMILES中の文字を数字にマッピングするための辞書を作成します。全SMILESから重複しないように文字をリストに格納していきます。

def smiles2one_hot_chars(smi_list: list) -> list:
    """obtain character in SMILES.

    :param smi_list: SMILES list
    :return: Char list
    """
    char_lists = [list(smi) for smi in smi_list]
    chars = list(set([char for sub_list in char_lists for char in sub_list]))
    chars.append(' ')

    return chars


df = pd.read_table('./data/train.csv',  header=None)
df.columns = ['smiles']
smiles = list(df.smiles)

zinc_list = smiles2one_hot_chars(smiles)

>
[
    "7", "6", "o", "]", "3", "s", "(", "-", "S", "/", "B", "4", "[", ")", "#", "I", "l", "O", "H", "c", "1",
    "@", "=", "n", "P", "8", "C", "2", "F", "5", "r", "N", "+", "\\", " "
]

次に、文字をマッピングするための辞書を作ります。

char2id = dict((c, i) for i, c in enumerate(zinc_list))
id2char = dict((i, c) for i, c in enumerate(zinc_list))

これでSMILESを文字からid, idから文字に変換することができます。ニューラルネットに学習できるように(データ数×最大原子数×辞書数)となるように変形していきます。ZINCは低分子DBなので、最大原子数は120、辞書数は35となります。Chemical VAEでは、文字列を同じにするように空白でパディングしています。

def pad_smile(string: str, max_len: int, padding: str = 'right') -> str:
    if len(string) <= max_len:
        if padding == 'right':
            return string + " " * (max_len - len(string))
        elif padding == 'left':
            return " " * (max_len - len(string)) + string
        elif padding == 'none':
            return string


def smiles_to_hot(smiles: list, max_len: int = 120, padding: str = 'right',
                  char_to_id: Dict[str, int] = char2id) -> np.ndarray:
    """smiles list into one-hot tensors.

    :param smiles: SMILES list
    :param max_len: max length of the number of atoms
    :param padding: types of padding: ('right' or 'left' or None)
    :param char_indices: dictionary of SMILES characters
    :return: one-hot matrix (Batch × MAX_LEN × len(dict))
    """
    smiles = [pad_smile(smi, max_len, padding) for smi in smiles]
    hot_x = np.zeros((len(smiles), max_len, len(char_to_id)), dtype=np.float32)
    for list_id, smile in enumerate(smiles):
        for char_id, char in enumerate(smile):
            try:
                hot_x[list_id, char_id, char_to_id[char]] = 1
            except KeyError as e:
                print("ERROR: Check chars file. Bad SMILES:", smile)
                raise e
    return hot_x

idがあれば1を追加していくという具合で配列を作成してきます。これは文章生成で行われる前処理とほとんど同じだと思います。

phenol = 'Oc1ccccc1
vec = smiles_to_hot([phenol])
vec

>array([[[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 1.],
        [0., 0., 0., ..., 0., 0., 1.],
        [0., 0., 0., ..., 0., 0., 1.]]], dtype=float32)

このような(1×120×35)のテンソルが得られます。ここでの目標は、この値を入力して、元に戻るように学習していくことです。SMILESへの復元は対応する番号へのマッピング操作を行ってやればよいです。

def hot_to_smiles(hot_x: np.ndarray, id2char: Dict[int, str] = id2char) -> list:
    """one hot list to SMILES list.

    :param hot_x: smiles one hot (id, max_len, node_dict)
    :param id2char: map from node id to smiles char
    :return: smiles list
    """
    smiles = ["".join([id2char[np.argmax(j)] for j in x]) for x in hot_x]
    smiles = [re.sub(' ', '', smi) for smi in smiles]  # paddingを消す
    return smiles


hot_to_smiles(vec[0])
> ['Oc1ccccc1']

学習

VAEではEncoder、Decoderを用いたモデルですが、エンコーダで得られるmu、sigmaから正規分布によるサンプリング結果をデコーダーに渡すことがAEとの違いになります。VAEの学習の数理に関しては前回の記事を参照してください。簡単に書くとELBOを最大にすることで、対数尤度を最大にするように学習を進めていけばいいです。式変形を行っていくと、損失関数は、クロスエントロピー + カルバックライブラーダイバージェンスを計算すればよいという形になります。VAEの実装に関しては様々なレポジトリがあるので、そちらを参考に実装しています。PyTorchには便利な関数があるので、

def loss_function(recon_x, x, mu, logvar):
    BCE = F.binary_cross_entropy(recon_x, x, reduction='sum')
    KLD = -0.5 * torch.mean(1. + logvar - mu ** 2. - torch.exp(logvar))
    return BCE + KLD
  • Encoder エンコーダでは、1d convolutionを行っています。

  • Decoder デコーダーでは、サンプリングされたあと、GRUに渡せるように次元を(バッチ数×最大原子数×辞書数)となるようにRepeatさせています。その後は、線形変換→ソフトマックスという形で入力と同じ形になるようにしています。

サンプルコードではGANも実装しています。モデルは全く同じです。

結果

データ数、学習時間が少ないため、ほとんど元の分子に戻っていません。ネットワークも正しく決めないとダメなのかもしれません。ConvやGRU増やしたり、学習エポック数を増やしたり。GANほうが割と分子っぽい結果は出ています。 (実装も間違ってるかも)

所感

SMILESは化学構造を一意に表すことができません。最初に記載したフェノールはc1ccccc1Oやc1c(O)cccc1と書いても良いので、組み合わせ爆発を起きてしまいます。 改良法には文法規則を組み合わせたgrammar VAEがありますが、低分子ならJT-VAEやInventを使えば問題ないと思います。

chem VAEは分子生成には使うのは厳しいですが、特徴量抽出には使えるかもしれません。

Junction Treeアルゴリズム

ベイジアンネットワークの推論では良く用いられるアルゴリズムです。分子構造もグラフ構造なので本アルゴリズムは適用することができます。グラフを木分解(tree decomposition)することで得られる構造がjunction treeです。有向グラフの場合、無向グラフにする操作(モラル化など)が必要となりますが、分子構造の場合は無向グラフなので、この変換は必要ありません。

参考にしてる文献は下記です。

arxiv.org

JT-VAEを使ってみた結果など、こちらも参考になります。

qiita.com

グラフ理論

改めてグラフ理論の定義をおさらいしてみます。

グラフ Gは頂点もしくはノードの有限非空集合 (finite non-empty set)  Vと頂点Vの非順序ペア (unordered pairs of vertices)からなるエッジの集合 E \subseteq  V \times V からなるデータ構造です。

\begin{align} G=(V, E) \end{align}

分子構造だと、ノードは原子(C, N, O, F)などでエッジが結合(単結合、二重結合、・・・)などです。データ同士に関係性を持つようなデータセットはグラフ構造で表すことができます。

グラフのデータ構造は

  • 隣接リスト

  • 隣接行列

形式の2通りで表現できます。隣接行列の方がわかりやすいのですが、スパースになるため大きい分子構造には対応が難しいかもしれません。ただ、GNNと生成モデルを組み合わせればグラフ生成をワンショットで行えるので面白いです。隣接リストの方がプログラミングする上では簡単です。どちらが良いかは状況によるので使い分けるの吉です。

サブグラフとは

グラフ

\begin{align} H=(V_H, E_H) \end{align}

を考えます。 V_H \in V, E_H \in E ならば  G のサブグラフといいます。またG Hのスーパーグラフともいえます。頂点のサブセット V' \subseteq Vが与えられれば、他のサブグラフ G' = (V',E') V中の頂点間 G中に存在する全てのグラフから完全に構成されるので、より厳密に全ての頂点は

\begin{align} v_i, v_j \in V', \\ (v_i, v_j) \in E ⇔ (v_j, v_i) \in E \end{align}

です。言い換えれば、2ノードは G' に隣接し、Gに隣接してさえすればよいです。ノードの全ペア間のエッジが存在するならば、(サブ)グラフは完全(もしくはクリーク)と呼ばれます。このクリークが他のクリークに含まれないならばそれを極大クリーク (maximal clique )と呼びます。

Junction Tree(JT

木分解によりグラフ Gをjunction treeに写像できます。この操作によりサイクルフリーな構造(木構造)になります。 JT \mathcal{T}_G = (\mathcal{V, E, X})はノード集合が \mathcal{V} = \{C_1, \ldots, C_n\} でエッジ集合が \mathcal{E}と表せます。 \mathcal{X}はサブグラフ構造のラベルです。各ノード、すなわちクラスタ C_i = (V_i, E_i)となります。正確な定義は次の通りになります。

(1) 全クラスターのユニオンはGと等しい。すなわち

\begin{align} \cup_i V_i = V \text{かつ} \cup_i E_i = E \end{align}

(2) インターセクションがある。全クラスター間で以下となる。

\begin{align} V_i \cap V_j \subseteq V_k \end{align}

ジャンクションツリーにすることで、分子構造を生成する上で問題であった「環構造の発生」が問題なく行えるようになります。たとえばサブグラフのベンゼンが1つのクラスター(ノード)としてあらわすことができるようになります。

簡易的なコード

f:id:udnp:20191014120931p:plain

この化合物を木分解してみます。

https://github.com/masahiro-mochizuki/pd1_inhibitor_dataset

クリークの作成

まずはRDkitの関数によりSSSR (smallest set of smallest rings)によりリングの個数を計算します。とりあえず、上記のpd1_inihibitor_datasetからsmilesを取得しておきます。

import sys
import rdkit
from rdkit import Chem


def get_clique_mol(mol, atoms):
    smiles = Chem.MolFragmentToSmiles(mol, atoms, kekuleSmiles=True)
    new_mol = Chem.MolFromSmiles(smiles, sanitize=False)
    new_mol = copy_edit_mol(new_mol).GetMol()
    new_mol = sanitize(new_mol)
    return new_mol


def get_cliques(mol):
    cliques = []
    for bond in mol.GetBonds():
        a1 = bond.GetBeginAtom().GetIdx()
        a2 = bond.GetEndAtom().GetIdx()
        if not bond.IsInRing():
            cliques.append([a1,a2])
    return cliques

path = "data/pd1_inhibitor_dataset/"
df = pd.read_csv(path+"PD1_inhibitor_dataset.csv")

mol = Chem.MolFromSmiles(df.smiles[0])
n_atoms = mol.GetNumAtoms()
cliques = get_cliques(mol)
ssr = [list(x) for x in Chem.GetSymmSSSR(mol)]
cliques.extend(ssr)
print(cliques)
for i in range(len(ssr)):
    print('{}番目の環のサイズ:{}'.format(i, len(ssr[i])))

# 隣接リスト
nei_list = [[] for i in range(n_atoms)]
for i in range(len(cliques)):
    for atom in cliques[i]:
        nei_list[atom].append(i)

print(nei_list)

上記のコードを実行すると、クリークとなる番号が取得できます。

0番目の環のサイズ:45
1番目の環のサイズ:5
2番目の環のサイズ:5
3番目の環のサイズ:6
4番目の環のサイズ:5
5番目の環のサイズ:6
6番目の環のサイズ:5
7番目の環のサイズ:6

[[0, 1], [1, 2], [2, 3], [3, 4], [5, 6], [12, 13], [15, 16], [16, 17], [17, 18], [17, 19], [20, 21], [23, 24], [24, 25], [24, 26], [27, 28], [29, 30], [31, 32], [32, 33], [39, 40], [42, 43], [43, 44], [44, 45], [44, 46], [47, 48], [55, 56], [58, 59], [59, 60], [69, 70], [72, 73], [73, 74], [77, 78], [81, 82], [84, 85], [85, 86], [86, 87], [86, 88], [89, 90], [92, 93], [93, 94], [93, 95], [95, 96], [96, 97], [97, 98], [97, 99], [103, 104], [106, 107], [107, 108], [110, 111], [112, 113], [114, 115], [117, 118], [119, 120], [121, 122], [122, 123], [129, 130], [131, 132]]

スパニングツリー

スパニングツリーは、あるグラフの全ての頂点とそのグラフを構成する辺の一部分のみで構成される木です。クリークが得られれば、ジャンクションツリーにするのは簡単です。スパニングツリーにするためのアルゴリズムを用います。

もともとは最小経路を求めるアルゴリズムです。

有名どころだと

があります。プリム法はダイクストラ法とほぼ同じ方法なので実装は簡単です。ポイントはプライオリティキューを用いることです。最初の経路を選択していき、通った経路をヒープキューにプッシュしていきます。scipyにはクラスカル法が実装されています。クラスカル法は適当なエッジを選択していきますが、1つノードに3つ以上のエッジが接続したとき、最小となる重みを選択していきます。英語版のwikiをみるとわかりやすいです。

プリム法の例

import heapq
 
 
# minimum spanning tree
# graph[vertex] = (weight, end, start)
def prim(graph):
    n = len(graph)
    used = [True] * n
    edgelist = []
    for edge in graph[0]:
        heapq.heappush(edgelist, edge)
 
    used[0] = False
    res_edge = []
    res = 0
    while edgelist:
        minedge = heapq.heappop(edgelist)
        if not used[minedge[1]]:
            continue
 
        v = minedge[1]
        used[v] = False
        for edge in graph[v]:
            if used[edge[1]]:
                heapq.heappush(edgelist, edge)
        res += minedge[0]
        res_edge.append((minedge[2], minedge[1])
 
    return res, res_edge

※ 実際はscipy.sparse.csgraphminimum_spanning_tree()を使えばよいです。

JTを作成する場合の注意点はMinimum Spanning TreeではなくMaximal Spanning Treeの点です。最大の重みを定義しておき、引いておく(逆の操作)で最大のスパニングツリーを取得します。

木分解

from collections import defaultdict
from scipy.sparse import csr_matrix
from scipy.sparse.csgraph import minimum_spanning_tree
from collections import defaultdict
from rdkit.Chem.EnumerateStereoisomers import EnumerateStereoisomers, StereoEnumerationOptions


def tree_decomp(mol):
    n_atoms = mol.GetNumAtoms()
    cliques = get_cliques(mol)
    ssr = [list(x) for x in Chem.GetSymmSSSR(mol)]
    cliques.extend(ssr)

    nei_list = [[] for i in range(n_atoms)]
    for i in range(len(cliques)):
        for atom in cliques[i]:
            nei_list[atom].append(i)
    
    #Merge Rings with intersection > 2 atoms
    for i in range(len(cliques)):
        if len(cliques[i]) <= 2: 
            continue
        for atom in cliques[i]:
            for j in nei_list[atom]:
                if i >= j or len(cliques[j]) <= 2: 
                    continue
                inter = set(cliques[i]) & set(cliques[j])
                if len(inter) > 2:
                    cliques[i].extend(cliques[j])
                    cliques[i] = list(set(cliques[i]))
                    cliques[j] = []
    
    cliques = [c for c in cliques if len(c) > 0]
    nei_list = [[] for i in range(n_atoms)]
    for i in range(len(cliques)):
        for atom in cliques[i]:
            nei_list[atom].append(i)
    
    #Build edges and add singleton cliques
    edges = defaultdict(int)
    for atom in range(n_atoms):
        if len(nei_list[atom]) <= 1: 
            continue
        cnei = nei_list[atom]
        bonds = [c for c in cnei if len(cliques[c]) == 2]
        rings = [c for c in cnei if len(cliques[c]) > 4]
        if len(bonds) > 2 or (len(bonds) == 2 and len(cnei) > 2): 
      
            cliques.append([atom])
            c2 = len(cliques) - 1
            for c1 in cnei:
                edges[(c1,c2)] = 1
        elif len(rings) > 2: #Multiple (n>2) complex rings
            cliques.append([atom])
            c2 = len(cliques) - 1
            for c1 in cnei:
                edges[(c1,c2)] = MST_MAX_WEIGHT - 1
        else:
            for i in range(len(cnei)):
                for j in range(i + 1, len(cnei)):
                    c1,c2 = cnei[i],cnei[j]
                    inter = set(cliques[c1]) & set(cliques[c2])
                    if edges[(c1,c2)] < len(inter):
                        edges[(c1,c2)] = len(inter) #cnei[i] < cnei[j] by construction

    edges = [u + (MST_MAX_WEIGHT-v,) for u,v in edges.items()]
    
    if len(edges) == 0:
        return cliques, edges

    #Compute Maximum Spanning Tree
    row, col, data = zip(*edges)
    n_clique = len(cliques)
    clique_graph = csr_matrix( (data,(row,col)), shape=(n_clique,n_clique) )
    junc_tree = minimum_spanning_tree(clique_graph)
    row,col = junc_tree.nonzero()
    edges = [(row[i],col[i]) for i in range(len(row))]
    return (cliques, edges)

結果の可視化

import networkx as nx
import matplotlib.pyplot as plt
plt.figure(figsize=(8, 6), dpi= 80, facecolor='w', edgecolor='k')

graph = nx.Graph()
graph.add_edges_from(edges)
import matplotlib.pyplot as plt
nx.draw_networkx(graph)
plt.show()

networkxで可視化すると次のような図が得られました。各ノードにはそれぞれサブグラフであるので、分子を生成するためには、木構造を拡大していけばよいです。

f:id:udnp:20191014123920p:plain

Next Action?

木構造に対してニューラルネットワークを行っていく ⇒ JT-VAEのアプローチです。

  • リーフを拡張するのにRNN + MCTS
  • そもそもリーフをでなくて強化学習でノードを拡張していく。

複数のアプローチがまだまだ考えられますね。

JT-VAEを試してみると、分子生成部分にはグラフニューラルネットワークにあまり依存しない感じなのと、復元率が論文程よくないため、まだまだ改善点がありそうですね。最適化の工夫もできそうです。時間があれば改良していきたい。他の論文もウォッチしなければ。DGLも試したい・・・。

tf-idfを実装する

Qiitaにもupしています。

qiita.com

wikipediaから引用すると、

tf-idfは文章中に含まれる単語の重要度を評価する手法の1つであり、主に情報検索やトピック分析などの分野で用いられている。

と記載されています。文書中にどの単語が重要かを測定するために

  • TF (Term Frequency )
  • IDF (Inverse Document Frequenc)

を掛けることでその指標としています。デジタルライブラリの文章レコメンドシステムの83%がtf-idfを使用しているそうです。

ここでは、tf-idfの考え方と実際にコーディングをすることで理解を深めたいと思います。

目次

Term frequency (TF)

Term fequency (単語の出現頻度)は、文字通り文書中に単語がでてくる頻度です。tfだけでも重み付けにより複数の方法が提案されています。最も単純な選択はドキュメント中のタームのカウントです。ドキュメント、タームをそれぞれ d, tとします。

  • バイナリ

ブール型の頻度「frequencies」です。文章中に単語が存在するなら1にそれ以外なら0です。これは単純なfingerprintに対応します。

\begin{align} \text{tf}(t, d) = 1 \quad \text{if t occurs in d else 0 } \end{align}

  • カウント (raw count)

文章中に単語が存在したとき単純にカウントしていく方法です。sklearnのCountVectorizer()を使えば簡単に計算できます。

\begin{align} \text{tf}(t, d) = f_{t, d} \end{align}

  • 単語の出現頻度 (term frequency)

これがいわゆるtfです。上記のカウントを文章中のカウントの総和で割り出現頻度を求める方法です。

\begin{align} \text{tf}(t, d) = \frac{f_{t, d}}{\text{number of words in d}} = \frac{f_{t, d}}{\sum_{t' \in d} f_{t', d}} \end{align}

他の方法では以下のようなものが提案されています。

  • 対数正規化 (log normalization)

\begin{align} \text{tf}(t, d) = \log{(1 + f_{t, d})} \end{align}

  • 二重K正規化 (double normalization K)

K = 0.5とした場合が使用されるようです。

\begin{align} \text{tf}(t, d) = K + (1 - K) \cdot \frac{f_{t,d}}{\max\{f_{t', d} :t' \in d\} } \end{align}

Inverse document frequency (IDF)

inverse document frequency (逆文書頻度)は、単語が与える情報がどれくらいを測る指標です。よく出てくるもの値は下がり、レアなものに値が大きくなるように重みづけされます。一般語フィルタとして働きます。

\begin{align} \text{idf}(t, D) = \log \frac{N}{|\{ d \in D:t \in d\}|} \end{align}

ここで、Nはコーパス中の文書数です。分母はdf (document frequency)であり、文書全体の単語の出現頻度です(tfは文章中の出現頻度)。

  • IDF

良く紹介されるものは次の式

\begin{align} \text{idf}(t, D) = \log \frac{N}{df_{t}} \end{align}

です。

  • smooth IDF

sklearnのTfidfVectorizerではデフォルトでsmooth_idf = Trueと平滑化しており、

\begin{align} \text{idf}(t, D) = \log (\frac{1 + N}{1+ df_t}) + 1 \end{align}

となっています。他の方法として

  • idf-max

\begin{align} \text{idf}(t, D) = \log (\frac{\max_{t' \in d} df_t'}{1+ df_t}) \end{align}

  • probabilistic idf

\begin{align} \text{idf}(t, D) = \log \frac{N - df_t}{df_t} \end{align}

があります。ドキュメントの頻度が小さいときはほとんどIDFの値は変わりませんが、頻度が大きくなるとIDFの値はsmooth IDF > IDF > proba. IDF の順になります。proba. IDFだとdf = 50からIDFの値は負になります。

TF-IDF

上記の二つの指標を掛け合わせたものです。

\begin{align} \text{tfidf} (t, d, D) = \text{tf} (t,d) \cdot \text{idf}(t, D)
\end{align}

複数の組み合わせがあります。document term weightquery term weightがあります。

手法 dtw qtw
1  f_{t, d} \cdot \log \frac{N}{df_t}  \left(0.5 + 0.5 \cdot \frac{f_{t,d}}{\max\{f_{t', d} :t' \in d\} } \right) \cdot  \log \frac{N}{df_t}
2  1 + \log f_{t, d}  \log (1 + \frac{N}{df_t})
3  (1 + \log f_{t, d}) \cdot \log \frac{N}{df_t}  (1 + \log f_{t, d}) \cdot \log \frac{N}{df_t}

コード

例文としてsklearnにある文章を用います。

corpus = [
'This is the first document.',
'This document is the second document.',
'And this is the third one.',
'Is this the first document?']

sklearnを使う場合

こちらは非常に簡単で3行でtfidfの計算ができます。

from sklearn.feature_extraction.text import TfidfVectorizer

tfidf = TfidfVectorizer()
x = tfidf.fit_transform(corpus)

tfidfの結果をみるために、データフレームに変換します。

import pandas as pd

df_tfidf = pd.DataFrame(x.toarray(), columns=tfidf.get_feature_names())
print(df_tfidf)

自作した場合

sklearnのものと同じになるか試してみます。corpus内の単語を数えるためにCountVectorizerを用います。

from sklearn.feature_extraction.text import CountVectorizer
from sklearn.preprocessing import normalize
import numpy as np

smooth_idf = True
norm_idf = True

wc = CountVectorizer()
x = wc.fit_transform(corpus)
wcX = np.array(x.toarray())

# term frequency:
N = wcX.shape[0]
tf = np.array([wcX[i, :] / np.sum(wcX, axis=1)[i] for i in range(N)])

# inverse documents frequency
df = np.count_nonzero(wcX, axis=0)
idf = np.log((1 + N) / (1 + df)) + 1  if smooth_idf else np.log( N / df )  

# normalize
tfidf = normalize(tf*idf) if norm_idf else tf*idf
tfidf = pd.DataFrame(tfidf, columns=wc.get_feature_names())

※ 正規化しないと同じ結果になりません。

print(tfidf)
        and  document     first  ...       the     third      this
0  0.000000  0.469791  0.580286  ...  0.384085  0.000000  0.384085
1  0.000000  0.687624  0.000000  ...  0.281089  0.000000  0.281089
2  0.511849  0.000000  0.000000  ...  0.267104  0.511849  0.267104
3  0.000000  0.469791  0.580286  ...  0.384085  0.000000  0.384085

CountVectorizerも使わない場合

本来はscipy.sparse.csr_matrixを使って処理すべきです。

import re
from collections import defaultdict

documents = [re.sub('[.|?]', '', i.lower()) for i in corpus]
documents = [doc.split(' ') for doc in documents]

vocab = defaultdict()
vocab.default_factory = vocab.__len__
for doc in documents:
    feature_counter = {}
    for feature in doc:
        feature_idx = vocab[feature]
        if feature_idx not in feature_counter:
            feature_counter[feature_idx] = 1
        else:
            feature_counter[feature_idx] += 1


sorted_feature = sorted(vocab.items())
for new_val, term in enumerate(sorted_feature):
    vocab[term] = new_val

X = np.zeros(shape=(len(corpus), len(sorted_feature)), dtype=int)
for idx, doc in enumerate(documents):
    for word in doc:
        if word in vocab.keys():
            X[idx, vocab[word]] += 1

これでCountVectorizer()と同じ結果になります。

最後に

tf-idfの計算だけならsklearnを用いれば非常に簡単に行えることがわかります。テキストの前処理はnltkWordNetLemmatizerがありこちらも必須です。また特徴量を作るのにナイーブベイズも良く用いられますね。

前処理としては、tf-idfは使っていきたいとも思います。目標は化学言語のSMILESを生成させることですので、tf-idfは用途が違うかもしれません。SMILESでは文章のつながりが重要ですので、文章中の単語を予測するために用いられるword2vecBERTなどもこちらの方が良さそうな気がしています。こちらと生成モデルのVAE (Seq2Seq)も試していきたいところです。

参考文献

ベイズ最適化 (その3):ハイパーパラメータ探索

前回の続きです。

udnp.hatenablog.com

目次

ベイズ最適化が良く用いられる例だと思います。他にはハイパーパラメータの探索にはグリッドサーチや最適化計算を行う方法もあります。ハイパーパラメータの束縛条件を設定し、とりあえず試すだけなので簡単です。

xgboostのハイパーパラメータを下記手法

  • Randomized Search
  • Grid Search
  • Bayesian Optimization
  • Nevergrad

で探索し、結果を比較してみようと思います。データセットはsklearn.datasetsにあるdiabetesです。結論から言うとPSOが最も性能が良かったです。場合によるのでしょうが。

全体のソースコードは以下です。

github.com

baseline

まずは探索無しのデフォルトの場合です。

from sklearn import datasets
from sklearn.model _selection import GridSearchCV, cross_val_score
from sklearn.preprocessing import StandardScaler
from xgboost import XGBRegressor 
import numpy as np 
import matplotlib.pyplot as plt

X,y = datasets.load_diabetes(return_X_y=True)

scaler = StandardScaler()
xgb = XGBRegressor()
pipe = Pipeline([('scaler': scaler), ('xgb', xgb)])
 
scoring = 'neg_mean_squared_error'
baseline = cross_val_score(pipe, X, y, scoring=score, cv=5)

必要なパッケージを呼び出して、パイプラインに通します。前処理は特に必要ないのですが、正規化しています。 デフォルトの結果を示します。

print(baseline.mean())
> -3435.44

Algorithm

Document

適当に探索する方法です。

param = { 'xgb__gamma' range(0, 5), 'xgb__max_depth': range(1, 50), 'xgb__min_child_weight': range(1, 10), 'xgb__n_estimators' : range(0, 300)
rs = RandomizedSearchCV(pipeline, param_distributions=param, scoring=scoring, n_iter=25)

Document

ハイパーパラメータをグリッド上にして網羅的に計算しますが、ハイパーパラメータが多いと時間がかかるのとドメイン知識が必要となります。

param = { 'xgb__gamma': [0.01, 0.1] , 'xgb__max_depth': [1, 3, 5],  'xgb__min_child_weight': [1, 3, 10],  'xgb__n_estimators' : [50, 100, 500]}
gs = GridSearchCV(pipe, param_distributions=param, cv=5, scoring=scoring, iid=False)

GPyOpt

ベイズ最適化を用いてハイパーパラメータを探索します。GPでは、自作kernelも渡すことができます。獲得関数も複数から選択できます。 目的関数ですが、交差検証でのスコア値が最大になるようにベイズ最適化を行っていきます。

import GPy
import GPyOpt
from GPyOpt.methods import BayesianOptimization

bounds = [
    {'name': 'learning_rate', 'type': 'continuous', 'domain': (0, 0.1)}, 
    {'name': 'gamma', 'type': 'continuous', 'domain': (0, 5)}, 
    {'name': 'max_depth', 'type': 'discrete', 'domain': (1, 50)}, 
    {'name': 'n_estimators', 'type': 'discrete', 'domain': (1, 500)}, 
    {'name': 'min_child_weight', 'type': 'discrete', 'domain': (1, 10)}, 
]

def cv_score(*args):
    args = args[0]
    score = cross_val_score(
        XGBRegressor(learning_rate=args[0],
                                gamma=args[1],
                                max_depth = int(args[2]),
                                n_estimators=int(args[3]),
                                min_child_weight=int(args[4)),
        X, y, scoring=score, cv=5).mean()
        return score        

optimizer = BayesianOptimization(f=cv_score, domain=bounds, model_type='GP', acquisition_type='EI')
optimizer.run_optimization(max_iter=20)

Nevergrad

勾配を必要としない最適化アルゴリズムが多数実行できるようなパッケージになります。現時点で使えるアルゴリムは 100個ほどあり下記で確認できます。

from nevergrad.optimization import registry
print(sorted(registry.keys()))

BOも入っているようで、nervergradでもベイズ最適化できます。githubを見るとTwoPointsDEが良いみたいです。worker数の選択が大事になりそうです。まずはnevergrad用に目的関数を設計します。

def cv_score2(*args):
    lr, gamma, max_depth, n_est, min_child = args
    score = cross_val_score(
        XGBRegressor(learning_rate=lr[0],
                                gamma=gamma[0],
                                max_depth = int(max_depth),
                                n_estimators=int(n_est),
                                min_child_weight=int(min_child)),
        X, y, scoring=score, cv=5).mean()
        return -score                        

各ハイパーパラメータごとにサンプリング分布の設定や境界条件などを設定し最適化計算を行います。

from nevergrad.optimization import optimizerlib
from nevergrad import instrumentation as inst

# log distributed between 0.001 and 1
lr = inst.var.Array(1).bounded(0, 3).exponentiated(base=10, coeff=-1)
gamma = inst.var.Array(1).bounded(0, 3)
max_depth = inst.var.OrderedDiscrete(range(1, 50))
n_estimators = inst.var.OrderedDiscrete(range(1, 500))
min_child = inst.var.OrderedDisrete(range(1, 10))

instrumentation = inst.Instrumentation(lr, gamma, max_depth, n_estimators, min_child, num_workers=5)
optimizer = optimizerlib.TwoPointsDE(instrumentation, budget = 100)
args, kwargs = instrumentation.data_to_argument([0] * instrumentation.dimension)

for _ in range(optimizer.budget):
    x = optimizer.ask()
    y = cv_score2(*x.args)
    optimizer.tell(x, y)

recommendation = optimizer.provide_recommendation()
print(recommendation.args, cv_score2(*recommendation.args))

budgetがイテレーション数です。askで計算をし、tellで結果を返し、stepを更新していきます。

色々な最適化計算ができるのですが、ほかの手法を使うときは

optimizer = optimizerlib.PSO(instrumentation, budget = 100)

となります。

結果

baselineではneg MSEは -3435.44です。CV=5, イテレーションは30です。Grid searchは比較できないですが参考に載せています。

Algorithm time (sec) neg-MSE
Random 47.98 -3241.93
Grid 74.47 -3193.28
BO 65.08 -3218.64
RS 52.64 -3284.53
2ptsDE 47.62 -3338.47
FastGA1+1 35.08 -3196.43
PSO 56.65 -3161.58

基本的にベースラインは超えるので、ハイパーパラメータサーチは行った方がいいですが、手法間でそこまで大差があるとは思えませんね。

所感

多数手法があるので適当にアルゴリズムを選択しましたが、PSOが最もよかったです。最適化計算は大体PSO使ったときが性能が高い気がします。それにしても物凄く簡単にハイパーパラメータ探索ができました。分子生成の最適化計算時にも試してみたいです。

最適化計算の問題は、やはり計算量が膨大となることですね。TorqueなどMPIは必須でしょう。 アニーラーやゲート型コンピューターが発展してきたとき、どうなるか楽しみです。

引用・参考元

ベイズ最適化のところは、主に下記のブログから引用しています。

krasserm.github.io

github.com

ベイズ最適化 (その2):獲得関数

前回の記事ではベイズ最適化で使用されるガウス過程回帰(Gaussian Process Regression)についてまとめていきました。今回の記事では、ガウス過程を用いたベイズ最適化について行っていきたいと思います。よく用いられるのは、ハイパーパラメータの探索や、実験点の探索などです。利用と探索を考慮した最適化手法になります。

github.com

スターが多かったのでBaysianOptimizationのパッケージを使って試してみたいと思います。ベイズ最適化ではないですが、facebook reserachの最適化ライブラリNervergradも気になっています。

理論について下記を参考にしています。

[1807.02811] A Tutorial on Bayesian Optimization

以下の記事も大変よくまとまっております。ベイズ最適化の例では、主にこちらのコードを使わせていただきました。

krasserm.github.io

ベイズ最適化

ベイズ最適化は、機械学習ベースの最適化法の一種です。数式で記述すると次の問題を解くことと同一です。

\begin{align} x^* = \arg \max_{x \in A} f(x) \end{align}

すなわち、ある目的関数 f(x)を最大にするようなxを求めていくことになります。

ベイズ最適化を行う前提条件として、以下のことを想定しています。

  • 入力xは \mathcal{R}^dが多すぎない。 d \leq 20が成功しやすい。
  • パラメータセットAはシンプルである。Aはhyper rectangle  \{x \in \mathcal{R}^d : a_i \leq x_i \leq b_i\}もしくはd次元のsimplex  \{x \in  \mathcal{R}^d : \sum_i x_i = 1\}が良い。
  • 目的関数 fは連続である。ガウス過程回帰が良く用いられる。
  •  fは"評価することが困難"である。評価自体が難しい or 時間・コストがかかるなど。
  •  fブラックボックスである。凸や線形性などもわからない。
  •  fを評価するとき、f(x)だけあり、1次導関数や2次導関数がない。導関数があるならニュートン法準ニュートン法を用いられる。
  • ノイズは各評価に対して独立。
  • 局所解ではなく大域的最適化を見つける。

ブラックボックス微分不要な関数の最適化ならばベイズ最適化がよく用いられるようになりました。特にニューラルネットワークのハイパーパラメータの調整(Snoek et al, 2012)です。また、AlphaGoでもベイズ最適化が用いられています(Yutian et al., arxiv)。

アルゴリズム概要

f:id:udnp:20190420195143p:plain

アルゴリズムの考え方自体は、非常にわかりやすいです。ガウシアンプロセスなどを用いて値を観察します。すべてのデータを用いて事後確率分布を更新し、このときの分布から獲得関数が最大値となる x_nを求めた後に、 y_n = f(x_n)を求めます。この操作を繰り返していきます。 f(x)は平均 \mu_n (x)や分散 \sigma^{2}_n (x)の分布が得られます。

ガウス過程回帰

ガウス過程の理論は前回の記事

ベイズ最適化:ガウシアンプロセス - 創薬・材料探索のための機械学習

となりますが、簡単に復習すると

\begin{align} f(x_{k+1} | x_{1:k}) \sim \mathcal{N}(x_{k+1}|μ_0 (x_{1:k}), Σ_0 (x_{1:k}, x_{1:k})) \end{align}

となるモデルです。この平均と分散は各点 x_i, x_jにおけるカーネル関数(共分散関数)により得られます。カーネルの選択はデータの関係性を考慮し、決定します。

Mean Function と カーネル選択

カーネルは入力空間中で点同士が強い相関を示すという特性があります。近い点はより近く、遠くの点はより遠くに写像することができます。また、カーネルは半正定値であることが必要条件です。よく使用されるカーネルをいかに示します

  • Gaussian kernel (power exponential)

\begin{align} k(x, x') = \alpha \exp ( -||x- x'||^2) \end{align}

ここで、 ||x- x'||^2 = \sum_{i=1}^{d} \alpha_{i} (x-x')^{2}です。αはカーネルのパラメータとなります。

  • Matern kernel

\begin{align} k(x, x') = \alpha_{0} \frac{2^{1-\nu}}{\Gamma(\nu)} (\sqrt{2\nu}||x-x'||)^{\nu} K_{\nu} (\sqrt{2\nu}||x-x'||) \end{align}

ここで、 K_{\nu}は修正ベッセル関数です。いずれもscikit-learnで使用できます。

良く使用される平均関数は \mu_{0} = \muですが、何らかのトレンドを有するなどの事前知識を有しているとき、パラメトリックで表した方が良い場合もあります。

\begin{align} \mu_{0} = \mu + \sum_{i=1}^{p} \beta_{i} \Psi_{i}(x) \end{align}

ここで、 \Psiは基底関数です。よく低次元の多項式が使用されます。多項式の場合では正規方程式を解けばよいので次数が小さければ簡単に求まります。

ハイパーパラメータの選択

平均関数とカーネルにはパラメータがありますが、これをハイパーパラメータ(超パラメータ)と呼びます。これを求めるための方法の1つに最尤推定(maximum likelihood estimate, MLE)があります。

\begin{align} \hat{η} = \arg \max_{η} P(f(x) | η) \end{align}

こちらは前回の記事で行いました。scipyのminimizeを使えば簡単にMLEを行えます。もう一つの方法は、事後分布最大化(maximum a posterior, MAP)推定です。事後分布を最大化させるパラメータを決定する方法です。

\begin{align} \hat{η} = \arg \max_{η} P( η | f(x) ) = \arg \max_{η} P(f(x) | η ) P(η) \end{align}

ベイズの定理を用いているがポイントです。分母の正規化定数(エビデンス)はパラメータに依存しないため省略することができます。MLEはMAPの特別な場合です。点推定か分布推定かの違いです。MCMCなどのサンプリング法で求めます。個人的にはHMCが良いと思います。

\begin{align} P(f(x)=y | f(x_{1:n})) &= \int P(f(x)=y| f(x_{1:n}), η) P(η | f(x_{1:n}) dη \\ & \approx \frac{1}{J} \sum_{j=1}^{J} P(f(x)=y| f(x_{1:n}), η=\hat{η_{j}}) \end{align}

MAP推定は完全なベイズ推論とみなせます。

獲得関数 (Acquisition Functions)

効用関数とも呼ばれています。次の点を探すための関数です。ノイズの評価、パラレルな評価、導関数を用いる方法、"exotic " extensionがあります。最もよく用いられるのが、Expected improvementです。簡単なわりに、精度も良いです。

Upper Confidence Bounds, UCB

まずは、UCBについて記載します。平均値と信頼区間の上限値を考慮した獲得関数です。

\begin{align} a_{UCB} = μ + \sqrt{\frac{\log N}{N}} σ \end{align}

データ密度の低い所を探していきます。

def ucb(X, X_sample, y_sample, gpr):
    mu, sigma = gpr.predict(X, return_std=True)
    mu_sample = gpr.predict(X_sample)

    n_sample = X_sample.shape[0]
    mu_sample_opt = np.max(mu_sample)
    z = mu_sample_opt + np.sqrt(np.log(n_sample) / n_sample) * sigma
    return z

Expected Improvement

思考実験由来のものです。ベイズ最適化で最もよく使用されています。

\begin{align} EI(x) = \mathbb{E} [\max (f(x) - f(x^{+}), 0)] \end{align}

ここで f(x^{+}) は最も良いサンプル値です。 x^{+} = \arg \max_{x_{i} \in x_{1:t}} f(x_{i})です。GPのモデルは以下の式になります。

\begin{align} EI(x) = \begin{cases} (μ\mathbf{x}) - f(\mathbf{x}^{+}) - \xi) \Phi(Z) + σ(\mathbf{x})\phi(Z) &\text{if} \ σ(\mathbf{x}) > 0 \\ 0 & \text{if}\ σ(\mathbf{x}) = 0 \end{cases} \end{align}

ここで、Zは次のようになります。

\begin{align} Z = \begin{cases} \frac{μ(\mathbf{x}) - f(\mathbf{x}^{+}) - \xi}{σ(\mathbf{x})} &\text{if}\ σ(\mathbf{x}) > 0 \\ 0 & \text{if}\ σ(\mathbf{x}) = 0 \end{cases} \end{align}

また、 \Phi, \phiはそれぞれ正規分布のCDFとPDFです。

def expected_improvement(X, X_sample, Y_sample, gpr, xi=0.01):
    mu, sigma = gpr.predict(X, return_std=True)
    mu_sample = gpr.predict(X_sample)
    sigma = sigma.reshape(-1, X_sample.shape[1])

    mu_sample_opt = np.max(mu_sample)

    with np.errstate(divide='warn'):
        imp = mu - mu_sample_opt - xi
        Z = imp / sigma
        ei = imp * norm.cdf(Z) + sigma * norm.pdf(Z)
        ei[sigma == 0.0] = 0.0

    return ei

Knowledge Gradient, KG

knowledge-gradientな獲得関数は、以前の評価点が最終的な解としてEIを最大化させようとします。この考え方は、評価がノイズフリーでリスクを嫌うときは良い方法となります。ノイズが含まれているときは提案点にノイズが含まれるため評価が困難となります。Frazier先生が提案している方法です。

\begin{align} KG (x) = \mathbb{E} [ μ_{n+1} - μ_{n} | x_{n+1} = x] \end{align}

こちらのアルゴリズムの最初はFranzier, 2012です。

f:id:udnp:20190420235207p:plain

Entropy Search

Entropy serach (ES)は情報理論に基づく方法です。

\begin{align} ES(x) = H(P(x)) - E[H(P(x|f(x)))] \end{align}

となっています。KG, ESも実装して試してみたいです。

Exotic Bayesian Optimization

上記の方法論はスタンダードなベイズ最適化です。超長方形やシンプレックス、勾配情報の不足、ノイズフリーな評価法といった前提を仮定しています。実際にはこれらの前提条件が崩れることがあり、それを考慮した方法が提案されています。

獲得関数の設計が大事になってきそうです。KGやESは下記の条件を満たせるように獲得関数を設定できます。

  • Noisy Evaluations
  • Parallel Evaluations
  • Constraints
  • Multi-Fidelity & Multi-Information Source Evaluations
  • Random Environmental Conditions and Multi-Task Bayesian Optimization
  • Derivative Observations

より適切にベイズ最適化を行うためには、こちらも詳しくサーベイしていく必要がありそうです。

ソフト

RとPythonMatlabのものが紹介されています。わりとたくさんあります。

R

  • DiceKriging
  • laGP

Python

  • Metrics Optimization Engine
  • Spearmint
  • GpyOpt
  • GPFlow
  • GPyTorch

など。

UCBとEIの例

とりあえずスタンダードの例でテストしてみます。下記関数を考えます。

\begin{align} f(x) = e^{-(x - 2)^{2}} + e^{\frac{-(x - 6)^{2}}{10}} + \frac{1}{x^{2} + 1} + noise \end{align}

f:id:udnp:20190421134812g:plain

import numpy as np
import matplotlib.pyplot as plt
from bayes_opt import BayesianOptimization
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel
from scipy.stats import norm
from scipy.optimize import minimize

noise = 0
bounds = np.array([[-2.0, 10.0]])
def f(x, noise=noise):
    return np.exp(-(x - 2) ** 2) + np.exp(-(x - 6) ** 2 / 10) + 1 / (x ** 2 + 1) + noise * np.random.randn(*x.shape)

def plot_result(X, X_init, y_init, gpr, itr):
    y_mean, y_cov = gpr.predict(X, return_std=True)
    upr = y_mean - np.sqrt(y_cov.reshape(-1, 1))
    lwr = y_mean + np.sqrt(y_cov.reshape(-1, 1))

    plt.clf()
    plt.subplot(211)
    plt.title("Bayesian Optimization # {}".format(i+1))
    plt.plot(X_init, y_init, "ro", label="Initial samples")
    plt.plot(X, y_mean, label="Surrogate model")
    plt.plot(X, f(X), label="Objective")
    plt.xlim([-2.0, 10.0])
    plt.ylim([-0.5, 1.5])
    plt.fill_between(X.ravel(), upr.ravel(), lwr.ravel(), alpha=0.5)
    plt.legend()
    plt.subplot(212)
    #plt.plot(X, ucb(X, X_init, y_init, gpr), label="UCB")
    plt.plot(X, expected_improvement(X, X_init, y_init, gpr), label="EI")
    plt.xlim([-2.0, 10.0])
    plt.legend()
    plt.savefig("./BayesOptEI{}.png".format(itr))

def next_X(acq, X_sample, Y_sample, gpr):
    dim = X_sample.shape[1]
    # 獲得関数の最大値を求めるために負にしています。
    def min_obj(X):
        return -acq(X.reshape(-1, dim), X_sample, Y_sample, gpr)

    # 局所解に陥るので、初期値を変えます。PSOやGAなどの方が良いかもしれません。
    for x0 in np.random.uniform(bounds[:, 0], bounds[:, 1], size=(25, dim)):
        res = minimize(min_obj, x0=x0, bounds=bounds, method='L-BFGS-B')
        if res.fun < 1:
            min_val = res.fun[0]
            min_x = res.x
    X_next = min_x.reshape(-1, 1)
    return X_next

def main():
    X_init = np.array([[-0.9], [1.1]])
    y_init = f(X_init)
    kernel = ConstantKernel(1.0) + RBF(length_scale=2.0, length_scale_bounds=(0, 10))
    gpr = GaussianProcessRegressor(kernel=kernel, random_state=42).fit(X_init, y_init)
    for i in range(10):
        gpr = GaussianProcessRegressor(kernel=kernel, random_state=42).fit(X_init, y_init)
        plot_result(X, X_init, y_init, gpr, i)
        X_next = next_X(expected_improvement, X_init, y_init, gpr)
        X_init = np.r_[X_init, X_next]
        y_init = np.r_[y_init, f(X_next)]

なおパッケージを使えば簡単です。

def plot_bo(f, bo):
    x = np.linspace(-2, 10, 10000)
    mean, sigma = bo._gp.predict(x.reshape(-1, 1), return_std=True)

    plt.figure(figsize=(16, 9))
    plt.plot(x, f(x))
    plt.plot(x, mean)
    plt.fill_between(x, mean + sigma, mean - sigma, alpha=0.1)
    plt.scatter(bo.space.params.flatten(), bo.space.target, c="red", s=50, zorder=10)
    plt.show()

def BayesRun():
    bo = BayesianOptimization(
        f=f,
        pbounds={"x": (-2, 10)},
        verbose=0,
        random_state=987234,
    )
    bo.maximize(n_iter=10, acq="ucb", kappa=0.1)
    plot_bo(f, bo)

結果

|   iter    |  target   |     x     |
-------------------------------------
|  6        |  1.305    |  1.651    |
|  9        |  1.401    |  1.967    |
|  12       |  1.402    |  2.019    |
|  13       |  1.402    |  2.0      |
|  14       |  1.402    |  2.0      |
=====================================

f:id:udnp:20190421141307p:plain

bootstrapサンプリング(random forestなど)を用いたベイズ最適化も行ってみようかと思いましたが、単純な例だと悲惨な結果となりそうなので、諦めました。

次回は、KG, ESの獲得関数の実装およびベイズ最適化を用いたハイパーパラメータサーチを行っていきたいです。