Qiskit.aquaでお手軽に増幅演算子Qを作る

以前のQAE(量子振幅推定)に関する記事
Qiskitで量子振幅推定(QAE)
では、頑張って増幅演算子を自作しました。
でも毎回やっているとめんどそうです。
そこでqiskit.aquaを使って手抜きで作ってみたいと思います。¹

qiskit.aquaでのコード

増幅演算子Qを作る

from qiskit import QuantumCircuit
from qiskit.aqua.algorithms import Grover

nqubit = 3

# initial state preparation 
state_preparation = QuantumCircuit(nqubit, name='State preparation')
state_preparation.h([1,2])
state_preparation.x(2)
state_preparation.ccx(2,1,0)
state_preparation.x(2)

# specify the oracle that marks the state '011' as a good solution
oracle = QuantumCircuit(nqubit)
oracle.z(0)

# define Grover's algorithm
grover = Grover(oracle=oracle, good_state=good_state, state_preparation=state_preparation)

これで回路が出来ています。

回路全体、あるいは一部（例えば$Q$)を取るには以下のようにします。
なお、グローバル位相$\pi$を勝手に入れてくれるので、教科書どおりの符号の振幅増幅演算子$Q$になってます。

grover.construct_circuit() #回路全体G=QA
grover.grover_operator.state_preparation #初期状態準備A
grover.grover_operator.oracle #オラクル部分の回路Sf
grover.grover_operator #増幅演算子Q

よって、$Q$の行列表示を取って古典的に固有値を出したいときは、

from qiskit.quantum_info.operators import Operator
circ_Q = grover.grover_operator
Unitary = Operator(circ_Q).data
a_eig = np.linalg.eig(+1*Unitary)
print("固有値\n {}\n".format(a_eig[0]))

でOKです。

こんな感じになります。

QへQPEをする

QPEで固有値を出したい場合、回路は以下の方法で手に入ります。

from qiskit.circuit.library import PhaseEstimation
circ_QPE = PhaseEstimation(num_evaluation_qubits=3,unitary=circ_Q)

ただ、どこを測定して、測定結果からどう位相を取り出すのか理解していないと使えません。
うーん・・・。。

続けます。
QPEをするには、回路に固有状態、またはその重ね合わせ状態を入れる必要があります。
上記の回路だと$|000 \rangle$が入っています。
指定した状態を入れるには、circ_QPE.initialize　のようにすればよいです。

例えば以下の回路だと、古典的に求めた固有ベクトルを入力として使います。(対応する固有値の規格位相は1/6です)

from qiskit.circuit.library import PhaseEstimation
num_evaluation_qubits = 3
circ_QPE = QuantumCircuit(num_evaluation_qubits+nqubit, num_evaluation_qubits)
circ_QPE.initialize(a_eig[1][:,1],range(num_evaluation_qubits,num_evaluation_qubits+nqubit))
circ_QPE.append(PhaseEstimation(num_evaluation_qubits=num_evaluation_qubits,unitary=circ_Q),range(nqubit+num_evaluation_qubits))
circ_QPE.measure(range(num_evaluation_qubits),range(num_evaluation_qubits-1,-1,-1))
circ_QPE.draw('mpl')

001 = 1/8 で、期待した1/6とまぁまぁ近いです。量子ビット数を増やすと1/6に漸近していきます。

もう一つの固有ベクトルを試してみます。initializeのときに　a_eig[1][:,2] を渡せばよいです。対応する固有値は5/6です。

111 = 1/2+1/4+1/8 = 0.875 です。5/6 = 0.833... にまぁまぁ近いですね。

一般には固有ベクトルは未知です。QAEではstate preparationのところで使った状態を入れるらしいです。
(実はこれは先の固有ベクトル2つの線形結合になっているんですが）

from qiskit.circuit.library import PhaseEstimation
num_evaluation_qubits = 3
circ_QPE = QuantumCircuit(num_evaluation_qubits+nqubit, num_evaluation_qubits)
circ_QPE.append(state_preparation,range(num_evaluation_qubits,num_evaluation_qubits+nqubit))
circ_QPE.append(PhaseEstimation(num_evaluation_qubits=num_evaluation_qubits,unitary=circ_Q),range(nqubit+num_evaluation_qubits))
circ_QPE.measure(range(num_evaluation_qubits),range(num_evaluation_qubits-1,-1,-1))
circ_QPE.draw('mpl')

はい、当然固有値も重ね合わせで出てきましたね。
あとは、$Q$の固有値の位相（QPEで求めたもの）が$2\theta/(2\pi)$であることに注意して、求めたい振幅$a=\sin^{2}\theta$を計算すればQAEが完了です。固有値は2つあるうちのどちらを入れてもOKです。なぜならば、

$a=\sin^{2}\theta = (\sin(\pi/6))^{2} = (\sin(5\pi/6))^{2} = 1/4 = 0.25 $

実際にQPEで（近似的に）求めたものを入れてみると

$a~(\sin(\pi/8))^{2} = 0.146$
または
$a~(\sin(0.875\pi))^{2} = 0.146$

です。

結論

qiskit.aquaは便利。²

1. 実は前の記事(QAA,Grover,QAEを勉強する。)でしれっとこの方法を使っています。 ↩

2. aquaにはQPEのライブラリもあります。が、受け付けるのがエルミート行列$H$のみになっているようでした。$H$の固有値を$U=e^{iH}$とおいて$U$へのQPEで取り出す流れを QPE と定義しているためではないかと思います。ただ、振幅推定のように$H$は与えられず$U$が与えられるケースも大変多いので、どうしたものか。 ↩