Qiitaにポチポチと投稿しはじめていつのまにか1年たってました。定年後のセカンドキャリアを見据えて、狩猟免許を取ったり(なんで？)、量子エンジニア（アニーリング式）の認定試験を受けたりしてました。将来は猟ができるエンジニア、みたいな。

量子コンピューティングで線形方程式

簡単なものはいわゆる連立方程式です。もともとCAEに携わっていたこともあったので、大規模の方程式を解くことには興味がありました。量子コンピューティング見たときは、コレだ、っと思ったんですが、うまいアイデアが無くて放置してました。

最近、上記の認定試験を受けたこともあって、HOBOで3次以上にも対応できることがわかって、ちょっとアイデアを試すことにしました。例によって下手クソなソースはこちらです。

数値をどう表わすか

浮動小数の表現には、仮数部と指数部を使います。このうち指数部に量子ビットを使って$10^{q1+2*q2+4*q3}$ とかしてしまうと、ダメな感じ(本当にダメかどうかは確認してない)なのでどうしようか悩んだ末、以下のようにしました。

epart = 10**1 * q1 + 10**0 * q2 + 10**(-1) * q3
epart_onehot = (q1 + q2 + q3 -1)**2

10のべき乗を必要な範囲で全部書き出して、各量子ビットをワンホットで制約するという感じです。
仮数部は $\pm 1$ を必要な量子ビットで表現してもよかったのですが、とりあえずある範囲を適当な量子ビットで表現することにしました。合わせると以下な感じです。

create_double_qubo()

def create_double_qubo(dim: int, m_upper: float, m_lower: float, num_m_bit: int, digits_upper: int, digits_lower: int):
  m_bits = symbols_list([dim, num_m_bit], 'm{}_{}')
  e_bits = symbols_list([dim, digits_upper-digits_lower+1], 'e{}_{}')

  qdouble = []

  Honehot = 0

  # 仮数部を1.0に
  norm = 1.0/(2**num_m_bit)

  for row in range(dim):
    # 仮数部
    qm = 0
    for col in range(num_m_bit):
      qm += m_bits[row][col] * (2**col)
    qmantissa = lower + (upper-lower)*qm*norm

    # 指数部
    k = 1
    qe = 0
    ho = 0
    for col, ex in enumerate(range(digits_lower, digits_upper+1)):
      qe += e_bits[row][col] * (10**ex)
      ho += e_bits[row][col] # ワンホット

    #print('指数 H',qe)

    qdouble.append(qmantissa*qe)
    Honehot += k*(ho -1)**2

  return qdouble, Honehot

効果を確認する

数値を表現できれば、方程式を組んで解くだけです。

比較にあたり、適当な反復法で収束しにくい方程式を使ってみたかったのですが、ほどほどの大きさのものを見つけるのに苦労しそうだったので、ちょっと手を抜きました。
ヤコビ法をこちらで見つけたのと、初期値を設定できそうだったので利用させて頂きました。

初期値をゼロと量子コンピューティングで求めたものを、収束までに要した回数で比較しています。さらに、解が合っていることをnumpy と比較しています。

実際に実行してみると、反復回数がそこそこ減っています(130 → 78回)。簡単な未知数3個の方程式なので効果は感じられません。大規模になって、反復法では収束が難しいものでも、適切な初期値があれば収束に貢献できそうです。

実行環境

HOBOを使うにあたり、GPUがあると計算が速くなります。
google colaboratory でもできますが、ランタイムはT4にしてください。

おわりに

収束を早めるのに貢献できそうな初期値が得られそうだという感じがしました。適当な方程式を利用して GNU Octave あたりで比較してみたいと思っています。

おまけ

量子エンジニア（アニーリング式）。エントリーは合格しました。アドバンスは受けたんですがこの時点では未発表です。合格してるといいなぁ。