量子位相推定(Quantum Phase Estimation: QPE)

pip install blueqat qiskit pylatexenc &>/dev/null

from blueqat import Circuit
import math

位相キックバックの分解能を上げる

前回に利用した位相キックバックアルゴリズム
を改良していきます。

位相キックバックは、位相シフト操作e^{i\theta}（ただし\theta=0 or \pi)が与えられた時、それを別の1量子ビットの|0 \rangle と |1 \rangleの位相差として移し取り、Hゲートの作用によりそれが0なのか\piなのか検出する1量子ビットのアルゴリズムでした。

これを\theta=0 or \piよりも高い分解能となるように拡張します。

(再掲) \theta=0 or \pi を検出する位相キックバック

Circuit().h[0].x[1].cu1(math.pi)[0,1].h[0].m[0].run_with_ibmq(returns='draw', output='mpl')

<Figure size 327.252x204.68 with 1 Axes>

具体的に \theta = \pi/4 の場合を考えます。

このとき、位相角を4倍すると、4\theta = \pi となります。
位相が\piに変換されていますので、位相キックバックができる形です。

つまり、位相シフト操作e^{i\theta}を4回行ったものに対して位相キックバックを行い、
得られた位相を後で（古典的に）1/4倍すればよいことになります。

同様に、位相角が\theta = 2\pi/2^{N}のときは、
位相シフト操作e^{i\theta}を2^{N-1}回行えばよいです。

問題は、位相角が半端な値で　 \theta = 2\pi(1/8+1/4+1/2) = 14\pi/8 のような時です。

この場合、少し工夫してやると1/8,1/4,1/2の部分を順繰り推定することができます。

1. まず、1/8の部分を取り出すために、位相シフト操作e^{i\theta}を2^2回数行います。

すると、4\theta = 2\pi(1/2+1+2) となります。

e^{i2n\pi}=1ですから、元の1/4,1/2の部分が何であったとしても、無視され、1/8の部分（つまり\pi/4)が位相\piに移されており、位相キックバックができる形です。

2. 次に1/4の部分を取り出します。

位相シフト操作e^{i\theta}を2^1回数行います。

2\theta = 2\pi(1/4+1/2+1) となります。

再びe^{i2n\pi}=1より、元の1/2の部分は無視できています。

取り出したい1/4の部分は位相\piに移されますので位相キックバックができる形です。

しかし、要らない1/8の部分が\pi/2という形で残ってしまっています。

だめじゃないか！　と思うわけですが、

1/8の部分というのは先程推定した位相なので、実は既知です。

なので、先程の推定結果が0、即ち0[rad]だった場合は0*\pi/2=0なので何もする必要がなく、

先程の推定結果が1、即ち\pi[rad]だった場合は、今は\pi/2という形で残っているはずですのでe^{-i\pi/2}を（量子ゲートとして）掛けることでキャンセルします。

3. 最後に1/2の部分を取り出します。

位相シフト操作e^{i\theta}を2^0回数行います。

1\theta = 2\pi(1/8+1/4+1/2) となります。

やはり 1/8および1/4の部分というのは先程推定した位相なので、既知です。

検出結果に応じてe^{-i\pi/4}およびe^{-i\pi/2}の項を相殺します。

このように、先に確定した桁の情報を使って古典的なif処理を行い、不要な部分を打ち消しながら順繰り1桁ずつ推定し、確定させていきます。

この操作により、位相は2\pi/2^{N}の精度で推定できます。

反復量子位相推定(Iterative QPE: IQPE)

上記のような方法で、決めた精度で位相の推定が可能となります。
応用例として、ユニタリ行列の固有値の位相推定に使うことが出来ます。

これを反復量子位相推定と言います。

量子位相推定

反復量子位相推定の場合、量子ビットは制御用に1bit,標的用に1bitの合計2bitあればよく、代わりにこの回路と測定をN度繰り返しています。

量子ビット数が一度にN+1個用意できる場合は、この操作を同時に（並列に）行うができ、たった1つの回路で終了します。

これを 量子位相推定 と呼びます。

量子位相推定の回路

blueqatで量子位相推定の回路を作ってみます。
この回路を覚える必要はなく、上述のような””考え方”がわかっていれば良いかと思います。

U = [1, 0 ; 0, e^{i\pi/4}] とし、入力する固有ベクトルは|1 \rangle, 対応する固有値の位相\pi/4=2\pi(1/8)を読み出してみます

def qft_dagger(circ, n):
    """n-qubit QFTdagger the first n qubits in circ"""
    # Don't forget the Swaps!
    for qubit in range(n//2):
        circ.swap[qubit, n-qubit-1]
    for j in range(n):
        for m in range(j):
            circ.cu1(-math.pi/float(2**(j-m)))[m,j]
        circ.h[j] 

def qpe(circ, n_iqft, theta):
  circ.h[0:n_iqft]
  repetitions = 1
  for counting_qubit in range(n_iqft):
    for i in range(repetitions):
        circ.cu1(theta)[counting_qubit, n_iqft] # This is Controlled-U        
    repetitions *= 2
  qft_dagger(circ, n_iqft)

circ_QPE = Circuit(4)
circ_QPE.x[3]
qpe(circ_QPE, 3, math.pi/4)
circ_QPE.m[:]
circ_QPE.run_with_ibmq(returns='draw', output='mpl', fold=30)

<Figure size 1772.2x325.08 with 1 Axes>

circ_QPE.run(shots=1024)

Counter({'1001': 1024})

'Counter({'1001': 1024})' と出ます。

'1001'のうち、右端の量子ビット'1'は、固有ベクトル入力|1 \rangleがそのまま出てきているだけです。

読み出し結果としては'100'で、左端の量子ビットが\pi/4=2\pi(1/8)の存在を意味しています。

その隣は、\pi/2=2\pi(1/4)、その隣は\pi=2\pi(1/2)の存在に相当します。

つまり　2\pi(1(1/8)+0(1/4)+0(1/2)) ということです。
少しUをかえて、U = [1, 0 ; 0, e^{i\pi/2}] とし、入力する固有ベクトルは|1 \rangle, 対応する固有値の位相\pi/2=2\pi(1/4)を読み出してみます。

circ_QPE = Circuit(4)
circ_QPE.x[3]
qpe(circ_QPE, 3, math.pi/2)
circ_QPE.m[:]
circ_QPE.run(shots=1024)

Counter({'0101': 1024})

Counter({'0101': 1024}) となります。つまり読み出し結果は010で、2\pi(0(1/8)+1(1/4)+0(1/2)) を意味しています。
以上、量子位相推定アルゴリズムの紹介でした。