kivantium活動日記

プログラムを使っていろいろやります

化合物でもDeep Learningがしたい!

この記事は2017年12月15日に https://kivantium.net/deep-for-chem/ に投稿したものです。 情報が古くなっていますが、まだ参照されているようなので再掲します。

この記事はDeep Learningやっていき Advent Calendar 2017の15日目です。

Deep Learningの威力が有名になったのは画像認識コンテストでの圧勝がきっかけでしたが、今ではDeep Learningはあらゆる分野に応用され始めています。NIPS2017でもMachine Learning for Molecules and Materialsが開催されたように、物質化学における機械学習の存在感が高まりつつあります。この記事ではその一例として化学の研究にDeep Learningが使われている例を紹介していきます。

化学物質の研究に機械学習が使われる主なパターンには

  • 分子を入力するとその分子の性質を出力する
  • 分子の性質を入力するとその性質を持った分子を出力する
  • 分子を入力するとその反応を出力する

の3つがあります。それぞれについて詳しく説明します。

分子から性質を予測する

Deep Learning以前

Deep Learning以前の性質予測では、職人の温かい手作りによる特徴量が使われていました。分子の特徴ベクトルはmolecular fingerprintsと呼ばれます。molecular fingerprintsは化合物の特徴的な一部分(fragmentと呼ばれる)がその分子にあるかどうかを0/1で表したbitを並べて作られます。

(画像はFingerprints in the RDKit p.4より引用)

どのfragmentを用いるのが有効かはデータセット・問題に依存するので様々な種類のfingerprintが提案されてきました。

主なfingerprintを挙げると

などがあります。 fingerprintはRDKitなどのライブラリを使うと簡単に計算できます。(各ソフトで計算できるfingerprintのリスト

このようなfingerprintを使ってSVMやRandomforestなどでその分子がある性質を持つ/持たないを予測する研究がたくさんあります。化学の分野でDeep Learningが大きく注目されるきっかけになったのは、kaggleの薬の活性予測のコンペでHintonらのチームが優勝したことですが、論文を見ると特徴量には上のように設計されたものを使っており、ニューラルネットワークで設計されたものではなかったようです。

graph convolutionの登場

fingerprintの設計にニューラルネットワークが導入されたのが[Duvenaud+, 2015]です。この研究ではcircular fingerprint (上のECFPのこと)をもとにneural graph fingerprint (NFP)を提案しています。以下にアルゴリズムを示します。

従来のfingerprint設計でhashやmodになっていた部分が重みを調整できる演算に変更されています。これにより、予測にとって重要なfragmentの寄与は大きく、重要ではないfragmentの寄与は小さくなるような特徴量が設計できるようになりました。実際に分子の水への溶けやすさをNFPで予測したところ水への溶けやすさに影響するR-OHのような構造の重みが大きくなったことが報告されています。

NFPの他にも分子のグラフ構造に基づいたニューラルネットワークベースの特徴量設計の研究が行われています。これらはグラフ構造に注目したニューラルネットワークなので総称としてgraph convolutionと呼ばれています。一番有名なのはGoogle BrainのNeural Message Passing for Quantum Chemistryでしょう。この論文ではMessage Passing Neural Network (MPNN) というグラフ上のニューラルネットワークを提案し、分子のニューラルネットワークの先行研究の多くがMPNNで一般的に記述できることを主張した上で、MPNNが分子の性質を予測する上で高い性能を発揮すると主張しています。 MPNNは

という式で表されます。グラフ上で隣接するエッジからのメッセージ Mを足し合わせるような処理をしていることが分かります。 M, Uなどをうまく定めることで各種のgraph convolutionを表すことができます。詳細は論文を読んでください。Google ResearchによるブログPredicting Properties of Molecules with Machine Learningも役立つかもしれません。(ちなみにこの論文のラストオーサーは先述したkaggleコンペの論文のファーストオーサーです)

graph convolutionのイメージとしてよく使われる絵が[Han Altae-Tran+, 2017]にあります。一枚引用します。

Graph convolutionではない方法としては[Goh+, 2017]のような分子を画像にしてCNNで予測するようなものもあります。

ちなみに、同じ人がつい先日SMILES2Vecという文字列から化合物の性質を予測する論文も書いていました。

実装

分子に対するDeep Learningのライブラリで最も有名なのはDeepChemでしょう。DeepChemはTensorFlowでgraph convolutionを実装しています。Graph Convolutions For Tox21などのチュートリアルを読むとだいたい使い方が分かるのではないでしょうか(私も使ったことはないです)。ちなみに、なぜかPong in DeepChem with A3Cのようなチュートリアルもあり何がしたいのか謎です……

また、PFNが最近Chainer Chemistryを公開しました。NFP, GGNN, Weave, SchNetなどのgraph convolution手法が実装されているほか、QM9, Tox21などの有名どころのデータセットを使うコードも揃っており、普段Chainerを使っている人はこれを試してみるのもよいかもしれません。

性質から分子を作る

創薬などの応用においては、「タンパク質Xの動きを抑制する」などの特定の性質を持った分子を作ることが必要になります。化学物質の構造と生物学的な活性の関係のことをQSARと呼びますが、逆に活性から構造を予測する問題をinverse-QSARのように言うことがあります。

分子設計の難しさの一つは、可能な分子の数が非常にたくさんあることです。[Bohacek+, 1996]によれば、C,N,O,Sを30個以下しか持たない分子に限っても 10^{60}種類の分子が存在できるとされています。そのため全探索は不可能なので何らかの効率的な探索法を考える必要があります。

Deep Learning以外の方法

創薬は重要な研究分野なので以前から研究が行われていました。多くの手法は[Nishibata+, 1991][Pierce+, 2004]のように既に知られている部分構造を組み合わせることで分子を設計しています。最近の研究では[Kawai+, 2014]のように構造の組み合わせに遺伝的アルゴリズム構造を使ったり、[Podlewska+, 2017]のように目的関数を機械学習の予測値にしたりするなどの工夫がなされています。

Deep Learningによる方法

分子設計にDeep Leaningを持ち込んだ研究が[Gómez-Bombarelli+, 2016]です。この研究では分子の文字列表現であるSMILES記法をvariational autoencoder (VAE) を用いて実数ベクトルに変換し、ベイズ最適化で最適化したベクトルをSMILESに戻すことで分子を設計しています。この手法の問題点はVAE空間上で最適化ベクトルをSMILESに戻したときに生成される文字列が文法的に正しくないなどの理由で分子と対応しなくなる率が非常に高かったことです。

SMILES記法は、グラフ構造として表される化合物を環を切り開くなどして文字列として表現できるようにしています。OpenSMILES specificationのように文脈自由文法で規定される文法を持っており、文法に従わない文字列は分子を表しません。(なお、文法に従っていても対応する分子が化学的に存在できるかは別の問題です)。例えば下のような図で表される分子のSMILESはO1C=C[C@H]([C@H]1O2)c3c2cc(OC)c4c3OC(=O)C5=C4CCC(=O)5となります。同じ数字はそこで環を形成していることを表し、カッコは分岐を表しています。

文法的に正しくないSMILESの文字列が生成される問題を解決するために、VAEの入出力にSMILESの文字列をそのまま使うのではなくSMILESを生成する文脈自由文法の生成規則列を使うことにしたのが[Kusner+, 2017]のGrammar Variational Autoencoderです。この研究で技術的に面白いところはVAE表現から文字列を生成する際にプッシュダウンオートマトンを考えて、現在スタックの一番上にある文字から選択できない生成規則の確率を0にする工夫を導入しているところです。この工夫により生成される文字列はSMILESの文法的に正しいものに限定することができるためデコードの効率が上がるほか、潜在空間自体もよりよいものになったと主張されています。

これらのアプローチに影響されたのかは分かりませんが、分子の構造を直接設計するのではなく、分子を表すSMILESを生成する研究が盛んに行われています。

  • [Segler+, 2017] はChEMBLのSMILESを学習したLSTMで新しいSMILESを生成しています。また、薬の候補になりそうな分子を入力としたRNNのファインチューニングなども行っています。

分子から反応を予測する

分子からの反応予測には、複数の分子を入力して反応結果を出力するものと、一つの分子を入力してその分子を作るのに必要な反応を予測するものがあります。

反応結果の予測

反応予測をコンピュータで行う試みは1960年代から行われていますが、従来の手法では専門家がルールをたくさん記述することで実現しています。この分野にもDeep Learningの波が来ています。

[Schwaller+, 2017] は反応物のSMILESを入力に生成物のSMILESを出力する言語モデルを用いて反応の予測を行っています。 SMILESによる反応の記述は、反応物の文字列を入力して生成物の文字を出力する処理なので、英語を入力してフランス語を入力する処理に似ていると彼らは考えましたアメリカの特許にある反応のデータベースから入力と出力のペアを作り、seq2seqという翻訳に使われるRNNモデルを適用して反応の予測を行いました。

結果としてtop-1で80%という先行研究を上回る精度の予測ができるようになったと主張されています。

逆合成の予測

目的の化合物を合成するための反応経路を求めることをretrosynthesisといいますが、実際に化合物を生産する上では非常に重要な技術です。この研究でもDeep Learningを使った論文が出ています。

[Segler+, 2017]ではAlphaGoと似た手法でretrosynthesisを行っています。(図は論文のFigure 1) (a)は目的の化合物(図ではIbuprofen)からはじめて分子をばらしていき、全てが既知の入手可能な分子(図では赤で示されている)にまで還元できたら逆合成が完了するというコンセプトを示しています。 (b)は(a)で用いられた既知の反応を示しています。 (c)は(a)の結果得られた反応経路から実際に目的の化合物を合成する過程を示しています。 (d)がこの論文の中心となるアイデアを表しています。現在の分子をばらすのに使える既知の反応はいくつもあります。反応の各段階を一つの状態ととらえると反応はグラフ上の状態遷移と考えることができ、逆合成はグラフ上の最適な経路を探す問題と解釈できます。そこでゲームの状態を表す木から最適な手を探すのと同じような方法を用いて、最適な次の反応を選ぶことで逆合成を解くことができると考えられます。 (e)のように分子の状態を入力すると良さそうな反応を返すDNNの確率をガイドにしたモンテカルロ木探索を実行することで逆合成を行うことができそうです。

論文の実験ではモンテカルロ木探索を用いた提案手法が先行研究よりも高い性能を示したと主張されています。

私が知っている主な研究はこれくらいですが、他にも面白い研究を知っている方がいらっしゃったらコメントなどで教えて下さい。

Illustration2VecをONNX経由で使う

趣味プロジェクトでIllustration2Vecを使いたくなったのですが、これは2015年の論文なのでモデルをCaffeかChainerで使うことになっています。 github.com

残念ながらCaffeもChainerも既に開発が終了しているため、Illustration2VecのモデルをONNXという共通フォーマットに変換して今後も使えるようにしました。 利用方法だけ知りたい人は「モデルの変換」を飛ばして「使い方」を見てください。

モデルの変換

まずはオリジナルのIllustration2Vecのモデルをダウンロードします。以下を実行するとCaffeのモデルがダウンロードできます。

git clone https://github.com/rezoo/illustration2vec.git
cd illustration2vec
./get_models.sh

このスクリプトでは特徴抽出モデルのprototxtがダウンロードできなかったので、Illustration2VecのInternet ArchiveからNetwork configuration file (feature vectors) illust2vec.prototxtを追加でダウンロードしました。

必要なライブラリをインストールします。

pip install onnx coremltools onnxmltools

以下のPythonスクリプトを実行すると、タグ予測モデルのillust2vec_tag_ver200.onnxと特徴ベクトル抽出モデルのillust2vec_ver200.onnxが作成されます。

import os
import onnx
import coremltools
import onnxmltools

# CaffeモデルをONNX形式で読み込む関数
# https://github.com/onnx/onnx-docker/blob/master/onnx-ecosystem/converter_scripts/caffe_coreml_onnx.ipynb


def caffe_to_onnx(proto_file, input_caffe_path):
    output_coreml_model = 'model.mlmodel'  # 中間ファイル名
    # 中間ファイルが既に存在したら例外を送出
    if os.path.exists(output_coreml_model):
        raise FileExistsError('model.mlmodel already exists')

    # CaffeのモデルをCore MLに変換
    coreml_model = coremltools.converters.caffe.convert(
        (input_caffe_path, proto_file))
    # Core MLモデルを保存
    coreml_model.save(output_coreml_model)
    # Core MLモデルを読み込む
    coreml_model = coremltools.utils.load_spec(output_coreml_model)
    # Core MLモデルをONNXに変換
    onnx_model = onnxmltools.convert_coreml(coreml_model)
    # Core MLモデルを削除
    os.remove(output_coreml_model)
    return onnx_model


# タグ予測モデルの変換・保存
onnx_tag_model = caffe_to_onnx(
    'illust2vec_tag.prototxt', 'illust2vec_tag_ver200.caffemodel')
onnxmltools.utils.save_model(onnx_tag_model, 'illust2vec_tag_ver200.onnx')

# 特徴ベクトル抽出モデルの変換・保存
onnx_model = caffe_to_onnx('illust2vec.prototxt',
                           'illust2vec_ver200.caffemodel')
# encode1レイヤーをONNXから利用できるようにする
# https://github.com/microsoft/onnxruntime/issues/2119
intermediate_tensor_name = 'encode1'
intermediate_layer_value_info = onnx.helper.ValueInfoProto()
intermediate_layer_value_info.name = intermediate_tensor_name
onnx_model.graph.output.extend([intermediate_layer_value_info])
onnx.save(onnx_model, 'illust2vec_ver200.onnx')

使い方

上のようにしてONNX形式に変換したモデルと、それを利用するためのコードを用意しました。 github.com

ONNX形式を他のフレームワークで読み込んで実行してもいいのですが、ONNX RuntimeというMicrosoft製のパフォーマンスを重視した推論専用のライブラリがあったのでこれを使うことにしました。 UbuntuでONNX RuntimeをCPU向けにインストールするコマンドは以下の通りです。

sudo apt install libgomp1
pip install onnxruntime

例を実行する前にコードとpre-trainedモデルのダウンロードを行ってください。

git clone https://github.com/kivantium/illustration2vec.git
cd illustration2vec
./get_onnx_models.sh

タグ予測

コード

import i2v
from PIL import Image
from pprint import pprint

illust2vec = i2v.make_i2v_with_onnx(
    "illust2vec_tag_ver200.onnx", "tag_list.json")

img = Image.open("images/miku.jpg")
pprint(illust2vec.estimate_plausible_tags([img], threshold=0.5))

入力

https://github.com/kivantium/illustration2vec/blob/master/images/miku.jpg?raw=true

Hatsune Miku (初音ミク), © Crypton Future Media, INC., http://piapro.net/en_for_creators.html. This image is licensed under the Creative Commons - Attribution-NonCommercial, 3.0 Unported (CC BY-NC).

出力

[{'character': [('hatsune miku', 0.9999994039535522)],
  'copyright': [('vocaloid', 0.9999999403953552)],
  'general': [('thighhighs', 0.9956372976303101),
              ('1girl', 0.9873461723327637),
              ('twintails', 0.9812833666801453),
              ('solo', 0.9632900953292847),
              ('aqua hair', 0.9167952537536621),
              ('long hair', 0.8817101716995239),
              ('very long hair', 0.8326565027236938),
              ('detached sleeves', 0.7448851466178894),
              ('skirt', 0.6780778169631958),
              ('necktie', 0.560835063457489),
              ('aqua eyes', 0.5527758598327637)],
  'rating': [('safe', 0.9785730242729187),
             ('questionable', 0.02053523063659668),
             ('explicit', 0.0006299614906311035)]}]

Chainer版とほとんど同じ結果が出力されました。Chainerではこの処理に6秒かかっていましたが、onnx-runtimeだと2秒で実行できたのでたしかにパフォーマンスにも優れているようです(ChainerではCaffeのモデルを変換する手間が掛かっているので1枚を処理する時間で比較するのは公平ではないですが)。

特徴ベクトルの抽出

コード

import i2v
from PIL import Image

illust2vec = i2v.make_i2v_with_onnx("illust2vec_ver200.onnx")

img = Image.open("images/miku.jpg")

# extract a 4,096-dimensional feature vector
result_real = illust2vec.extract_feature([img])
print("shape: {}, dtype: {}".format(result_real.shape, result_real.dtype))
print(result_real)

# i2v also supports a 4,096-bit binary feature vector
result_binary = illust2vec.extract_binary_feature([img])
print("shape: {}, dtype: {}".format(result_binary.shape, result_binary.dtype))
print(result_binary)

先ほどと同じ入力に対する出力

shape: (1, 4096), dtype: float32
[[ 7.474596    3.6860986   0.537967   ... -0.14563629  2.7182112
   7.3140917 ]]
shape: (1, 512), dtype: uint8
[[246 215  87 107 249 190 101  32 187  18 124  90  57 233 245 243 245  54
  229  47 188 147 161 149 149 232  59 217 117 112 243  78  78  39  71  45
  235  53  49  77  49 211  93 136 235  22 150 195 131 172 141 253 220 104
  163 220 110  30  59 182 252 253  70 178 148 152 119 239 167 226 202  58
  179 198  67 117 226  13 204 246 215 163  45 150 158  21 244 214 245 251
  124 155  86 250 183  96 182  90 199  56  31 111 123 123 190  79 247  99
   89 233  61 105  58  13 215 159 198  92 121  39 170 223  79 245  83 143
  175 229 119 127 194 217 207 242  27 251 226  38 204 217 125 175 215 165
  251 197 234  94 221 188 147 247 143 247 124 230 239  34  47 195  36  39
  111 244  43 166 118  15  81 177   7  56 132  50 239 134  78 207 232 188
  194 122 169 215 124 152 187 150  14  45 245  27 198 120 146 108 120 250
  199 178  22  86 175 102   6 237 111 254 214 107 219  37 102 104 255 226
  206 172  75 109 239 189 211  48 105  62 199 238 211 254 255 228 178 189
  116  86 135 224   6 253  98  54 252 168  62  23 163 177 255  58  84 173
  156  84  95 205 140  33 176 150 210 231 221  32  43 201  73 126   4 127
  190 123 115 154 223  79 229 123 241 154  94 250   8 236  76 175 253 247
  240 191 120 174 116 229  37 117 222 214 232 175 255 176 154 207 135 183
  158 136 189  84 155  20  64  76 201  28 109  79 141 188  21 222  71 197
  228 155  94  47 137 250  91 195 201 235 249 255 176 245 112 228 207 229
  111 232 157   6 216 228  55 153 202 249 164  76  65 184 191 188 175  83
  231 174 158  45 128  61 246 191 210 189 120 110 198 126  98 227  94 127
  104 214  77 237  91 235 249  11 246 247  30 152  19 118 142 223   9 245
  196 249 255   0 113   2 115 149 196  59 157 117 252 190 120  93 213  77
  222 215  43 223 222 106 138 251  68 213 163  57  54 252 177 250 172  27
   92 115 104 231  54 240 231  74  60 247  23 242 238 176 136 188  23 165
  118  10 197 183  89 199 220  95 231  61 214  49  19  85  93  41 199  21
  254  28 205 181 118 153 170 155 187  60  90 148 189 218 187 172  95 182
  250 255 147 137 157 225 127 127  42  55 191 114  45 238 228 222  53  94
   42 181  38 254 177 232 150  99]]

Chainer版と同じbinary vectorが出力されていました。

次回はこれを使ってイラストの機械学習をします。

閉殻Hartree-Fock法によるHeH+のエネルギー計算

量子化学計算を勉強するために新しい量子化学―電子構造の理論入門〈上〉を読んでいます。
サンプルコードがFortranだったので勉強がてらC言語に移植しました。

新しい量子化学―電子構造の理論入門〈上〉

新しい量子化学―電子構造の理論入門〈上〉

理論

本で200ページくらいかけて説明されている内容をブログに書くのは大変なので省略します。

ちなみに、以下のツイートの「100ページかけて説明してる別の本」はこの本を指します。
難しい内容を易しく簡潔に説明しようとするのは大変で、難しい内容は難しいまま長い時間をかけて説明する必要があるのだなぁという学びがありました。


実装

本のFortran実装を素直にC言語で書き直しました。
本のコードを打ち込むのが面倒だったので、検索して見つけたページにあったコードと見比べながら書きました。

最終的に得られたエネルギーはFortran実装の結果と小数第5位まで一致しました。
この差がどこで生じているのか謎ですが、実装は正しいんじゃないでしょうか……

#include <math.h>
#include <stdio.h>

// プロトタイプ宣言
void HFCALC(int, int, double, double, double, double, double);
void INTGRL(int, int, double, double, double, double, double);
void COLECT(int, int, double, double, double, double, double);
void SCF(int, int, double, double, double, double, double);
double calcS(double, double, double);
double calcT(double, double, double);
double F0(double);
double calcV(double, double, double, double, double);
double TWOE(double, double, double, double, double, double, double);
void FORMG(void);
void DIAG(double [][2], double [][2], double [][2]);
void MULT(double [][2], double [][2], double [][2], int, int);
void MATOUT(double [][2], int, int, int, int, char*);

double S12, T11, T12, T22, V11A, V12A, V22A, V11B, V12B, V22B,
       V1111, V2111, V2121, V2211, V2221, V2222;

// Slater型関数に対する最小2乗近似の短縮係数と軌道指数
double COEF[3][3] = {{1.0, 0.678914, 0.444635}, 
                     {0.0, 0.430129, 0.535328},
                     {0.0, 0.0, 0.154329}};
double EXPON[3][3] = {{0.27095, 0.151623, 0.109818},
                      {0.0, 0.851819, 0.405771},
                      {0.0, 0.0, 2.22766}};
double D1[3], A1[3], D2[3], A2[3];
double PI = 3.1415926535898;

double S[2][2], X[2][2], XT[2][2], H[2][2], F[2][2], G[2][2], C[2][2],
       FPRIME[2][2], CPRIME[2][2], P[2][2], OLDP[2][2], TT[2][2][2][2], E[2][2];

int main(void) {
  int IOP = 2;
  // 基底関数STO-NG
  int N = 3;
  // 基底状態の平衡結合長(a.u.)
  double R = 1.4632;
  // Heに対するSTO-3G軌道指数値
  // He原子に対する最良の軌道指数は1.6875
  // HeH+では電子雲が縮んでいることを反映して1.24倍
  double ZETA1 = 2.0925;
  // Hに対するSTO-3G軌道指数値
  // H原子に対する最良の軌道指数は1.0
  // HeH+では電子雲が縮んでいることを反映して1.24倍
  double ZETA2 = 1.24;
  double ZA = 2.0;
  double ZB = 1.0;
  HFCALC(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
}

/*
 * 1s基底関数に対応するSTO-NG最小基底を用いて2電子系2原子分子に対する
 * Hartree-Fock計算を行う関数
 *
 * IOP=0 何も表示しない
 * IOP=1 収束の結果のみ表示
 * IOP=2 全繰り返しの結果を表示
 * N STO-NG を使う (N=1, 2, 3)
 * R 結合長 (a.u.)
 * ZETA1 原子Aの基底関数のSlater軌道指数値
 * ZETA2 原子Bの基底関数のSlater軌道指数値
 * ZA 原子Aの原子番号
 * ZB 原子Bの原子番号
 */
void HFCALC(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
            double ZB) {
  if (IOP != 0) {
    printf("   STO-%dG for atomic numbers %5.2lf and %5.2lf\n\n\n", N, ZA, ZB);
  }
  INTGRL(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
  COLECT(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
  SCF(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
}

/* SCF計算に必要な積分を計算する関数*/
void INTGRL(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
            double ZB) {
  double R2 = R * R;
  // 軌道指数のスケーリングと短縮係数の規格化
  for (int i = 0; i < N; ++i) {
    A1[i] = EXPON[i][N-1] * pow(ZETA1, 2.0);
    D1[i] = COEF[i][N-1] * pow(2.0 * A1[i] / PI, 0.75);
    A2[i] = EXPON[i][N-1] * pow(ZETA2, 2.0);
    D2[i] = COEF[i][N-1] * pow(2.0 * A2[i] / PI, 0.75);
  }
  S12 = 0.0;
  T11 = 0.0;
  T12 = 0.0;
  T22 = 0.0;
  V11A = 0.0;
  V12A = 0.0;
  V22A = 0.0;
  V11B = 0.0;
  V12B = 0.0;
  V22B = 0.0;
  V1111 = 0.0;
  V2111 = 0.0;
  V2121 = 0.0;
  V2211 = 0.0;
  V2221 = 0.0;
  V2222 = 0.0;

  // 1電子積分の計算
  // V12A = 原子核Aへの各引力の非対角項
  // RAP2 = 中心Aと中心Pの距離の2乗
  for (int i = 0; i < N; ++i) {
    for (int j = 0; j < N; ++j) {
      double RAP = A2[j] * R / (A1[i] + A2[j]);
      double RAP2 = pow(RAP, 2.0);
      double RBP2 = pow(R - RAP, 2.0);
      S12 += calcS(A1[i], A2[j], R2) * D1[i] * D2[j];
      T11 += calcT(A1[i], A1[j], 0.0) * D1[i] * D1[j];
      T12 += calcT(A1[i], A2[j], R2) * D1[i] * D2[j];
      T22 += calcT(A2[i], A2[j], 0.0) * D2[i] * D2[j];
      V11A += calcV(A1[i], A1[j], 0.0, 0.0, ZA) * D1[i] * D1[j];
      V12A += calcV(A1[i], A2[j], R2, RAP2, ZA) * D1[i] * D2[j];
      V22A += calcV(A2[i], A2[j], 0.0, R2, ZA) * D2[i] * D2[j];
      V11B += calcV(A1[i], A1[j], 0.0, R2, ZB) * D1[i] * D1[j];
      V12B += calcV(A1[i], A2[j], R2, RBP2, ZB) * D1[i] * D2[j];
      V22B += calcV(A2[i], A2[j], 0.0, 0.0, ZB) * D2[i] * D2[j];
    }
  }
  // 2電子積分の計算
  for (int i = 0; i < N; ++i) {
    for (int j = 0; j < N; ++j) {
      for (int k = 0; k < N; ++k) {
        for (int l = 0; l < N; ++l) {
          double RAP = A2[i] * R / (A2[i] + A1[j]);
          //double RBP = R - RAP;
          double RAQ = A2[k] * R / (A2[k] + A1[l]);
          double RBQ = R - RAQ;
          double RPQ = RAP - RAQ;
          double RAP2 = RAP * RAP;
          //double RBP2 = RBP * RBP;
          //double RAQ2 = RAQ * RAQ;
          double RBQ2 = RBQ * RBQ;
          double RPQ2 = RPQ * RPQ;
          V1111 += TWOE(A1[i], A1[j], A1[k], A1[l], 0.0, 0.0, 0.0) 
                   * D1[i] * D1[j] * D1[k] * D1[l];
          V2111 += TWOE(A2[i], A1[j], A1[k], A1[l], R2, 0.0, RAP2)
                   * D2[i] * D1[j] * D1[k] * D1[l];
          V2121 += TWOE(A2[i], A1[j], A2[k], A1[l], R2, R2, RPQ2)
                   * D2[i] * D1[j] * D2[k] * D1[l];
          V2211 += TWOE(A2[i], A2[j], A1[k], A1[l], 0.0, 0.0, R2)
                   * D2[i] * D2[j] * D1[k] * D1[l];
          V2221 += TWOE(A2[i], A2[j], A2[k], A1[l], 0.0, R2, RBQ2)
                   * D2[i] * D2[j] * D2[k] * D1[l];
          V2222 += TWOE(A2[i], A2[j], A2[k], A2[l], 0.0, 0.0, 0.0) * D2[i] *
                   D2[j] * D2[k] * D2[l];
        }
      }
    }
  }
  if (IOP != 0) {
    printf("   R          ZETA1      ZETA2      S12       T11\n\n");
    printf("%11.6lf%11.6lf%11.6lf%11.6lf%11.6lf\n\n\n", R, ZETA1, ZETA2, S12, T11);
    printf("   T12        T22        V11A       V12A        V22A\n\n");
    printf("%11.6lf%11.6lf%11.6lf%11.6lf%11.6lf\n\n\n", T12, T22, V11A, V12A, V22A);
    printf("   V11B       V12B       V22B       V1111        V2111\n\n");
    printf("%11.6lf%11.6lf%11.6lf%11.6lf%11.6lf\n\n\n", V11B, V12B, V22B, V1111,
           V2111);
    printf("   V2121      V2211       V2221     V2222\n\n");
    printf("%11.6lf%11.6lf%11.6lf%11.6lf\n", V2121, V2211, V2221, V2222);
  }
}

// 重なり積分
double calcS(double A, double B, double RAB2) {
  return pow(PI / (A + B), 1.5) * exp(-A * B * RAB2 / (A + B));
}

// 運動エネルギー積分
double calcT(double A, double B, double RAB2) {
  return A * B / (A + B) * (3.0 - 2.0 * A * B * RAB2 / (A + B)) *
         pow(PI / (A + B), 1.5) * exp(-A*B*RAB2/(A+B));
}

// F0関数
double F0(double ARG) {
  double ret;
  if (ARG < 1.0e-6) {
    ret = 1.0 - ARG / 3.0;
  } else {
    ret = sqrt(PI / ARG) * erf(sqrt(ARG)) / 2.0;
  }
  return ret;
}

// 核引力積分
double calcV(double A, double B, double RAB2, double RCP2, double ZC) {
  double v = 2.0 * PI / (A + B) * F0((A + B) * RCP2) * exp(-A * B * RAB2 / (A + B));
  return -v * ZC;
}

// 2電子積分
double TWOE(double A, double B, double C, double D, double RAB2, double RCD2,
            double RPQ2) {
  return 2.0 * pow(PI, 2.5) / ((A + B) * (C + D) * sqrt(A + B + C + D))
         * F0((A + B) * (C + D) * RPQ2 / (A + B + C + D))
         * exp(-A * B * RAB2 / (A + B) - C * D * RCD2 / (C + D));
}

// 積分の結果から行列を求める
void COLECT(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
            double ZB) {
  // 核ハミルトニアン cf. 式(3.153)
  H[0][0] = T11 + V11A + V11B;
  H[0][1] = T12 + V12A + V12B;
  H[1][0] = H[0][1];
  H[1][1] = T22 + V22A + V22B;

  // 重なり行列 cf. 式(3.136)
  S[0][0] = 1.0;
  S[0][1] = S12;
  S[1][0] = S[0][1];
  S[1][1] = 1.0;

  // 正準直交化 cf. 式(3.169), p.191
  X[0][0] = 1.0 / sqrt(2.0 * (1.0 + S12));
  X[1][0] = X[0][0];
  X[0][1] = 1.0 / sqrt(2.0 * (1.0 - S12));
  X[1][1] = -X[0][1];

  // 変換行列の転置
  XT[0][0] = X[0][0];
  XT[0][1] = X[1][0];
  XT[1][0] = X[0][1];
  XT[1][1] = X[1][1];

  // 2電子積分の行列
  TT[0][0][0][0] = V1111;
  TT[1][0][0][0] = V2111;
  TT[0][1][0][0] = V2111;
  TT[0][0][1][0] = V2111;
  TT[0][0][0][1] = V2111;
  TT[1][0][1][0] = V2121;
  TT[0][1][1][0] = V2121;
  TT[1][0][0][1] = V2121;
  TT[0][1][0][1] = V2121;
  TT[1][1][0][0] = V2211;
  TT[0][0][1][1] = V2211;
  TT[1][1][1][0] = V2221;
  TT[1][1][0][1] = V2221;
  TT[1][0][1][1] = V2221;
  TT[0][1][1][1] = V2221;
  TT[1][1][1][1] = V2222;

  if (IOP != 0) {
    MATOUT(S, 2, 2, 2, 2, "S   ");
    MATOUT(X, 2, 2, 2, 2, "X   ");
    MATOUT(H, 2, 2, 2, 2, "H   ");
    printf("\n");
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        for (int k = 0; k < 2; ++k) {
          for (int l = 0; l < 2; ++l) {
            printf("   (%2d%2d%2d%2d )%10.6lf\n", i+1, j+1, k+1, l+1, TT[i][j][k][l]);
          }
        }
      }
    }
  }
}

// SCFを実行する
void SCF(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
         double ZB) {
  // 密度行列の収束規準
  double CRIT = 1.0e-4;
  // 繰り返しの最大回数
  int MAXIT = 25;
  // 核ハミルトニアンを最初の推定に用いる(つまりPを零行列で初期化)
  for (int i = 0; i < 2; ++i) {
    for (int j = 0; j < 2; ++j) {
      P[i][j] = 0.0;
    }
  }
  if (IOP == 2) {
    MATOUT(P, 2, 2, 2, 2, "P    ");
  }
  int ITER = 0;
  double EN;
  for (;;) {
    ITER++;
    if (IOP == 2) {
      printf("\n    START OF ITERATION NUMBER = %2d\n", ITER);
    }
    // Fock行列の2電子部分を計算
    FORMG();

    if (IOP == 2) {
      MATOUT(G, 2, 2, 2, 2, "G    ");
    }

    // 核ハミルトニアンを足してFock行列を得る
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        F[i][j] = H[i][j] + G[i][j];
      }
    }

    // 電子エネルギーを計算する
    EN = 0.0;
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        EN += 0.5 * P[i][j] * (H[i][j] + F[i][j]);
      }
    }

    if (IOP == 2) {
      MATOUT(F, 2, 2, 2, 2, "F   ");
      printf("\n\n\n    ELECTRONIC ENERGY = %20.12lf\n", EN);
    }

    // Gを用いてFock行列を変換
    MULT(F, X, G, 2, 2);
    MULT(XT, G, FPRIME, 2, 2);
    // Fock行列の対角化
    DIAG(FPRIME, CPRIME, E);
    // 固有ベクトルから行列Cを得る
    MULT(X, CPRIME, C, 2, 2);
    // 新しい密度行列を計算する
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        OLDP[i][j] = P[i][j];
        P[i][j] = 0.0;
        for (int k = 0; k < 1; ++k) { // ここはk<1で本当にいいのか不明
          P[i][j] += 2.0 * C[i][k] * C[j][k];
        }
      }
    }
    if (IOP == 2) {
      MATOUT(FPRIME, 2, 2, 2, 2, "F'   ");
      MATOUT(CPRIME, 2, 2, 2, 2, "C'   ");
      MATOUT(E, 2, 2, 2, 2, "E   ");
      MATOUT(C, 2, 2, 2, 2, "C   ");
      MATOUT(P, 2, 2, 2, 2, "P   ");
    }
    // deltaを計算する
    double DELTA = 0.0;
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        DELTA += pow(P[i][j] - OLDP[i][j], 2.0);
      }
    }
    DELTA = sqrt(DELTA / 4.0);
    if (IOP != 0) {
      printf("\n    DELTA(CONVERGENCE OF DENSITY MATRIX) = %10.6lf\n", DELTA);
    }

    // 収束判定
    if (DELTA < CRIT) break;
    if (ITER < MAXIT) continue;
    // 収束しなかった場合
    printf("    NO CONVERGENCE IN SCF");
    return;
  }

  double ENT = EN + ZA * ZB / R;
  if (IOP != 0) {
    printf("\n\n    CALCULATION CONVERGED"
           "\n\n    ELECTORONIC ENERGY = %20.12lf"
           "\n\n    TOTAL ENERGY =       %20.12f", EN, ENT);
  }
  if (IOP == 1) {
    MATOUT(G, 2, 2, 2, 2, "G   ");
    MATOUT(F, 2, 2, 2, 2, "F   ");
    MATOUT(E, 2, 2, 2, 2, "E   ");
    MATOUT(C, 2, 2, 2, 2, "C   ");
    MATOUT(P, 2, 2, 2, 2, "P   ");
  }
  MULT(P, S, OLDP, 2, 2);
  if (IOP != 0) {
    MATOUT(OLDP, 2, 2, 2, 2, "PS  ");
  }
}

void FORMG(void) {
  for (int i = 0; i < 2; ++i) {
    for (int j = 0; j < 2; ++j) {
      G[i][j] = 0.0;
      for (int k = 0; k < 2; ++k) {
        for (int l = 0; l < 2; ++l) {
          G[i][j] += P[k][l] * (TT[i][j][k][l] - 0.5 * TT[i][l][k][j]);
        }
      }
    }
  }
}

// Fを対角化し、Cに固有ベクトル、Eに固有値を入れる
void DIAG(double F[][2], double C[][2], double E[][2]) {
  double theta;
  if (fabs(F[0][0] - F[1][1]) > 1.0e-20) {
    theta = 0.5 * atan(2.0 * F[0][1] / (F[0][0] - F[1][1]));
  } else {
    theta = PI / 4.0;
  }
  C[0][0] = cos(theta);
  C[1][0] = sin(theta);
  C[0][1] = sin(theta);
  C[1][1] = -cos(theta);
  E[0][0] = F[0][0] * pow(cos(theta), 2.0) + F[1][1] * pow(sin(theta), 2.0)
            + F[0][1] * sin(2.0 * theta);
  E[1][1] = F[1][1] * pow(cos(theta), 2.0) + F[0][0] * pow(sin(theta), 2.0)
            - F[0][1] * sin(2.0 * theta);
  E[1][0] = 0.0;
  E[0][1] = 0.0;
  // 固有値と固有ベクトルを整列する
  if (E[1][1] > E[0][0]) return;
  double temp = E[1][1];
  E[1][1] = E[0][0];
  E[0][0] = temp;
  temp = C[0][1];
  C[0][1] = C[0][0];
  C[0][0] = temp;
  temp = C[1][1];
  C[1][1] = C[1][0];
  C[1][0] = temp;
}

// AとBの積をCに入れる
void MULT(double A[][2], double B[][2], double C[][2], int IM, int M) {
  for (int i = 0; i < M; ++i) {
    for (int j = 0; j < M; ++j) {
      C[i][j] = 0.0;
      for (int k = 0; k < M; ++k) {
        C[i][j] += A[i][k] * B[k][j];
      }
    }
  }
}

void MATOUT(double A[][2], int IM, int IN, int M, int N, char *LABEL) {
  int IHIGH = 0;
  for (;;) {
    int LOW = IHIGH + 1;
    IHIGH = IHIGH + 5;
    IHIGH = (IHIGH < N) ? IHIGH : N;
    printf("\n\n\n    THE %s ARRAY\n", LABEL);
    printf("               ");
    for (int i = LOW; i <= IHIGH; ++i) {
      printf("          %3d      ", i);
    }
    printf("\n");
    for (int i = 1; i <= M; ++i) {
      printf("%10d     ", i);
      for (int j = LOW; j <= IHIGH; ++j) {
        printf(" %18.10lf", A[i - 1][j - 1]);
      }
      printf("\n");
    }
    if (N <= IHIGH) break;
  }
}

最終的な出力は

    CALCULATION CONVERGED

    ELECTORONIC ENERGY =      -4.227529304014

    TOTAL ENERGY =            -2.860662163500


    THE PS   ARRAY
                           1                  2      
         1            1.5296370381       1.1199277981
         2            0.6424383868       0.4703629619

でした。参照したFortranコードの出力は

    CALCULATION CONVERGED

    ELECTRONIC ENERGY =  -0.422752913203D+01

    TOTAL ENERGY =       -0.286066199152D+01



    THE PS   ARRAY
                           1                  2
         1        0.1529637121D+01   0.1119927796D+01
         2        0.6424383099D+00   0.4703628793D+00

でした。

FortranからCに移植するときのtips

詰まったことを挙げておきます。

配列の添字

Fortranでは添字が1から始まるので配列にアクセスするときに添字をうまく調整する必要があります。

多次元配列の並び順

このページに詳しく書いてありますが、C言語の多次元配列はメモリ上でa[0][0], a[0][1], a[1][0], a[1][1], a[2][0], a[2][1]のような順番で並びますが、Fortranではa(1,1), a(2,1), a(3,1), a(1,2), a(2,2), a(3,2)のような順番に並びます。多次元配列の初期化を書き直すときに気をつける必要があります。

整合配列

Fortranではサブルーチンに配列を渡した配列のサイズをサブルーチン側で定義することができます。(整合配列という)

      dimension a(6)
      a(1) = 1
      a(2) = 2
      a(3) = 3
      a(4) = 4
      a(5) = 5
      a(6) = 6
      call sub(a,3,2)
      end

      subroutine sub(a,n,m)
      dimension a(n,m) <---- a(n,*) でも同じ結果を得る
      do i=1,n
         write(6,*) (a(i,j),j=1,m)
      end do
      return
      end

コードはFortran プログラミングの基礎知識から引用しました。

C99で導入された可変長配列の機能を使うと同じことができるそうです。cf. 3.3 整合配列、可変長配列
今回は配列のサイズが決め打ちできたのでvoid MULT(double A[][2], double B[][2], double C[][2], int IM, int M) { のように関数を定義しました。

C言語のabs

特にFortranとは関係がないのですが、C言語のabsはint型に対して使う関数で、double型に対してはfabsという別の関数が用意されています。(これはC言語で関数オーバーロードができないことが原因です。C11で導入されたGenericという機能を使うとそれっぽいことが実現できるらしいです。cf. Generic - C言語入門

TwitterにMP4動画をアップロードするにはyuv420pを使う必要がある(らしい)

以前の記事でMP4をアニメーションGIFに変換する方法を紹介しました。
kivantium.hateblo.jp

その後Twitterの機能変更でMP4の動画がそのままアップロードできるようになったのですが、アップロードに失敗することが非常に多かったです。

Twitterの公式ドキュメント(How to share and watch videos on Twitter)には

モバイルアプリではMP4とMOVの動画形式をサポートしています。

ブラウザではMP4(H264形式、AACオーディオ)をサポートしています。アップロードできる動画のサイズは最大512MBです。ただし長さは2分20秒間以下にしてください。

最小解像度: 32 x 32
最大解像度: 1920 x 1200(および1200 x 1900)
縦横比: 1:2.39~2.39:1の範囲(両方の値を含む)
最大フレームレート: 40fps
最大ビットレート: 25Mbps

としか書いていないのですが、実際にはフォーマットをyuv420pにしないとアップロードできないらしいです。[要出典]
(そのように書いてあるTwitter公式のドキュメントは見当たらないのですが、多くのサイトでそう指摘されていました。)

今まで失敗していたのはffmpegのデフォルトだとyuv444になってしまうからのようです。

H.264で作ったのにアップロードに失敗した動画をyuv420pにするのは簡単で、

ffmpeg -i input.mp4 -pix_fmt yuv420p output.mp4

とするだけです。

コーデックの変換も含めて行うには

ffmpeg -i input.webm -vcodec libx264 -pix_fmt yuv420p -strict -2 -acodec aac output.mp4

のようにすればいいようです。

ボトムアップ型自動微分の実験

Deep Learningの関係で自動微分が脚光を浴びつつあるような気がしますが、自動微分を解説したページは少なくまだまだマイナーな分野だと思います。先日ようやく「アルゴリズムの自動微分と応用」を一通り眺めたのでいろいろ実験して遊んでいます。今日は自動微分のうち、ボトムアップ型自動微分をオペレータオーバーロードを用いて実現する方法について書きます。

自動微分とは

自動微分は、アルゴリズムによって定義された関数からその関数の偏導関数値を計算するアルゴリズムを導出するための技術です。一般的にはy=x^2のような関数が先にあって、その関数を計算するアルゴリズムやプログラムがあるというように考えますが、自動微分ではどちらかというとアルゴリズムが先にあってアルゴリズムが表現する関数が生まれるというような考え方をします。

プログラムで微分を扱う上でよく知られている技術には数値微分と数式微分があります。

数値微分

数値微分は、小さな値hを用いて{\displaystyle \frac{f(a+h)-f(a)}{h}}{\displaystyle \frac{f(a+h)-f(a-h)}{2h}}を計算することで微分の近似値を得る技術です。数値微分はあくまで微分の近似値を計算する技術ですが、自動微分ではアルゴリズムで定義された関数から解析的に得た偏導関数の値を計算します。また、数値微分では摂動幅hを設定する必要がありますが、hが大きすぎると近似精度が悪くなり、小さすぎると丸め誤差の影響で精度が悪くなるという問題があります。

この問題を理解するために簡単なプログラムを書きました。{\displaystyle x=\frac{\pi}{4}}における\sin(x)微分{\displaystyle \frac{f(a+h)-f(a)}{h}}を計算することで求めます。解析的には{\displaystyle \frac{\sqrt{2}}{2}}が解なので、この2倍を出力するようにしてよく知っている\sqrt{2}の値にどれくらい近づくかを見てみます。

#include <iostream>
#include <iomanip>
#include <cmath>

using namespace std;

int main(void){
    // 小数点以下15桁表示する
    cout << fixed << setprecision(15);
    
    // 刻み幅を小さくして数値微分
    double h = 0.1;
    for(int i=0; i<15; ++i) {
        cout << 2 * (sin(M_PI/4 + h) - sin(M_PI/4))/h << " (h=" << h << ")" << endl;
        h /= 10;
    }
}

結果は

1.341205945807979 (h=0.100000000000000)
1.407118983378419 (h=0.010000000000000)
1.413506219948735 (h=0.001000000000000)
1.414142849338607 (h=0.000100000000000)
1.414206491290315 (h=0.000010000000000)
1.414212855488372 (h=0.000001000000000)
1.414213490757987 (h=0.000000100000000)
1.414213568473599 (h=0.000000010000000)
1.414213635086980 (h=0.000000001000000)
1.414215411443819 (h=0.000000000100000)
1.414224293228016 (h=0.000000000010000)
1.414202088767524 (h=0.000000000001000)
1.414424133372449 (h=0.000000000000100)
1.421085471520200 (h=0.000000000000010)
1.554312234475219 (h=0.000000000000001)

のようになりました。大きすぎても小さすぎても良くなくて、h=0.000000010000000のときの結果が最も良いようです。

一般に、{\displaystyle \frac{f(a+h)-f(a)}{h}}の計算誤差をなるべく小さくするためには関数の計算値に含まれる計算誤差を\deltaとして{\displaystyle h \simeq 2\sqrt{\frac{\delta}{|f''(a)|}}}とするのがよいことが知られています。しかし、f'を求めるのにf''の値が必要になってしまうのでこの式を使うのは現実的ではありません。

また、入力変数が多い関数に対しては入力変数の数だけ関数を計算する必要があり計算量が増えてしまいますが、自動微分を使うことで計算量を抑えることができます。

数式微分

MathematicaMaximaのようなソフトでは数式を処理して偏導関数を求める数式微分を行いますが、自動微分では繰り返しや分岐を含むようなアルゴリズムに対しても偏導関数を求めるアルゴリズム偏導関数そのものではない)を求めます。

自動微分の種類

自動微分にはボトムアップ型(フォーワードモードとも)とトップダウン型(リバースモードとも)の2種類があります。どちらも微分の連鎖律を用いて関数を単純な演算に分解することで偏導関数値を求めますが、分解の方法が異なります。

ボトムアップ型自動微分

ボトムアップ型自動微分では一つの入力変数に対して、全ての中間変数の偏導関数値を計算していく手法です。

例としてWikipediaの自動微分のページにあるf(x_1, x_2)=x_1 x_2+\sin{x_1}の例を考えます。この関数の計算をステップごとに中間変数に保存するとすれば、

  •  w_1 \leftarrow x_1
  •  w_2 \leftarrow x_2
  •  w_3 \leftarrow w_1 \cdot w_2
  •  w_4 \leftarrow \sin{w_1}
  •  w_5 \leftarrow w_3 + w_4

のように5つの中間変数を用いて表すことができます。

ここで入力変数x_1について各中間変数の微分を計算すると、

  •  {\displaystyle \frac{\partial w_1}{\partial x_1} = 1}
  •  {\displaystyle \frac{\partial w_2}{\partial x_1} = 0}
  •  {\displaystyle \frac{\partial w_3}{\partial x_1}= \frac{\partial (w_1 \cdot w_2 )}{\partial x_1} = \frac{\partial w_1}{\partial x_1}\cdot w_2 + w_1 \cdot \frac{\partial w_2}{\partial x_1}}
  •  {\displaystyle \frac{\partial w_4}{\partial x_1} = \frac{\partial \sin{w_1}}{\partial x_1}  = \frac{\partial w_1}{\partial x_1} \cdot \cos{w_1}}
  •  {\displaystyle \frac{\partial w_5}{\partial x_1} = \frac{\partial (w_3+w_4)}{\partial x_1} = \frac{\partial w_3}{\partial x_1} + \frac{\partial w_4}{\partial x_1}}

のように分解することができます。これは関数値を求める際に必要な各中間変数の計算と同時に行うことができます。

計算グラフを使って表すと以下のようになります。(Wikipediaから引用)
f:id:kivantium:20160325000216p:plain:w600

トップダウン型自動微分

トップダウン型自動微分では一つの出力変数について、全ての中間変数に対する偏導関数値を計算していく手法です。

先ほどと同じ例を考えます。出力変数w_5について各中間変数に対する微分を計算すると、

  •  {\displaystyle \frac{\partial w_5}{\partial w_5} = 1}
  •  {\displaystyle \frac{\partial w_5}{\partial w_4} = \frac{\partial w_5}{\partial w_5}\cdot \frac{\partial w_5}{\partial w_4} = \frac{\partial w_5}{\partial w_5}\cdot\frac{\partial (w_3+w_4)}{\partial w_4}=\frac{\partial w_5}{\partial w_5}}
  •  {\displaystyle \frac{\partial w_5}{\partial w_3} = \frac{\partial w_5}{\partial w_5}\cdot \frac{\partial (w_3+w_4)}{\partial w_3} = \frac{\partial w_5}{\partial w_5}}
  •  {\displaystyle \frac{\partial w_5}{\partial w_2} = \frac{\partial w_5}{\partial w_3} \cdot \frac{\partial w_3}{\partial w_2} = \frac{\partial w_5}{\partial w_3} \cdot \frac{\partial (w_1 \cdot w_2)}{\partial w_2} = \frac{\partial w_5}{\partial w_3}\cdot w_1}
  •  {\displaystyle \frac{\partial w_5}{\partial w_1} = \frac{\partial w_5}{\partial w_3} \cdot \frac{\partial w_3}{\partial w_1} + \frac{\partial w_5}{\partial w_4} \cdot \frac{\partial w_4}{\partial w_1}= \frac{\partial w_5}{\partial w_3} \cdot \frac{\partial (w_1 \cdot w_2)}{\partial w_1}+ \frac{\partial w_5}{\partial w_4} \cdot \frac{\partial \sin{w_1}}{\partial w_1} = \frac{\partial w_5}{\partial w_3} \cdot w_2 + \frac{\partial w_5}{\partial w_4} \cdot \cos{w_1}}

のように分解することができます。

w_1, \ w_2についての分解が分かりにくいですが、w_5の変化は中間変数w_3, \ w_4の変化を経由してしか起こらないことに注目すると納得できるかもしれません。計算グラフは以下の通りです。
f:id:kivantium:20160325002547p:plain:w600

トップダウン型自動微分では一度計算した後で各中間変数がどのように計算されたのかの履歴を逆向きに辿って偏導関数値を計算する必要があるため、計算グラフを保存するためのメモリが必要になります。

勾配を求める際には、入力変数の方が多いときはトップダウン型、出力変数の方が多い時はボトムアップ型を使うのが計算量の点で有利になります。ニューラルネットワークバックプロパゲーショントップダウン型の自動微分の一種です。

自動微分の実現方法

自動微分の実現方法には、ソースコードを解析して偏導関数値を求めるプログラムを作る方法とオペレータオーバーロードによる方法があります。ソースコードを解析する方法はコンパイラを作るのと同じような労力が必要になるので難しいですが、オペレータオーバーロードによる方法は比較的容易なので今回はこちらを採用します。

オペレータオーバーロードによるボトムアップ型自動微分

ボトムアップ型自動微分では関数値を計算するのと同時に偏導関数値を計算することができます。そのためdoubleを自動微分用に拡張したBUdoubleというクラスをつくって計算を行うと同時に偏導関数値を保存するようにします。それを実現したのが以下のコードです。

#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;
 
// ボトムアップ型微分積分用doubleクラス
class BU_double {
    // 変数値
    double val;
    // 偏導関数値
    double d_val;
public:
    // コンストラクタ
    BU_double(double v=0, double dv=0) {
        val = v;
        d_val = dv;
    }
    // 微分する入力変数として選択する関数
    void select(void) {
        d_val = 1.0;
    }
    // 変数値を返す
    double get_value(void) {
        return val;
    }
    // 偏導関数値を返す
    double get_d_value(void) {
        return d_val;
    }

    // 各種演算子の定義
    friend BU_double operator + (BU_double x, BU_double y) {
        return BU_double(x.val+y.val, x.d_val+y.d_val);
    }
    friend BU_double operator - (BU_double x, BU_double y) {
        return BU_double(x.val-y.val, x.d_val-y.d_val);
    }
    friend BU_double operator * (BU_double x, BU_double y) {
        return BU_double(x.val*y.val, x.d_val*y.val+x.val*y.d_val);
    }
    friend BU_double operator / (BU_double x, BU_double y) {
        double w = x.val/y.val;
        return BU_double(w, (x.d_val-w*y.d_val)/y.val);
    }
    friend BU_double operator + (BU_double x) {
        return BU_double(x.val, x.d_val);
    }
    friend BU_double operator - (BU_double x) {
        return BU_double(-x.val, -x.d_val);
    }
    friend bool operator < (BU_double x, BU_double y) {
        return x.val < y.val;
    }
    friend bool operator <= (BU_double x, BU_double y) {
        return x.val <= y.val;
    }
    friend bool operator > (BU_double x, BU_double y) {
        return x.val > y.val;
    }
    friend bool operator >= (BU_double x, BU_double y) {
        return x.val >= y.val;
    }
    
    // 基本関数の定義
    friend BU_double sqrt(BU_double x) {
        return BU_double(sqrt(x.val), 0.5*x.d_val/sqrt(x.val));
    }
    friend BU_double exp(BU_double x) {
        return BU_double(exp(x.val), x.d_val*exp(x.val));
    }
    friend BU_double log(BU_double x) {
        return BU_double(log(x.val), x.d_val/x.val);
    }
    friend BU_double sin(BU_double x) {
        return BU_double(sin(x.val), cos(x.val));
    }
    friend BU_double cos(BU_double x) {
        return BU_double(cos(x.val), -sin(x.val));
    }

    // coutに出力するフォーマットの定義
    friend ostream& operator<<(ostream &s, BU_double x) {
        return s << "BU_double("<< x.val << ", " << x.d_val << ") ";
    }

};

int main(void){
    cout << setprecision(15);
    BU_double x, y;

    // sin(x)の微分 
    x = M_PI / 4;
    x.select();
    y = sin(x);
    cout << 2 * y.get_d_value() << endl; // 2016-07-13に修正

    
    // y = 10*x^2 のfor文を使ったアホな定義
    x = 5;
    y = 0;
    x.select();
    for(int i=0; i<10; ++i) y = y + x*x;
    // x = 5.0のときのyとdy/dxを出力
    cout << y << endl;

    // ニュートン法の実行
    x = 2.0; // 適当な初期値
    for(int i=0; i<10000; ++i) {
        x.select();
        // y = (x-sqrt(2))*(x^3+1)
        y = (x-sqrt(2))*(x*x+x+1);
        // yの値がdoubleの丸め誤差と同程度に小さければ終了
        if (y.get_value() < 1e-15) break;
        // ニュートン法の更新式で次のxを求める
        x = x - y.get_value()/y.get_d_value();
    }
    // x(=sqrt(2))を出力
    cout << x.get_value() << endl;
}

出力は以下の通りです。

1.41421356237309
BU_double(250, 100) 
1.4142135623731

出力の1行目は{\displaystyle \sin{x}}微分{\displaystyle x=\frac{\pi}{4}}のときの値の2倍(すなわち\sqrt{2})を表しています。最初に示した数値微分より高い精度で計算できていることが分かります。
2行目はy=10x^2という関数にx=5.0を入れて計算したときの関数値と導関数値です。for文を使った繰り返し文として定義しているにも関わらずきちんと導関数値が求まっていることが分かります
3行目はニュートン法を用いて求めた(x-\sqrt{2})(x^2+x+1)の解を表しています。ニュートン法f(x)=0の解を {\displaystyle x_{n+1} = x_n - \frac{f(x)}{f'(x)}}という更新式を用いて反復することで求めるというアルゴリズムです。導関数の計算が自動的に行われるため簡潔に記述することができます。


次回はトップダウン型自動微分を実装して実験してみようと思います。

参考文献

アルゴリズムの自動微分と応用 (現代非線形科学シリーズ)

アルゴリズムの自動微分と応用 (現代非線形科学シリーズ)

ここで使った自動微分クラスはこの本に書いてあるものをそのまま使っています。日本語で自動微分をテーマにした本はこの一冊しかなさそうです。Amazonの在庫は中古しかありませんが出版社のサイトから直接注文したら定価で入手することができました。在庫は少ないようです。


数値計算の常識

数値計算の常識

ニュートン法停止条件などで参考にしました。

広告コーナー