「IBM Quantumで学ぶ量子コンピュータ」のQAOAのコードのアップデート

derwind

2023/05/05 17:02

目的

先日 blueqat club standard認定試験に向けての各種ポイントという記事が公開されました。この中に「QAOA」が含まれているのですが、恐らく書籍「IBM Quantumで学ぶ量子コンピュータ」に掲載されたコードは新しいQiskitでは動きません。Deprecated になってしまったパッケージを使っているからです。個人的な感想で恐縮ですが、QAOA 関係の論文

A Quantum Approximate Optimization Algorithm

XY-mixers: analytical and numerical results for QAOA

などを読んだだけでピンと来るかと言うとそうではないように思います。少なくとも私はそうでした。と言う事で、やはり手も動かしたいところかと思いますので、コードをQiskit 0.43.0 で動くところまでアップデートしたいと思います。

また、最近大規模アニーリングサービスTYTAN SDKの使い方などにあるように TYTAN SDK が開発されるなど、まだまだアニーリングの需要もあるようです。私はアニーリングは全然分かっておりませんので、今回併せて遊んでみたいと思います。

QAOA の実問題を解く (IBM Quantumで学ぶ量子コンピュータ pp.203-210 より)

まずは TYTAN SDK で解いてみる

ここからは書籍が手元にある前提とします。記号類も極力書籍に合わせます。まずはTYTAN SDKをインストールします。


!pip install -q git+https://github.com/tytansdk/tytan.git

チュートリアルは tytan_tutorial に色々ありますので、適宜参照したら良いと思います。日本語混じりのファイル名のやつが優しめな気がします。適当に真似をしながら実装をしてみます。


import tytan
import sympy as sym

q = sym.symbols("q0:4")
q

(q0, q1, q2, q3)


# コスト関数。
C = (q[0] + q[1] + q[2]) ** 2 + \
    q[3] ** 2 + \
    (q[0] + q[1] + q[2]) ** 2 + \
    q[0] ** 2 + \
    q[3] ** 2 + \
    (q[1] + q[2]) ** 2 + \
    q[0] ** 2 + \
    q[3] ** 2 + \
    (q[1] + q[2] + q[3]) ** 2

print(C)

2*q0**2 + 3*q3**2 + (q1 + q2)**2 + 2*(q0 + q1 + q2)**2 + (q1 + q2 + q3)**2

これだけだと、car1 が候補1と候補2を同時に通る、または両方とも通らないような解もでてきてしまうので、必ず片方通るように強制します。イジング変数 $\sigma_i \in \{-1, 1\}$ で考えると、 $\sigma_i \sigma_j$ のような項を加えると、イジング変数が同時に同じ値をとる時に正の値をとる、つまりコストが大きくなるような変数が作れます。イジング変数とバイナリ変数 $x \in \{0, 1\}$ の間には例えば $\sigma = 2x - 1$ のような関係が考えられるので、これを利用します。


# 車ごとに片方の道だけ通って欲しい。 
C_car_unique_path = (2*q[0] - 1) * (2*q[1] - 1) + (2*q[2] - 1) * (2*q[3] - 1)

weight = 3

C_total = C + weight * C_car_unique_path

C_total.expand(force=True)

という結果になりました。ところで、バイナリ変数なので $q_i^2 = q_i$ が成立しています。その処理を含めたほうが良いのでしょうか？念の為試してみましょう。


C1 = C_total.expand(force=True)
for i in range(4):
    C1 = C1.subs(q[i]**2, q[i])

C1

となりました。それぞれコンパイルしてみます。


Q, offset = tytan.qubo.Compile(C_total).get_qubo()

print(Q)
print(offset)

{('q0', 'q0'): -2, ('q1', 'q1'): -2, ('q2', 'q2'): -2, ('q3', 'q3'): -2, ('q1', 'q3'): 2, ('q0', 'q2'): 4, ('q1', 'q2'): 8, ('q2', 'q3'): 14, ('q0', 'q1'): 16}

6


Q1, offset1 = tytan.qubo.Compile(C1).get_qubo()

print(Q1)
print(Q1 == Q)
print(offset1)

{('q0', 'q0'): -2, ('q1', 'q1'): -2, ('q2', 'q2'): -2, ('q3', 'q3'): -2, ('q1', 'q3'): 2, ('q0', 'q2'): 4, ('q1', 'q2'): 8, ('q2', 'q3'): 14, ('q0', 'q1'): 16}

True

6

ログで出してしまいましたが、バイナリ変数を線形化してもしなくてもコンパイル結果には影響がないようです。ただ、ログで見える線形項と二次項の係数は線形化したQUBO式、つまり二番目の式のものになっているので注意が必要です。

QUBOへの定式化はできましたのでソルバを使って問題を解きたいと思います。


solver = tytan.sampler.SASampler()

result = solver.run(Q, shots=100)
print(result)

[[{'q0': 1.0, 'q1': 0.0, 'q2': 0.0, 'q3': 1.0}, -4.0, 79], [{'q0': 0.0, 'q1': 0.0, 'q2': 1.0, 'q3': 0.0}, -2.0, 14], [{'q0': 0.0, 'q1': 1.0, 'q2': 0.0, 'q3': 0.0}, -2.0, 7]]

と言う事で、書籍と同様に、car1は候補1を、car2は候補2を通ると良いそうです。

ちなみに、コード内では weight = 3 としていますが、これをもっと弱めた場合、例えば weight = 1 とした場合、全て 0 という解が得られます。これはどう言うことかについてはご自身で考えていただいたり、書籍 p.209 の上の部分を読んでいただいて、色々値を変更してみると面白いと思います。weight = 100 なども試すと面白いし、書籍に書いてある意味が分かると思います。

お次はQiskit (QAOA)で解いてみる

ここからが本命の「書籍のコードのアップデート」ではあります。まずは必要なモジュールをインストールします。


! pip install -q qiskit qiskit-optimization pylatexenc

基本的にはこのまま書籍の通りに進んで、書き方だけアップデートする形にしたいと思いますが、多少面倒臭かったところでは独自にヘルパー関数を定義したりもします。


from qiskit_optimization import QuadraticProgram
from typing import Tuple, Dict

LinearInfo = Dict[str, int | float]
QuadraticInfo = Dict[Tuple[str, str], int | float]

早速ここで、書籍にはなかったコードを追加します。折角TYTAN SDKをインストールしているので、これを使って線形項と二次項の係数の計算をさぼりたいと思います。コスト関数の部分にはcarが候補1か候補2のいずれかを通る事を強制する条件は含めません。量子回路でケアするからです。


q = sym.symbols("q0:4")
C = (q[0] + q[1] + q[2]) ** 2 + \
    q[3] ** 2 + \
    (q[0] + q[1] + q[2]) ** 2 + \
    q[0] ** 2 + \
    q[3] ** 2 + \
    (q[1] + q[2]) ** 2 + \
    q[0] ** 2 + \
    q[3] ** 2 + \
    (q[1] + q[2] + q[3]) ** 2

Q, offset = tytan.qubo.Compile(C).get_qubo()

def make_linear_and_quadratic(Q: dict) -> Tuple[LinearInfo, QuadraticInfo]:
    linear: LinearInfo = {}
    quadratic: QuadraticInfo = {}
    for key, coeff in Q.items():
        key0, key1 = key
        if key0 == key1:
            linear[key0] = coeff
        else:
            quadratic[(key0, key1)] = coeff
    return linear, quadratic

linear, quadratic = make_linear_and_quadratic(Q)
print(linear, quadratic)

{'q0': 4, 'q1': 4, 'q2': 4, 'q3': 4} {('q1', 'q3'): 2, ('q2', 'q3'): 2, ('q0', 'q1'): 4, ('q0', 'q2'): 4, ('q1', 'q2'): 8}

計算をさぼれました。再び書籍の書き方に戻ります。


qubo = QuadraticProgram()
qubo.binary_var('q0')
qubo.binary_var('q1')
qubo.binary_var('q2')
qubo.binary_var('q3')

qubo.minimize(linear=linear, quadratic=quadratic)
print(qubo)

minimize 4*q0*q1 + 4*q0*q2 + 8*q1*q2 + 2*q1*q3 + 2*q2*q3 + 4*q0 + 4*q1 + 4*q2 + 4*q3 (4 variables, 0 constraints, '')

上記はもうちょっと見やすいかもしれない以下のような出力もできます。この辺はお好みで。


print(qubo.export_as_lp_string())

\ This file has been generated by DOcplex

\ ENCODING=ISO-8859-1

\Problem name: CPLEX



Minimize

 obj: 4 q0 + 4 q1 + 4 q2 + 4 q3 + [ 8 q0*q1 + 8 q0*q2 + 16 q1*q2 + 4 q1*q3

      + 4 q2*q3 ]/2

Subject To



Bounds

続けて初期状態を用意します。API変更により、初期状態を用意する回路を定義すれば十分になりました。


from qiskit import QuantumCircuit

initial_state_circuit = QuantumCircuit(4)
initial_state_circuit.h(0)
initial_state_circuit.cx(0, 1)
initial_state_circuit.x(0)
initial_state_circuit.h(2)
initial_state_circuit.cx(2, 3)
initial_state_circuit.x(2)

display(initial_state_circuit.draw('mpl'))

この量子回路がどのような状態を用意するかというと

\begin{align*} |\psi_0\rangle = \frac{1}{2}(|0101\rangle + |0110\rangle + |1001\rangle + |1010\rangle) \end{align*}

です。続けて、QAOAを解いてみましょう。

以下、新しいQiskitを使っている場合warningが出るかもしれません。Qiskit Runtimeのプリミティブでアルゴリズム開発をより簡単になどにあるように、去年辺りから「Primitives」という概念を使ってサンプリングや期待値計算を楽にできるようにしようという流れになっているようです。そのため、APIなどが刷新されていっているのですが・・・今回、うまく噛み合わせきれなかったので、諦めて少し前のAPIで実装しています。Qiskitは結構APIが変わりやすいので、頑張ってください・・・。


from qiskit import Aer
# 旧版: https://github.com/Qiskit/qiskit-terra/blob/main/qiskit/algorithms/minimum_eigen_solvers/qaoa.py
from qiskit.algorithms import QAOA
from qiskit.opflow import I, X, Y
from qiskit_optimization.algorithms import MinimumEigenOptimizer

mixer = ((X^X^I^I) + (Y^Y^I^I)) / 2 + ((I^I^X^X) + (I^I^Y^Y)) / 2

backend = Aer.get_backend('aer_simulator')
step = 1
qaoa_mes = QAOA(quantum_instance=backend, reps=step, initial_state=initial_state_circuit, mixer=mixer)
qaoa = MinimumEigenOptimizer(qaoa_mes)
qaoa_result = qaoa.solve(qubo)
print(qaoa_result)

/tmp/ipykernel_1914/3785411303.py:11: DeprecationWarning: The class ``qiskit.algorithms.minimum_eigen_solvers.qaoa.QAOA`` is deprecated as of qiskit-terra 0.24.0. It will be removed no earlier than 3 months after the release date. Instead, use the class ``qiskit.algorithms.minimum_eigensolvers.QAOA``. See https://qisk.it/algo_migration for a migration guide.

  qaoa_mes = QAOA(quantum_instance=backend, reps=step, initial_state=initial_state_circuit, mixer=mixer)

fval=8.0, q0=1.0, q1=0.0, q2=0.0, q3=1.0, status=SUCCESS

ということで、アニーリング版と同様にcar1は候補1を、car2は候補2を通ると良いらしいことが分かりました。こちらもミキサー（ミキシングハミルトニアン）を渡さない場合にどうなるか考えてみて、書籍と照らし合わせると良いと思います。なお、アニーリングにしてもQAOAにしてもメタヒューリスティックなアルゴリズムですので、厳密解が得られるとは限らない・・・ことにご注意ください。

まとめ

単に書籍をそのまま書き直しただけですが、大体これでコードは動く状態にアップデートできると思います。念の為にQiskitのパッケージの各バージョンを出力しておきます。

QAOAに興味がわいた場合、最近ホットかもしれない？

Factoring integers with sublinear resources on a superconducting quantum processor

を読んでみるのも面白いかもしれません。


from qiskit import __qiskit_version__

__qiskit_version__

{'qiskit-terra': '0.24.0', 'qiskit-aer': '0.12.0', 'qiskit-ignis': None, 'qiskit-ibmq-provider': '0.20.2', 'qiskit': '0.43.0', 'qiskit-nature': None, 'qiskit-finance': None, 'qiskit-optimization': '0.5.0', 'qiskit-machine-learning': None}

「IBM Quantumで学ぶ量子コンピュータ」のQAOAのコードのアップデート

derwind

目的.css-1yssqf7{opacity:0.6;border:0;border-color:inherit;border-style:solid;border-bottom-width:1px;width:100%;margin-top:var(--chakra-space-2);}

QAOA の実問題を解く (IBM Quantumで学ぶ量子コンピュータ pp.203-210 より)

まずは TYTAN SDK で解いてみる

お次はQiskit (QAOA)で解いてみる

まとめ

目的