common.title
Overview
Service overview
Terms of service

Privacy policy

Contact
Sign in
Sign up
common.title

量子ゲート方式インスパイアード、量子アニーリングシミュレータ

Yuichiro Minato

2021/01/14 12:59

#量子アニーリング #量子ゲート

はじめに

量子コンピュータのソフトウェア技術は日進月歩です。量子アニーリングという量子効果を活用しながら磁場でアニーリンをかける技術がありますが、量子ゲートの方ではそのアルゴリズム自体の改善が進んでいます。

その中でも現在各社から発売されている量子アニーリングや量子インスパイアードのアルゴリズムやハードウェアの弱点を補う新しいタイプのシミュレータを趣味で開発してみました。

課題

量子アニーリングはイジングモデルと呼ばれる物理モデルに特化した最適化アルゴリズムをハードウェアに起こした専用マシンですので、イジングモデルと呼ばれる物理モデルに特化した形で問題を設定します。アルゴリズムの上にアルゴリズムを作る感じです。

イジングモデルのアルゴリズムを作る上で設定できるのは、

１、局所磁場と呼ばれる量子ビット自体にかかる値
２、相互作用と呼ばれる量子ビット間にかかる値

で、これをQUBOやイジングハミルトニアンと呼ばれる式に落とし込んでいきます。問題を作っていく最中でやったことある人ならわかると思いますが、課題がいくつかあります。

１、多体相互作用という量子ビットが複数かかっている物は余剰量子ビットを割り当てて数式で分解する必要がある。
２、制約条件と呼ばれる満たすべき条件はソフトウェアで実装する必要があるソフト制約で実装する必要があり、面倒で解を間違える上、調整パラメータの調整が難しい。

それらを部分的に解決する量子インスパイアードアルゴリズムをイジングモデルで実装して実際に使えるものを作ってみたいと思います。

課題１：多体相互作用

こちらの解決は簡単です。wildqatのシミュレータはモンテカルロで動いているので、多体のままモンテカルロを動かして、エネルギーの評価をします。これをすることで、人間がわざわざ数式を作って苦労する必要がありません。ゲート方式の時間発展もわざわざ分解をしませんので。

素因数分解は複数の量子ビットの相互作用が出てきますが、これを分解せずにそのまま解いてみます。量子インスパイアードが可能なら、下記のような計算も可能ではないでしょうか。下記はイジングモデルで素因数分解をする際の数式になります。

H=128q0q1q2+16q0q1–56q0q2–52q0–48q1q2–96q1–52q2+196H = 128q_0q_1q_2 + 16q_0q_1 – 56q_0q_2 – 52q_0 – 48q_1q_2 – 96q_1 – 52q_2 + 196
import numpy as np
import copy

q = np.random.choice([0,1],3)
T = 5

while T>0.02:
  x_list = np.random.randint(0, 3, 1000)
  for x in x_list:
    q2 = copy.copy(q)
    q2[x] = 1-q[x]
    dE = 128*(q2[0]*q2[1]*q2[2]-q[0]*q[1]*q[2]) +16*(q2[0]*q2[1]-q[0]*q[1]) -56*(q2[0]*q2[2]-q[0]*q[2]) -48*(q2[1]*q2[2]-q[1]*q[2]) -52*(q2[0]-q[0]) -96*(q2[1]-q[1]) - 52*(q2[2]-q[2])
    if np.exp(-dE/T) > np.random.random_sample():
      q[x] = 1-q[x]
  T *= 0.99

print(q)
[0 1 1]

これで、答えが出ました。より一般化するには、

import numpy as np
import copy

#最初に定式化を行います。
def Hamiltonian(x):
  return 128*x[0]*x[1]*x[2] +16*x[0]*x[1] -56*x[0]*x[2] -48*x[1]*x[2] -52*x[0] -96*x[1] -52*x[2] +196

def run(H):

  #パラメータ初期化
  T = 100
  ite = 1000
  N = 3
  targetT = 0.02
  red = 0.95

  #量子ビットの初期化
  q = np.random.choice([0,1], N)

  while T>targetT:
    x_list = np.random.randint(0, N, ite)
    for x in x_list:
      q2 = copy.copy(q)
      q2[x] = 1-q[x]
      dE = H(q2) - H(q)
      if np.exp(-dE/T) > np.random.random_sample():
        q[x] = 1-q[x]
    T *= red

  return q

#とりあえず入れてrunすれば動きます
run(Hamiltonian)
array([0, 1, 1])

きちんと答えが出ました。

array([0, 1, 1])

ハミルトニアンを関数として定義して、runすると答えが出ます。これでなんでも解けそうです。タンパク質折りたたみでやってみます。

元となる文献は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

ハミルトニアンは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の後の数字は量子ビットの通し番号、各項の先頭の数字は係数です。

今回のシミュレータでは、QUBOをイジングに直さずにそのままコスト計算してしまっています。

def Hamiltonian(x):
  return -x[1] + 8*x[0]*x[1] + 15*x[1]*x[2] - 18*x[0]*x[1]*x[2] - 3*x[0]*x[3] + 12*x[0]*x[1]*x[3] + 4*x[2]*x[3] + 3*x[0]*x[2]*x[3] - 6*x[1]*x[2]*x[3] - 12*x[0]*x[1]*x[2]*x[3] + 4*x[1]*x[4] + 3*x[0]*x[1]*x[4] - 15*x[1]*x[2]*x[4] + 15*x[3]*x[4] + 3*x[0]*x[3]*x[4] - 6*x[1]*x[3]*x[4] - 12*x[0]*x[1]*x[3]*x[4] - 15*x[2]*x[3]*x[4] + 28*x[1]*x[2]*x[3]*x[4] - 2*x[0]*x[1]*x[5] - 4*x[2]*x[5] + 2*x[1]*x[2]*x[5] + 13*x[0]*x[1]*x[2]*x[5] - 2*x[0]*x[3]*x[5] + 4*x[0]*x[1]*x[3]*x[5] + 2*x[2]*x[3]*x[5] + 13*x[0]*x[2]*x[3]*x[5] + 4*x[1]*x[2]*x[3]*x[5] - 37*x[0]*x[1]*x[2]*x[3]*x[5] + 7*x[4]*x[5] + 2*x[1]*x[4]*x[5] + 13*x[0]*x[1]*x[4]*x[5] + 4*x[2]*x[4]*x[5] + 9*x[1]*x[2]*x[4]*x[5] - 33*x[0]*x[1]*x[2]*x[4]*x[5] - 20*x[3]*x[4]*x[5] + 13*x[0]*x[3]*x[4]*x[5] + 4*x[1]*x[3]*x[4]*x[5] - 37*x[0]*x[1]*x[3]*x[4]*x[5] + 9*x[2]*x[3]*x[4]*x[5] - 33*x[0]*x[2]*x[3]*x[4]*x[5] - 37*x[1]*x[2]*x[3]*x[4]*x[5] + 99*x[0]*x[1]*x[2]*x[3]*x[4]*x[5] - 4*x[1]*x[6] + 4*x[1]*x[2]*x[6] + 7*x[3]*x[6] + 2*x[1]*x[3]*x[6] + 13*x[0]*x[1]*x[3]*x[6] + 4*x[2]*x[3]*x[6] + 9*x[1]*x[2]*x[3]*x[6] - 33*x[0]*x[1]*x[2]*x[3]*x[6] + 4*x[1]*x[4]*x[6] - 18*x[3]*x[4]*x[6] + 9*x[1]*x[3]*x[4]*x[6] - 33*x[0]*x[1]*x[3]*x[4]*x[6] - 33*x[1]*x[2]*x[3]*x[4]*x[6] + 62*x[0]*x[1]*x[2]*x[3]*x[4]*x[6] + 7*x[5]*x[6] + 2*x[1]*x[5]*x[6] + 13*x[0]*x[1]*x[5]*x[6] + 4*x[2]*x[5]*x[6] + 9*x[1]*x[2]*x[5]*x[6] - 33*x[0]*x[1]*x[2]*x[5]*x[6] - 20*x[3]*x[5]*x[6] + 13*x[0]*x[3]*x[5]*x[6] + 4*x[1]*x[3]*x[5]*x[6] - 37*x[0]*x[1]*x[3]*x[5]*x[6] + 9*x[2]*x[3]*x[5]*x[6] - 33*x[0]*x[2]*x[3]*x[5]*x[6] - 37*x[1]*x[2]*x[3]*x[5]*x[6] + 99*x[0]*x[1]*x[2]*x[3]*x[5]*x[6] - 18*x[4]*x[5]*x[6] + 9*x[1]*x[4]*x[5]*x[6] - 33*x[0]*x[1]*x[4]*x[5]*x[6] - 33*x[1]*x[2]*x[4]*x[5]*x[6] + 62*x[0]*x[1]*x[2]*x[4]*x[5]*x[6] + 53*x[3]*x[4]*x[5]*x[6] - 33*x[0]*x[3]*x[4]*x[5]*x[6] - 37*x[1]*x[3]*x[4]*x[5]*x[6] + 99*x[0]*x[1]*x[3]*x[4]*x[5]*x[6] - 33*x[2]*x[3]*x[4]*x[5]*x[6] + 62*x[0]*x[2]*x[3]*x[4]*x[5]*x[6] + 99*x[1]*x[2]*x[3]*x[4]*x[5]*x[6] - 190*x[0]*x[1]*x[2]*x[3]*x[4]*x[5]*x[6]

def run(H):
  #パラメータ初期化
  T = 100
  ite = 5000
  N = 7
  targetT = 0.02
  red = 0.95

  #量子ビットの初期化
  q = np.random.choice([0,1], N)

  while T>targetT:
    x_list = np.random.randint(0, N, ite)
    for x in x_list:
      q2 = copy.copy(q)
      q2[x] = 1-q[x]
      dE = H(q2) - H(q)
      if np.exp(-dE/T) > np.random.random_sample():
        q[x] = 1-q[x]
    T *= red

  return q

run(Hamiltonian)
array([0, 0, 0, 1, 0, 1, 1])

結構計算量あるので時間かかりますが、そのうち計算時間については考えるとして、、、

array([0, 0, 0, 1, 0, 1, 1])

こちらも答えが出ました。

これまでは多くの問題で多体問題のシミュレーションは、数学分解によって多くの量子ビットが必要になるだけでなく、数式自体のパラメータ調整が新しく出てきてしまっていましたので、計算量が増えてしまうというトレードオフはありますが、大幅に人的なコストと時間の節約になります。

before:数学的な分解、数式のパラメータ調整が必要
after:分解不要、パラメータ調整不要

となりました。こちらは全く並列化されていないので、このまま並列化ができれば結構つかやすいのではないかと思います。

課題２：制約条件

イジングモデルを通常と臭いには、制約条件というものを設定する必要があるのですが、不幸なことにこちらがソフト制約と言って、ハミルトニアンの中に入れなくてはいけません。これは結構大きなネックです。

こちらも量子ゲート方式の量子コンピュータでは一部解決ができていますので、そちらの考え方を導入してみます。具体的にはN量子ビットからK量子ビットを選ぶような問題では制約条件を式の中に入れる必要がありません。拡張すれば巡回セールスマン問題も制約条件なしで解けそうです。

量子アニーリング方式では通常ランダムで量子ビットの値がスタートしますが、量子ゲート方式は全て0からスタートし、それから回路を使って値を制御していきます。ansatzと呼ばれる量子ビットの状態を事前に準備する段階で、指定の量子ビット数だけを1にすることができます。

さらに、量子ビット数を調整した上でswapゲートなどの量子ビット同士を入れ替える操作を時間発展という量子シミュレーションに組み込むことで、事前に準備した量子ビットの入れ替えをすることで、古典計算機のモンテカルロ方に対応したような操作を行うことができます。

参考：「[QAOA]SWAPの時間発展で|10>と|01>とを入れ替える」
https://qiita.com/gyu-don/items/c51a9e3d5d16a6d5baf6

ということで、巡回セールスマン問題のように、量子ビットのうち1になる数が決まっているようなものだったりは、swapゲートを適用することで大幅に利便性が上がる可能性があります。

https://arxiv.org/abs/1709.03489

ということで、ゲートインスパイアードで、アニーリングでも量子ビットをランダムにせず、入れ替えを許容することで制約条件を外せそうです。

H=∑i=0Nhixi+B(∑i=0Nxi−K)2H = \sum_{i=0}^N h_i x_i + B(\sum_{i=0}^N x_i -K)^2

というハミルトニアンを解いてみますが、前半はコスト関数と呼ばれるコストを最小化したいものを実装した式、後半のBがついた部分は制約条件項という量子ビット数を2に指定する制約条件を実装したものです。

今回作るシミュレータは後半のBがついた制約をなくそうというものです。

問題は、h0=1,h1=2,h3=3h_0=1,h_1=2,h_3=3としてみます。そして、量子ビットは1量子ビットだけを選ぶように条件を決めます。通常のアニーリングマシンのQUBOは、

H=h0+2x1+3x2+B(x0+x1+x2−1)2H = h_0 + 2x_1 + 3x_2 + B(x_0 + x_1 + x_2 -1)^2

を解きます（Bはハイパーパラメータなので探す）が、新しいシミュレータでは、

H=h0+2x1+3x2H = h_0 + 2x_1 + 3x_2

のみを解きます。今回工夫すべきは、これまでランダムで設定していた量子ビットをランダムから指定をします。

import numpy as np
import copy

#最初に定式化を行います。
def Hamiltonian(x):
  return 1*x[0] + 2*x[1] + 3*x[2]

def run(H):

  #パラメータ初期化
  T = 10
  ite = 1000
  N = 3
  targetT = 0.02
  red = 0.95

  #量子ビットの初期化
  q = [0,0,1]

  while T>targetT:
    x_list = np.random.randint(0, N, ite)
    for x in x_list:
      q2 = copy.copy(q)
      y = np.random.randint(0, N)
      q2[x] = q[y]
      q2[y] = q[x]
      dE = H(q2) - H(q)
      if np.exp(-dE/T) > np.random.random_sample():
        q[x] = q2[x]
        q[y] = q2[y]
    T *= red

  return q

run(Hamiltonian)
[1, 0, 0]

これを実行すると、

array([1, 0, 0])

上手い感じに解けました。念のためもう一つやってみます。

import numpy as np
import copy

#最初に定式化を行います。
def Hamiltonian(x):
  return 1*x[0] + 2*x[1] + 3*x[2] + 5*x[0]*x[1]

def run(H):

  #パラメータ初期化
  T = 10
  ite = 1000
  N = 3
  targetT = 0.02
  red = 0.95

  #量子ビットの初期化
  q = np.array([0,1,1])

  while T>targetT:
    x_list = np.random.randint(0, N, ite)
    for x in x_list:
      q2 = copy.copy(q)
      y = np.random.randint(0, N)
      q2[x] = q[y]
      q2[y] = q[x]
      dE = H(q2) - H(q)
      if np.exp(-dE/T) > np.random.random_sample():
        q[x] = q2[x]
        q[y] = q2[y]
    T *= red

  return q

run(Hamiltonian)
array([1, 0, 1])

答えは、

array([1, 0, 1])

となりました。うまく計算できています。

まとめ

このようにシミュレータを見てみました。

ということでできました。量子ゲート方式は2019年段階では制約条件を減らす方向にトレンドが進んでいますので、それを踏襲した内容を見直してみました。また進捗に合わせて色々実装してみたいと思います。

© 2024, blueqat Inc. All rights reserved