量子断熱時間発展(QAA)

量子断熱時間発展は、例えば組み合わせ最適化問題の解（最低エネルギー状態＝基底状態）を
求めるアルゴリズムです。
自明な問題の解から初めて、じわじわと系を解きたい問題にすり替えていくことで、求めたい問題の解へ収束させます。

QAAの実装は、QAOAの実装とほとんど同じですが、意外とサンプルコードがなかったので、まとめてみます。
本記事の貢献として、いろんなPython SDKで書いてみます。

解く問題は4ノードのMax cut です。
https://dojo.qulacs.org/ja/latest/notebooks/5.3_quantum_approximate_optimazation_algorithm.html

qiskit

qiskit実装として、ライブラリを使うものと、愚直に書くものがあります。

from qiskit.aqua.operators import WeightedPauliOperator
from qiskit import QuantumCircuit
from qiskit import execute, Aer
from qiskit.quantum_info.operators import Operator
from qiskit.tools.visualization import plot_histogram
import numpy as np

ライブラリ利用

n_qubits = 4
pauli_dict = {'paulis': [{"coeff": {"imag": 0.0, "real": 1}, "label": "XIII"}, {"coeff": {"imag": 0.0, "real": 1}, "label": "IXII"}, {"coeff": {"imag": 0.0, "real": 1}, "label": "IIXI"},{"coeff": {"imag": 0.0, "real": 1}, "label": "IIIX"}]}
H_ref = WeightedPauliOperator.from_dict(pauli_dict).to_opflow()

pauli_dict = {'paulis': [{"coeff": {"imag": 0.0, "real": 1}, "label": "ZZII"}, {"coeff": {"imag": 0.0, "real": 1}, "label": "IZZI"}, {"coeff": {"imag": 0.0, "real": 1}, "label": "IIZZ"},{"coeff": {"imag": 0.0, "real": 1}, "label": "ZIIZ"}]}
H_cost = WeightedPauliOperator.from_dict(pauli_dict).to_opflow()

T = 2
N = 8
dt = T/N
t = np.arange(0,T,dt)


circ = QuantumCircuit(n_qubits)
circ.barrier()
circ.h([0,1,2,3])
circ.barrier()
for tt in t:
    circ.hamiltonian(H_ref,(1-tt/T)*dt,[0,1,2,3])  # time evolution 
    circ.barrier()
    circ.hamiltonian(-1*H_cost,tt/T*dt*3,[0,1,2,3])  # time evolution 
circ.barrier()
circ.measure_all()

results = execute(circ,backend=Aer.get_backend('qasm_simulator'),shots=1024).result()
counts = results.get_counts()

# zero-padding
for i in range (2**n_qubits):
    counts.setdefault(format(i, '0'+str(n_qubits)+'b'),0)
plot_histogram(counts)

0101と1010が解です。

ちなみに実機では以下の通り。

辛い。

愚直に書く


circ = QuantumCircuit(n_qubits)
circ.barrier()
circ.h([0,1,2,3])
circ.barrier()
for tt in t:
    for i in range(n_qubits):        
        circ.rx(-2*(1-tt/T)*dt,i) #ref
    circ.barrier()

    for i in range(n_qubits):
        circ.cnot(i,(i+1)%n_qubits)  # time evolution 
        circ.rz(+2*tt/T*dt*3,(i+1)%n_qubits) #cost 
        circ.cnot(i,(i+1)%n_qubits)  # time evolution 
    circ.barrier()    

circ.barrier()
circ.measure_all()

braket

from braket.circuits import Circuit
from braket.devices import LocalSimulator
import numpy as np

n_qubits = 4
T = 2
N = 8
dt = T/N
t = np.arange(0,T,dt)

circ = Circuit()
circ.h(range(n_qubits))

for tt in t:
    for i in range(n_qubits):        
        circ.rx(i,+2*(1-tt/T)*dt) #ref

    for i in range(n_qubits):
        circ.cnot(i,(i+1)%n_qubits)  # time evolution 
        circ.rz((i+1)%n_qubits,-2*tt/T*dt*3) #cost 
        circ.cnot(i,(i+1)%n_qubits)  # time evolution 

device = LocalSimulator()
result = device.run(circ, shots=1024).result()
counts = result.measurement_counts
# 量子回路の実行結果を表示
# zero-padding
for i in range (2**n_qubits):
    counts.setdefault(format(i, '0'+str(n_qubits)+'b'),0)
print(counts)

Counter({'1010': 282, '0101': 244, '0011': 94, '0110': 80, '1001': 74, '1100': 68, '0100': 32, '1101': 25, '0001': 24, '1110': 20, '1011': 19, '0111': 18, '1000': 18, '0010': 16, '0000': 7, '1111': 3})

import matplotlib.pyplot as plt
plt.bar(counts.keys(), counts.values());
plt.xlabel('bitstrings');
plt.ylabel('counts');
plt.xticks(rotation=90);

出来ました。
braketの場合、RigettiやIonQにもかんたんに投げられます。
blueqat社のサービスを使うと、個人でアカウント登録→APIキーget→オンデマンド課金　がスムーズに出来ました。

PennyLane

ライブラリも使っちゃいますが。

import pennylane as qml
from pennylane import numpy as np
from pennylane.templates import ApproxTimeEvolution

T = 2
N = 8
dt = T/N
t = np.arange(0,T,dt)

n_wires = 4
wires = range(n_wires)

dev = qml.device('default.qubit', wires=n_wires)

coeffs = [1, 1, 1, 1]
op_ref = [qml.PauliX(0), qml.PauliX(1), qml.PauliX(2), qml.PauliX(3)]
op_cost = [qml.PauliZ(0) @ qml.PauliZ(1), qml.PauliZ(1) @ qml.PauliZ(2), qml.PauliZ(2) @ qml.PauliZ(3), qml.PauliZ(3) @ qml.PauliZ(0)]

hamiltonian_ref = qml.Hamiltonian(coeffs, op_ref)
hamiltonian_cost = qml.Hamiltonian(coeffs, op_cost)

@qml.qnode(dev)
def circuit():
    for i in range(n_wires):
        qml.Hadamard(wires=i)
    for tt in t:
        ApproxTimeEvolution(hamiltonian_ref, (1-tt/T)*dt, 1)
        ApproxTimeEvolution(-1*hamiltonian_cost, tt/T*dt*3, 1)
    return [qml.probs(wires = range(n_wires))]

probs = circuit()
probs = probs.numpy()
probs = probs.flatten()
probs

array([0.24672301, 0.02254855, 0.02254855, 0.07918272, 0.02254855, 0.00471733, 0.07918272, 0.02254855, 0.02254855, 0.07918272, 0.00471733, 0.02254855, 0.07918272, 0.02254855, 0.02254855, 0.24672301])

PennyLaneの場合、テンソル積はアットマークで表現しています。
probs の戻り値はテンソルなので、numpy配列に変換が必要です。

# プロットする
import matplotlib.pyplot as plt
%matplotlib inline

## z方向に射影測定した時に得られる可能性があるビット列
z_basis = [format(i,"b").zfill(n_wires) for i in range(probs.size)]

plt.figure(figsize=(10, 5))
plt.xlabel("states")
plt.ylabel("probability(%)")
plt.bar(z_basis, probs*100)
plt.show()

出来ました。
まぁPennyLaneでやるなら、こいつは最適化向きなので、QAOAをやったほうが良いですけど。

まとめ

いろんなPython SDK でQAA書いてみました。
やっていることはどれも同じですが、入出力の仕方などが微妙に違いますよね。

いろんなPython SDK で量子断熱時間発展(QAA)