解説1

機械学習研究とオリジナリティ

研究にはオリジナリティ（独自性）が欠かせない。他人の結果の追試は数件必要だが、それ以上は単に「車輪の再発明」である。機械学習、特にお金が絡む画像や動画の研究では、独自な結果を得るのは大変難しいと思う。大学や研究所に加えて、大小多くの企業が、優秀な人材と豊富な資金を投入している。かなりの研究結果が、学術論文誌ではなくweb上で自由に読めたり、プログラムやデータが手に入る。更に、誰でも参加可能なコンテストも開催されている。試しに “image video cup competition machine learning”で検索すると大量に見つかる。有名なのは Kaggle だが、ここには機械学習による予測コンテストが約250件ある。画像や動画が多いが、表情や身振りなども、6年前から行われている。有望な結果を参考に企業も改良をしているから、審査中を含めて多くの特許があると思われる。

研究を始める際には、これらを良く調べた上、例えば修士論文がきれいにまとまる2-3年後の状況を予想して、競争に勝てそうなら始める訳だ。画像や動画の研究は、非常に混んでいる上に理論が完成しつつあるので、私はやろうとは思わない。これらに比べると分子の機械学習はまだましだが、製薬会社を始め強力な競争相手もいる。私達が機械学習に物理や化学の知識を組み合わせようとする理由は、オリジナリティを得るためである。又、既にあるプログラムでデータを分析する研究なら、データをコンテストに出せば終わりではないか？データ分析の技術は急速にコモディティ化しそうだし、長期間価値を保つデータサイエンスの知識とは何か悩ましい。

先日聞いた修論発表にも、画像や動画の機械学習が数件あった。良いプログラムが公開や販売されていて簡単に手に入る上、過去の研究結果や手に入るプログラムとの比較も余りなされなかったので、これらの研究の価値に疑問を持った。

分子構造の生成器

ResNetなどの深い畳み込みニューラルネットで、画像分類が非常に上手にできるようになった頃、潜在変数から逆に画像を作り出す研究が盛んになった。敵対的生成ネット(GANs)が代表例であり、セレブ画像を学習して、現実にはいないセレブ画像を生成できたりする。プログラムとデータが公開されており、誰でも実験できる。

分子科学では長い間、望みの性質を持つ分子を設計したいと思ってきた。量子力学の発見以来「分子の構造が与えられた時、それがどんな性質か」予想する努力がされてきた。また現実的には量子力学がお手上げな複雑な系では、構造を変えると性質がどう変わるか実験で調べ、それを分析して設計をしてきた。大量のデータを機械学習できる時代になり「データを学習したモデルを逆転して、良い分子構造を出力する」可能性もでてきた。

そこで陸明浩さん（修士学生）とこのような方法を開発した。我々の方法では、分子構造を高次元特徴ベクトルで表現し、その高次元空間に適切な基底を取り、その部分空間のほぼ全ての格子点で、分子構造の原像を持つようにした。また特徴ベクトルを使い、教師データを良く学習できる事や、特徴ベクトルから分子構造を再構成するアルゴリズムを提案した。

ニューラルネットで発がん性を予想

2016年の秋に中国浙江省から、陸明浩さんが研究生として来日した。半年後の2月入試を受験し、情報学研究科に入学するためである。彼は動画から車線を検出する卒業研究をしていたので、畳み込みニューラルネットで化合物の性質を機械学習する事にした。画像や動画の機械学習は完成に近づいており、狭い分野でも世界一になるのは非常に難しい。それに対しグラフとして表される分子構造の機械学習には、当時決定版が無かった。我々は先行研究を参考に、畳み込みニューラルネットを用いる事にした。

各原子には近傍の原子番号を補い、分子の局所構造（官能基）を表した。これを並べて作った特徴ベクトルを左図のようなニューラルネットに入力した。芳香族ニトロ化合物の発がん性を教師データとして学習させた。畳み込み第1層で重要な局所構造を選択し、第2層でその組み合わせを選択することを狙った。

標準的方法で過学習を避け、交差検証で予測値を求めると、古典的構造活性相関の予測精度と同等以上の結果が得られた（左図参照）。構造活性相関では、疎水性を表すlogP値や軌道エネルギーで発がん性を予測したが、これら物理量も機械学習が可能な事を暗示している。他方、適切なpoolingで長距離相関を拾うことや、入力の置換に対して不変性を保証する事に課題が残った。

陸さんはこれらの結果をまとめ、大学院入試で試問を受け、良い成績で合格した。半年間我々が研究する間にも、幾つかの研究グループがグラフの機械学習法を学会やwebで報告している。機械学習分野では、学会のproceedingsやwebのarchivesに成果を発表し、伝統的な学術誌に発表しないことが多い。そのためweb of scienceなどの伝統的データベースで見つからない研究が多数あるし、計画-遂行-発表のサイクルが短く、人気分野の発展は非常に速い。

共鳴状態の量子シミュレーション 4

この方法で、独立電子（つまり軌道）模型の範囲で、共鳴状態を求めることができた。しかし、励起状態は占有軌道から作ったSlater行列式の線形結合で表され、軌道模型は荒い近似しか与えない。そこで、共鳴励起状態を簡潔に記述するため、我々は多体摂動論のBethe-Salpeter方程式を解くことにした。これは軌道模型に対応する無摂動1粒子Green関数から出発し、摂動効果を表すファインマン図形和を求めて、分極伝播関数を与える方程式である。我々の問題では、電子相関と境界条件効果という二重摂動が加わる事になる。これに問題を解き易くする色々な工夫をした。例えば、無摂動Green関数を密度汎関数理論から求めた軌道で作り、Bethe-Salpeter方程式の固有分極部分を交換相関核で代用し、境界条件効果を表す主要な自己エネルギーを残すなど。その結果、時間依存密度汎関数理論の状態に対する補正として、共鳴状態を求めることが出来た。

次にこの理論を、色素増感型太陽電池の電子移動過程に用いた。ペロブスカイト型太陽電池が出現する数年前までは、安価な太陽電池候補として色素増感型が熱心に研究され、1万件近くの実験結果が蓄積されていたので、理論検証に適していたためである。単体色素の理論計算の研究も500件以上あり、電子移動をシミュレートする新しい方法も数種類報告されていた。我々は（長い話を短くすると）過去の研究を参考にしつつ、anatase(101)面という酸化チタン表面に吸着した、有機ルテニウム(Ru)色素の良いモデルを探った。そして高性能色素の共鳴励起状態の性質を探った。

左図はblack dyeと呼ばれる色素が吸着した際の、吸光度と各励起状態 (単位femto秒)の寿命を示すが、これらは実験値と整合し、固体表面上の共鳴状態が上手にシミュレートできる事が分かった。また、色々な構造のRu色素で、構造と性質の関係を調べた。今回は実験が多い色素増感型太陽電池を例としたが、実験がない例で、理論の真価が発揮できるだろう。

図2: 共鳴状態を持つ1次元系

共鳴状態の量子シミュレーション 3

では固体表面に吸着した分子が、光エネルギーを吸収して励起状態になった時、どんな波動関数になるのか。この分子の「イオン化」エネルギーを、固体の伝導体に電子が入りうる最低エネルギーとしよう。分子がそれ以上のエネルギーを持つなら、ある速さで自動的にイオン化するが、これを超励起状態とみなそう。この場合の境界条件は明らかに発散球面波ではないし、孤立分子でよく使われる複素座標法を用いる根拠もない。

そこで我々は、適切な境界条件を持つ表面Green関数から複素数の光学ポテンシャルを求め、方程式にこの項を加えて境界条件を満たすようにした。これはメゾスコピック系の伝導度を計算するやり方と似ている。量子論が少し分かる皆さんに、我々のやり方を説明しよう。

簡単のため1次元系を考え、図2の左端の分子と右側の固体が障壁ポテンシャルで隔てられているとする。波動関数には左端でゼロ、右側で進行波となる境界条件を課すと、分子に局在化した準安定状態を表す。これを基底関数展開したり、サイト表示をすると、次式のHamiltonian Hとなる。

β は分子と表面の、t は固体原子間の飛び移り積分である。この行列を対角化するとシュレーディンガー方程式を解いたことになるが、有限の大きさの行列を使うと右端でゼロの境界条件を課すことになり、準安定状態が得られない事に注意しよう。表面Green関数を g(E)=(E−h)^−1 と定義し、ブロック行列の逆行列を求めると、g(E)が2次方程式 g=(E−c−tgt^T)^−1 を満たすことが分かる。

そしてg(E) を用いるとシュレーディンガー方程式は見かけ上2×2の方程式となる。赤の項がエネルギーに依存する複素数のポテンシャルで、自己エネルギーとか光学ポテンシャルと呼ばれる。g(E) は2次方程式の解なので2個あり、片方が右側で進行波となる解を与える。

共鳴状態の量子シミュレーション 2

私達はこの問題を、準安定な量子状態の波動関数を求める事で解決した。シュレーディンガー方程式を解く際には、波動関数が無限遠でゼロになるという孤立境界条件を課すか、格子定数の何倍かで出発点と同じ値になる周期境界条件を課す。それに対し準安定状態は、波動関数が無限遠で発散球面波となる境界条件を課す。これは有限寿命を持つ、崩壊していく系を表す波動関数である。これはGamovが原子核の研究に昔使ったもので、その後原子や分子の超励起状態にも使われてきたアイデアである。十分遠方でポテンシャルがゼロになるなら、発散球面波が散乱問題の解になる。

超励起状態とはイオン化エネルギーより高エネルギーの状態である。分子にイオン化エネルギーより高エネルギーの光を当てると、分子はまず光を吸収し励起状態となり、その後電子を放出する場合が見られる。これを量子力学で扱うには、時間依存シュレーディンガー方程式とか、非弾性散乱を解く方法が考えられるが、準安定状態が簡潔で扱いやすい。日本人では、慶応大学の薮下（元）教授や、分子科学研究所の江原教授が、分子の超励起状態の代表的な理論家である。この波動関数は規格化できない（二乗可積分でない）、bra-ketで期待値が得られないなど、普通の量子力学と違う点がある。

また波動関数に課す発散波境界条件も明らかではない。例えば水分子は10個の電子を持つので、30変数の偏微分方程式（シュレーディンガー方程式）を解き、それに適切なスピン関数をかけて足すと波動関数が得られる。イオン化する電子の座標に発散波境界条件を課す方法や、その関数を展開する基底関数の選び方などを考える必要がある。これまでの研究では、全電子座標に複素位相因子をかける複素座標法が良く使われる。

色素増感型太陽電池の構造：光励起した色素分子の電子が、酸化チタンに移動して電荷分離する。

共鳴状態の量子シミュレーション 1

固体表面上の分子などの波動関数には、最初は分子に局在しているが、やがて固体に広がっていく例がある。このような準安定状態の量子シミュレーションの方法を開発し、色素増感型太陽電池に応用した。この研究はポスドク研究員のDavid Sulzer博士や、名古屋大学情報学研究科の井内先生と共に行った。固体のバンドという、構造のある連続スペクトルに埋もれた波動関数を効率的に求める理論を、共鳴理論と多体摂動論から組み立てた。これまで孤立系に限られていた準安定状態の量子シミュレーションが、分子-固体複合系でも可能になり、光エネルギー変換過程の解明などに応用が可能である。またこの研究では、理論物理の色々な手法を組み合わせているので、物理好きの学生の皆さんにも楽しんで頂きたい。

全ての物質は量子力学に従うので、シュレーディンガー方程式を解くことで、物質の性質を予言したり、設計を試みたりできるはずだが、現在でも難しい対象が幾つかある。その例として広がった部分と小部分からなるもの、例えば固体表面上の分子や結晶中の欠陥などの、励起状態がある。この種の状態は光エネルギー変換でよく見られる。例えば色素増感型太陽電池では、透明な酸化チタン結晶表面に色素分子が吸着しており、この色素分子が光を吸収し、酸化チタン半導体に電子を渡すことで発電する。

酸化チタン表面の光学的性質は、バンド計算、つまり並進対称性下で密度汎関数理論を解く事、で理解できる。バンド計算では基底状態（最安定状態）で電子が入る占有軌道、及び空軌道と軌道エネルギーが分かり、軌道エネルギーは光学吸収を割と説明する。固体物理の基礎を学んでいれば、VASPやPWSCDなどのソフトで、バンド計算は手軽にできるようになった。厳密には密度汎関数理論の軌道エネルギーと光学吸収は異なり、時間依存密度汎関数理論を使うか、Green関数に基づくBethe-Salpeter方程式を解く必要があるが、詳細は省略する。

色素分子の光学的性質は、量子化学計算、つまり孤立系のシュレーディンガー方程式を解き励起状態の波動関数を求めるか、時間依存密度汎関数理論を解けば、励起エネルギーや光学吸収強度が分かる。量子化学の基礎を学んでいれば、GaussianやGAMESSなどのソフトで、この種の計算は簡単にできるようになった。ただ電子相関効果や溶媒効果など気を付ける点が幾つかある。

GPUによる量子化学シミュレーション 4

量子化学計算には、 Hartree-Fock交換項という計算が面倒な項が、しばしば必要になる。長代新治さん（当時修士学生）は、この複雑で難しい項を、超並列で計算するアルゴリズムを開発し、プログラムを実装した。先行研究と同等以上の性能が得られ、量子化学計算が相当速くなる事が分かった。この研究は2014-16年に行い、2016年アメリカ化学会年会などで発表済である。GPUによる高速プログラミングには、良いアルゴリズムと高度な実装技術が必要である。海外や国内で発表される研究には、比較用CPUプログラムがしょぼいだけという、見せかけだけの研究が多い。が本物のgeekを目指す学生の皆さんには、一度は超並列計算に挑戦して貰いたいと思う。

GPUによる量子化学シミュレーション 3

実用的計算ではd基底関数を使わなければならないが、最悪の場合、各threadは63個以下の変数を使い、約500個の積分を計算せねばならない。例えるなら9個の原材料から、2個を選んで色々な中間財を作り、原材料と作った中間財から別の中間財を作り、という手順を繰り返し、165個の最終製品を作る。この際、原材料や中間財の保管容器は30個程度しかない、という状況である。またある中間財を作るレシピは3種類あり、どれを使うか決める必要がある。

このスケジューリング問題をコンパイラーは上手に解けなかったので、丸岡寛典さん（当時学部4年生）は、良い計算手順を探索するプログラムを作成した。手持ちの中間財から作れる新しい中間財を調べ、それらの生産順序を全て調べる、というものである。この巨大な探索木を、枝刈りや確率などで小さく刈り込み、良さそうな解を見つけた。これらの解を元に十数種類のkernel関数を作成した。GPUでの計算速度を測定したところ、理論性能の40-70%が得られ、先行研究の5-10%の値を大幅に改善した。この結果を国際論文で発表した。

交換相関項は、排他原理と電子相関の効果を表すが、これは計算時間の4割を占める。交換相関項は電子密度の複雑な関数なので、3次元グリッドで求積する必要が有る。グリッド上での基底関数値の計算、グリッド上での電子密度の計算、交換相関ポテンシャルの行列要素の計算が必要である。無視できる程小さい項は計算しないので、不規則なデータが多く、SIMD processorで効率良く処理する方法が問題だった。細粒度並列性を確保するため、代表例のデータを解析し、original programを相当書き換えた。また短縮長や基底関数の数など可変長ループが複数あり、作業領域も足りない。色々なアイデアでプログラムを作り比較した。論文に書き難いノウハウが沢山あった。最終的にはCPUプログラムの5-10倍程度に速くなった。

密度汎関数理論の方程式。静電、交換相関、HF交換ポテンシャルの計算に時間がかかる。

GPUによる量子化学シミュレーション 2

密度汎関数計算には多数の複雑なステップがあるが、最も良く使われる500〜1000基底関数で計算時間を食う部分は、(1) 近距離Coulombポテンシャルの計算（40%）、(2) 交換相関ポテンシャルの計算（40%）、(2) Fock行列対角化（10%）なので、この部分のみ加速すれば良い。

この内streaming processorが苦手なのは(1)である。ここではGauss関数間の電子間反発積分を計算する。漸化公式を使い短縮Gauss関数を同時処理する方法が、理論計算量が最小で、汎用CPUでは最も高速である。この方法では1個の電子間反発積分を計算するのに、数千語の作業領域が必要だが、streaming processorには数十語の作業領域しかない。いろいろ検討した結果、direct J法で静電ポテンシャルを計算する時に限り、primitive Hermite Gauss 関数を使っても計算速度は下がらず、かつ作業領域も何とかなりそうだと分かった。

密度汎関数計算には多数の複雑なステップがあるが、最も良く使われる500〜1000基底関数で計算時間を食う部分は、(1) 近距離静電ポテンシャルの計算（40%）、(2) 交換相関ポテンシャルの計算（40%）、(2) Fock行列対角化（10%）なので、これらをまず加速する。

この内streaming processorが苦手なのは(1)である。ここではGauss関数間の電子間反発積分を計算する。漸化公式を使い短縮Gauss関数を同時処理する方法が、理論計算量が最小で、汎用CPUでは最も高速である。この方法では1組の電子間反発積分を計算するのに、数千語の作業領域が必要だが、streaming processorには数十語の作業領域しかない。検討した結果、Hermite Gauss基底で静電ポテンシャルを計算すれば、計算量が数割増えるだけで、作業領域も100語程度になりそうだと分かった。

2つめの問題は計算精度である。量子化学計算のは倍精度計算が必要だと言われていたが、描画装置であるGPUは、単精度計算しか速くない。そこで、大型分子で計算される電子間反発積分を調べると、ごく少数で大きな値の積分が、誤差の大部分を作る事が分かった。現在の積分計算アルゴリズムでは、Schwarz不等式を用いて積分上限値を求め、十分小さい値は計算しない。そこで単精度計算での誤差が無視できる項のみGPUで計算し、残りはCPUで計算した。9割の積分を単精度で計算しても、精度に問題が生じない事が分かった。積分計算の各ステップで必要精度も調べた。

3つめの問題はCPU-グラフィックボード間の通信の遅さである。平均すると数十演算で 1積分が計算されるため、これらをグラフィックボードから回収するには通信が遅すぎる。そこで静電ポテンシャル行列に足し合わせて回収した。これは古典2体力計算を加速するGRAPEシリーズと同じアイデアである。また実装上の問題として、十分な細粒度並列性を確保するため、ホストプログラムを相当修正する必要があった。またNVIDIAのコンパイラーの吐くkernel codeは最適化が不十分だったので、手探りで最適化した。

Page updated

Google Sites

Report abuse

解説1