@yuichirominato 2018.11.24更新 156views

タンパク質折りたたみ問題のイジング多体問題を効率的(?)にBlueqatで解く

Blueqat QAOA QUBO VQE イジング ザパタ タンパク質折りたたみ 固有値 量子ゲート 量子コンピュータ

はじめに

以前タンパク質折りたたみ問題の簡単な問題をイジングの量子アニーリングでといてみました。

D-Waveとwildqat.jsでタンパク質折りたたみ問題を解いて、アプリも作ってみた
http://blog.mdrft.com/post/414

解き方はタンパク質の折りたたみ方向を00,01,10,11の2量子ビットを利用したイジング実装によって組合せ最適化問題でときました。
元となる文献は2012年にハーバード大学のアランアスプルグジック先生の研究室からのNature論文でした。

Finding low-energy conformations of lattice protein models by quantum annealing
Alejandro Perdomo-Ortiz, Neil Dickson, Marshall Drew-Brook, Geordie Rose & Alán Aspuru-Guzik
Scientific Reports volume 2, Article number: 571 (2012)
https://www.nature.com/articles/srep00571

上記のブログでは簡単なタンパク質の折りたたみの例を取り出して最後の部分だけをD-Waveとウェブアプリ実装をしてみましたが、今回は6つの残基の折りたたみのハミルトニアンを頑張ってといてみたいと思います。ハミルトニアンは上記のNature論文のAppendixにあります。

ハミルトニアン

ハミルトニアンはQUBO表記されており、QUBOは最終的にイジングで解く{1,-1}ではなく、{0,1}のビット表現されています。そのビット表現をイジング表現に直す必要があります。

ハミルトニアンは下記です。

H = −q2 + 8q1q2 + 15q2q3 − 18q1q2q3 − 3q1q4 + 12q1q2q4 + 4q3q4 + 3q1q3q4 − 6q2q3q4 − 12q1q2q3q4 + 4q2q5 + 3q1q2q5 − 15q2q3q5 + 15q4q5 + 3q1q4q5 − 6q2q4q5 − 12q1q2q4q5 − 15q3q4q5 + 28q2q3q4q5 − 2q1q2q6 − 4q3q6 + 2q2q3q6 + 13q1q2q3q6 − 2q1q4q6 + 4q1q2q4q6 + 2q3q4q6 + 13q1q3q4q6 + 4q2q3q4q6 − 37q1q2q3q4q6 + 7q5q6 + 2q2q5q6 + 13q1q2q5q6 + 4q3q5q6 + 9q2q3q5q6 − 33q1q2q3q5q6 − 20q4q5q6 + 13q1q4q5q6 + 4q2q4q5q6 − 37q1q2q4q5q6 + 9q3q4q5q6 − 33q1q3q4q5q6 − 37q2q3q4q5q6 + 99q1q2q3q4q5q6 − 4q2q7 + 4q2q3q7 + 7q4q7 + 2q2q4q7 + 13q1q2q4q7 + 4q3q4q7 + 9q2q3q4q7 − 33q1q2q3q4q7 + 4q2q5q7 − 18q4q5q7 + 9q2q4q5q7 − 33q1q2q4q5q7 − 33q2q3q4q5q7 + 62q1q2q3q4q5q7 + 7q6q7 + 2q2q6q7 + 13q1q2q6q7 + 4q3q6q7 + 9q2q3q6q7 − 33q1q2q3q6q7 − 20q4q6q7 + 13q1q4q6q7 + 4q2q4q6q7 − 37q1q2q4q6q7 + 9q3q4q6q7 − 33q1q3q4q6q7 − 37q2q3q4q6q7 + 99q1q2q3q4q6q7 − 18q5q6q7 + 9q2q5q6q7 − 33q1q2q5q6q7 − 33q2q3q5q6q7 + 62q1q2q3q5q6q7 + 53q4q5q6q7 − 33q1q4q5q6q7 − 37q2q4q5q6q7 + 99q1q2q4q5q6q7 − 33q3q4q5q6q7 + 62q1q3q4q5q6q7 + 99q2q3q4q5q6q7 − 190q1q2q3q4q5q6q7

こちらの式の最低基底状態を求めてみたいと思います。ちなみにqの後の数字は量子ビットの通し番号、各項の先頭の数字は係数です。

解き方

今回はQAOAと呼ばれる汎用型量子コンピュータ向けの量子断熱計算のシミュレーションのようなアルゴリズムで組合せ最適化問題を解いてみたいと思います。QAOAに関しては下記を参照してください。

量子ゲートで組合せ最適化問題を解くQAOAの実装
http://blog.mdrft.com/post/229

QAOAではイジングのハミルトニアンをユニタリオペレータの時間発展シミュレーションを使って量子断熱計算の手法を利用して、|+>の重ね合わせの基底状態からスタートし、求めたいハミルトニアンに交換していくという手法で計算をします。

使うツール

使うのは手前味噌でMDRの配布している無料汎用型量子コンピュータ向けの開発キットSDKのBlueqat。pipで簡単にインストールできます。

https://github.com/mdrft/Blueqat

まずやること

まず最初の前処理が大変ですが、最初はQUBOをイジングに直します。この部分は社内のエンジニアの機転により自動化されました。今回はあまり深くは触れません。。。

このハミルトニアンを実際にQAOAでといていきます。時間発展の量子シミュレーションでは、各々のハミルトニアンの最小基底状態をVQEという手法で計算をします。また、QAOAでは精度を高めるために断熱量子計算のシミュレーションステップ数を決めることができます。

QAOAステップ数は8、ハミルトニアンの入れ替えに角度のパラメータが2あるので、合計2*8=16パラメータの最適化をします。最適化はBlueqatのVQEの呼び出す古典最適化アルゴリズムに依存したいと思います。実行コードは開発者により簡略化されており、


from blueqat.vqe import *
from blueqat.pauli import *
 
hamiltonian = 5.890625*I - 0.15625*Z[0]*Z[1] + 2.890625*Z[0]*Z[1]*Z[2] - 0.359375*Z[0]*Z[1]*Z[2]*Z[3] - 1.03125*Z[0]*Z[1]*Z[2]*Z[3]*Z[4] + 0.0625*Z[0]*Z[1]*Z[2]*Z[3]*Z[4]*Z[5] + 1.484375*Z[0]*Z[1]*Z[2]*Z[3]*Z[4]*Z[5]*Z[6] - 0.515625*Z[0]*Z[1]*Z[2]*Z[3]*Z[4]*Z[6] - 0.453125*Z[0]*Z[1]*Z[2]*Z[3]*Z[5] + 0.0625*Z[0]*Z[1]*Z[2]*Z[3]*Z[5]*Z[6] + 0.96875*Z[0]*Z[1]*Z[2]*Z[4] - 0.515625*Z[0]*Z[1]*Z[2]*Z[4]*Z[5]*Z[6] - 0.453125*Z[0]*Z[1]*Z[2]*Z[4]*Z[6] + 0.171875*Z[0]*Z[1]*Z[2]*Z[5] - 0.0625*Z[0]*Z[1]*Z[2]*Z[6] + 0.34375*Z[0]*Z[1]*Z[3] - 0.359375*Z[0]*Z[1]*Z[3]*Z[4] - 0.453125*Z[0]*Z[1]*Z[3]*Z[4]*Z[5] + 0.0625*Z[0]*Z[1]*Z[3]*Z[4]*Z[5]*Z[6] - 0.0625*Z[0]*Z[1]*Z[3]*Z[5] - 0.453125*Z[0]*Z[1]*Z[3]*Z[5]*Z[6] + 0.171875*Z[0]*Z[1]*Z[3]*Z[6] + 0.265625*Z[0]*Z[1]*Z[4] + 0.171875*Z[0]*Z[1]*Z[4]*Z[5] - 0.0625*Z[0]*Z[1]*Z[4]*Z[6] + 0.171875*Z[0]*Z[1]*Z[5]*Z[6] + 0.109375*Z[0]*Z[1]*Z[6] - 2.796875*Z[0]*Z[2] + 0.265625*Z[0]*Z[2]*Z[3] + 0.96875*Z[0]*Z[2]*Z[3]*Z[4] - 0.515625*Z[0]*Z[2]*Z[3]*Z[4]*Z[5]*Z[6] - 0.453125*Z[0]*Z[2]*Z[3]*Z[4]*Z[6] + 0.171875*Z[0]*Z[2]*Z[3]*Z[5] - 0.0625*Z[0]*Z[2]*Z[3]*Z[6] - 0.90625*Z[0]*Z[2]*Z[4] - 0.0625*Z[0]*Z[2]*Z[4]*Z[5] - 0.453125*Z[0]*Z[2]*Z[4]*Z[5]*Z[6] + 1.421875*Z[0]*Z[2]*Z[4]*Z[6] + 0.109375*Z[0]*Z[2]*Z[5] - 0.0625*Z[0]*Z[2]*Z[5]*Z[6] + 0.125*Z[0]*Z[2]*Z[6] - 0.28125*Z[0]*Z[3] + 0.265625*Z[0]*Z[3]*Z[4] + 0.171875*Z[0]*Z[3]*Z[4]*Z[5] - 0.0625*Z[0]*Z[3]*Z[4]*Z[6] + 0.171875*Z[0]*Z[3]*Z[5]*Z[6] + 0.109375*Z[0]*Z[3]*Z[6] - 0.171875*Z[0]*Z[4] + 0.109375*Z[0]*Z[4]*Z[5] - 0.0625*Z[0]*Z[4]*Z[5]*Z[6] + 0.125*Z[0]*Z[4]*Z[6] + 0.0625*Z[0]*Z[5] + 0.109375*Z[0]*Z[5]*Z[6] - 0.390625*Z[0]*Z[6] + 0.09375*Z[0] - 0.15625*Z[1]*Z[2] + 0.34375*Z[1]*Z[2]*Z[3] + 2.140625*Z[1]*Z[2]*Z[3]*Z[4] - 0.453125*Z[1]*Z[2]*Z[3]*Z[4]*Z[5] + 0.0625*Z[1]*Z[2]*Z[3]*Z[4]*Z[5]*Z[6] - 0.0625*Z[1]*Z[2]*Z[3]*Z[5] - 0.453125*Z[1]*Z[2]*Z[3]*Z[5]*Z[6] - 0.078125*Z[1]*Z[2]*Z[3]*Z[6] + 0.265625*Z[1]*Z[2]*Z[4] - 0.078125*Z[1]*Z[2]*Z[4]*Z[5] - 0.0625*Z[1]*Z[2]*Z[4]*Z[6] - 0.078125*Z[1]*Z[2]*Z[5]*Z[6] + 0.109375*Z[1]*Z[2]*Z[6] - 0.921875*Z[1]*Z[3] + 0.34375*Z[1]*Z[3]*Z[4] - 0.0625*Z[1]*Z[3]*Z[4]*Z[5] - 0.453125*Z[1]*Z[3]*Z[4]*Z[5]*Z[6] - 0.078125*Z[1]*Z[3]*Z[4]*Z[6] + 1.234375*Z[1]*Z[3]*Z[5] - 0.0625*Z[1]*Z[3]*Z[5]*Z[6] - 0.28125*Z[1]*Z[4] - 0.078125*Z[1]*Z[4]*Z[5]*Z[6] + 0.109375*Z[1]*Z[4]*Z[6] + 0.234375*Z[1]*Z[5] + 0.0625*Z[1]*Z[6] - 3.046875*Z[1] - 0.28125*Z[2]*Z[3] + 0.265625*Z[2]*Z[3]*Z[4] - 0.078125*Z[2]*Z[3]*Z[4]*Z[5] - 0.0625*Z[2]*Z[3]*Z[4]*Z[6] - 0.078125*Z[2]*Z[3]*Z[5]*Z[6] + 0.109375*Z[2]*Z[3]*Z[6] - 2.171875*Z[2]*Z[4] + 0.109375*Z[2]*Z[4]*Z[5] - 0.0625*Z[2]*Z[4]*Z[5]*Z[6] + 0.125*Z[2]*Z[4]*Z[6] + 0.0625*Z[2]*Z[5] + 0.109375*Z[2]*Z[5]*Z[6] + 0.359375*Z[2]*Z[6] + 0.09375*Z[2] - 0.28125*Z[3]*Z[4] + 2.671875*Z[3]*Z[4]*Z[5]*Z[6] + 0.109375*Z[3]*Z[4]*Z[6] - 2.515625*Z[3]*Z[5] + 0.0625*Z[3]*Z[6] - 0.671875*Z[3] + 0.0625*Z[4]*Z[5] + 0.109375*Z[4]*Z[5]*Z[6] - 2.390625*Z[4]*Z[6] + 0.21875*Z[4] + 0.0625*Z[5]*Z[6] - 0.203125*Z[5] - 0.125*Z[6]

result = Vqe(QaoaAnsatz(hamiltonian,8)).run()    
print(result.most_common(10))  

実行結果

実行結果は下記の通りです。計算は繰り返し計算とステップ数がありますのでかなり時間がかかりますが、16個のパラメータの最適化過程と、最終のハミルトニアンの期待値エネルギーが表示されます。


(略)

params: [0.88655361 0.41914504 0.96542813 0.20695551 0.77300442 0.12601914
 0.11899099 0.87285336 0.70352873 0.91151805 0.26923012 0.98264302
 0.91944822 0.3097855  0.88654577 0.74455706] val: 5.131692398138156

(略)

一定の最適化が終わると解が得られます。高速なシミュレータでないとなかなか大変ですが、数分から数十分かかることもあります。今回はHPCを使い、313秒で終わりました。状態ベクトルはギリギリ僅差で答えに到達しています。


(((0, 0, 0, 1, 0, 1, 1), 0.0526207005085115), ((0, 1, 0, 1, 0, 1, 1), 0.0509613913647074), ((0, 0, 1, 0, 1, 0, 0), 0.04083238349765797), ((0, 1, 0, 1, 0, 1, 0), 0.03858710306550349), ((1, 0, 1, 0, 0, 1, 0), 0.03797289965489858), ((1, 0, 0, 0, 0, 0, 1), 0.03159655384832289), ((1, 1, 1, 0, 1, 0, 0), 0.028299534909668748), ((1, 0, 0, 0, 0, 0, 0), 0.026228368428091498), ((1, 0, 1, 1, 1, 1, 0), 0.025016295668237147), ((1, 1, 1, 0, 1, 0, 1), 0.02392454173189753))
313.90521216392517

ちなみに答えは0100001011となり、最初の010は自明なので、0001011という配列が得られれば成功です。
QAOAは量子断熱計算で、エネルギーギャップの部分をきちんと通過しないといけないのと、初期のパラメータの設定次第で局所解に陥りやすくなるので、運もありますがある程度の精度で解くことができます。

あわせて読みたい

SERVICE

汎用量子コンピュータ向けアプリケーション開発SDK

詳しく見る Githubで入手する(無料)

汎用量子コンピュータ向け高速シミュレータ

詳しく見る

量子コンピュータ向けクラウドサービス(準備中)

詳しく見る

イジングマシン向けアプリケーション開発SDK

詳しく見る Githubで入手する(無料)

COMMUNITY

量子コンピュータのことを学びたい、仕事としたいなどの情報交換の場を設け、コミュニティの形成を進めています。オフラインの勉強会と、オンラインのチャットコミュニティの2種類あります。オフラインのConnpassは1400名、オンラインのSlackは880名を超える参加があります。どちらも無料ですのでお気軽にご参加ください。

CONNPASS SLACK

CONTACT

弊社あての連絡はinfo@mdrft.comより用件を明記の上、メールで連絡を頂けますと幸いです。

ブログトップに戻る

Recent Posts