化合物でも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を挙げると
- MACCS fingerprint [Durant+, 2002]
- Extended-connectivity fingerprints (ECFP) [Rogers+, 2010]
- Klekota-Roth fingerprint [Klekota+, 2008]
などがあります。 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は
という式で表されます。グラフ上で隣接するエッジからのメッセージを足し合わせるような処理をしていることが分かります。, などをうまく定めることで各種の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個以下しか持たない分子に限っても種類の分子が存在できるとされています。そのため全探索は不可能なので何らかの効率的な探索法を考える必要があります。
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のファインチューニングなども行っています。
- [Guimaraes+, 2017] はGANに強化学習を組み込むことで偏った性質を持つ分子のSMILESを生成しています
- [Yang+, 2017] ではモンテカルロ木探索と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))
入力
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言語に移植しました。
- 作者:Attila Szabo,Neil S. Ostlund
- 出版社/メーカー: 東京大学出版会
- 発売日: 1987/07/01
- メディア: 単行本
理論
本で200ページくらいかけて説明されている内容をブログに書くのは大変なので省略します。
ちなみに、以下のツイートの「100ページかけて説明してる別の本」はこの本を指します。
難しい内容を易しく簡潔に説明しようとするのは大変で、難しい内容は難しいまま長い時間をかけて説明する必要があるのだなぁという学びがありました。
ある本で30ページで説明されてて何も分からなかった内容を100ページで説明してる別の本読んだら完全に理解して何が分からなかったのか分からなくなった
— きばん (@kivantium) 2019年3月21日
実装
本の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) {
のように関数を定義しました。
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の関係で自動微分が脚光を浴びつつあるような気がしますが、自動微分を解説したページは少なくまだまだマイナーな分野だと思います。先日ようやく「アルゴリズムの自動微分と応用」を一通り眺めたのでいろいろ実験して遊んでいます。今日は自動微分のうち、ボトムアップ型自動微分をオペレータオーバーロードを用いて実現する方法について書きます。
自動微分とは
自動微分は、アルゴリズムによって定義された関数からその関数の偏導関数値を計算するアルゴリズムを導出するための技術です。一般的にはのような関数が先にあって、その関数を計算するアルゴリズムやプログラムがあるというように考えますが、自動微分ではどちらかというとアルゴリズムが先にあってアルゴリズムが表現する関数が生まれるというような考え方をします。
プログラムで微分を扱う上でよく知られている技術には数値微分と数式微分があります。
数値微分
数値微分は、小さな値を用いてやを計算することで微分の近似値を得る技術です。数値微分はあくまで微分の近似値を計算する技術ですが、自動微分ではアルゴリズムで定義された関数から解析的に得た偏導関数の値を計算します。また、数値微分では摂動幅を設定する必要がありますが、が大きすぎると近似精度が悪くなり、小さすぎると丸め誤差の影響で精度が悪くなるという問題があります。
この問題を理解するために簡単なプログラムを書きました。におけるの微分をを計算することで求めます。解析的にはが解なので、この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のときの結果が最も良いようです。
一般に、の計算誤差をなるべく小さくするためには関数の計算値に含まれる計算誤差をとしてとするのがよいことが知られています。しかし、を求めるのにの値が必要になってしまうのでこの式を使うのは現実的ではありません。
また、入力変数が多い関数に対しては入力変数の数だけ関数を計算する必要があり計算量が増えてしまいますが、自動微分を使うことで計算量を抑えることができます。
自動微分の種類
自動微分にはボトムアップ型(フォーワードモードとも)とトップダウン型(リバースモードとも)の2種類があります。どちらも微分の連鎖律を用いて関数を単純な演算に分解することで偏導関数値を求めますが、分解の方法が異なります。
ボトムアップ型自動微分
ボトムアップ型自動微分では一つの入力変数に対して、全ての中間変数の偏導関数値を計算していく手法です。
例としてWikipediaの自動微分のページにあるの例を考えます。この関数の計算をステップごとに中間変数に保存するとすれば、
のように5つの中間変数を用いて表すことができます。
ここで入力変数について各中間変数の微分を計算すると、
のように分解することができます。これは関数値を求める際に必要な各中間変数の計算と同時に行うことができます。
計算グラフを使って表すと以下のようになります。(Wikipediaから引用)
トップダウン型自動微分
トップダウン型自動微分では一つの出力変数について、全ての中間変数に対する偏導関数値を計算していく手法です。
先ほどと同じ例を考えます。出力変数について各中間変数に対する微分を計算すると、
のように分解することができます。
についての分解が分かりにくいですが、の変化は中間変数の変化を経由してしか起こらないことに注目すると納得できるかもしれません。計算グラフは以下の通りです。
トップダウン型自動微分では一度計算した後で各中間変数がどのように計算されたのかの履歴を逆向きに辿って偏導関数値を計算する必要があるため、計算グラフを保存するためのメモリが必要になります。
勾配を求める際には、入力変数の方が多いときはトップダウン型、出力変数の方が多い時はボトムアップ型を使うのが計算量の点で有利になります。ニューラルネットワークのバックプロパゲーションはトップダウン型の自動微分の一種です。
自動微分の実現方法
自動微分の実現方法には、ソースコードを解析して偏導関数値を求めるプログラムを作る方法とオペレータオーバーロードによる方法があります。ソースコードを解析する方法はコンパイラを作るのと同じような労力が必要になるので難しいですが、オペレータオーバーロードによる方法は比較的容易なので今回はこちらを採用します。
オペレータオーバーロードによるボトムアップ型自動微分
ボトムアップ型自動微分では関数値を計算するのと同時に偏導関数値を計算することができます。そのため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行目はの微分ののときの値の2倍(すなわち)を表しています。最初に示した数値微分より高い精度で計算できていることが分かります。
2行目はという関数にを入れて計算したときの関数値と導関数値です。for文を使った繰り返し文として定義しているにも関わらずきちんと導関数値が求まっていることが分かります
3行目はニュートン法を用いて求めたの解を表しています。ニュートン法はの解をという更新式を用いて反復することで求めるというアルゴリズムです。導関数の計算が自動的に行われるため簡潔に記述することができます。
参考文献
- 作者: 久保田光一,伊理正夫
- 出版社/メーカー: コロナ社
- 発売日: 1998/06
- メディア: 単行本
- クリック: 1回
- この商品を含むブログを見る
- 作者: 伊理正夫,藤野和建
- 出版社/メーカー: 共立出版
- 発売日: 1985/06/03
- メディア: 単行本
- 購入: 9人 クリック: 41回
- この商品を含むブログ (10件) を見る