🌐

QiskitでQRAM付きのQuantum Countingを実装する方法

2021/09/20に公開
1

すでに量子計算の資料がかなりネット上に充実してきているので、quantum countingの実装はたくさん乗っていると思いますが、最近諸事情で改めてスクラッチで実装する機会があったので、僕なりの実装方法を共有します。

概要

今回は、Quantum Countingを使って「ステップbyステップで2つの選択肢から1つ選ぶ際に、選択した値の和がある定数以上の選択方法の数を数える」量子回路を実装します。

問題設定の具体的な状況としては、以下のような感じです。

list0 = [1,3,5,7]
list1 = [8,6,4,2]
list0, list1の対応するインデックスについて、いずれか一方から1つずつ選んで長さ4の新しいlistを作るとします
ステップ1: list0を選ぶ -> list0[0]を選択
ステップ2: list1を選ぶ -> list1[1]を選択
ステップ3: list1を選ぶ -> list1[2]を選択
ステップ4: list0を選ぶ -> list0[3]を選択
とすると、選び方は"0110"のように記述され、[1,6,4,7]のようなリストが作られます。
このリストの和は18になります。
このような選び方の全体は2^4 = 16通りあります。
一方、和が定数C=22より大きいような選び方は、"1100", "1110", "1011"の3通りのみです。

Quantum Countingについては、Qiskit Textbookの解説が割とわかりやすいと思うので、詳細はそちらに譲ります。
(僕は最初この教材で勉強しました、一部のコードはこの記事のものを改変しています: 使用箇所は以下に明記しています。)

また、実装について、量子回路のサブルーチンはqiskit.circuit.Gateにしてから使用することで、目視での確認を楽にすることができます(特にcontrolled gateの可視化)。
Gate型で使うメリットとして、他にも元のQuantumCircuitのqubit数やqubit配置に縛られず、コード上で楽に任意のqubitに対して任意のqubitの順序で量子操作を適用できる点や、楽にinverseを作成できる点が挙げられます。
具体的な使い方は、Qiskit Docs --- ControlledGateなどが詳しいです。

下準備

とりあえず、以下のライブラリをimportしました。

from typing import Union
import numpy as np
pi = np.pi

from qiskit import Aer
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, execute, assemble
from qiskit.compiler import transpile
from qiskit.circuit import Gate
from qiskit.circuit.library.standard_gates import ZGate
from qiskit.circuit.library import QFT
from qiskit.visualization import plot_histogram

量子定数加算減算回路

今回は加算、減算回路を連続的に何度も使用するので、Draper による、QFTを使用した定数加算手法を採用します。
量子定数加算減算回路は"さまざまな実装方法がありますが、古典の算術演算を量子回路上の可逆回路として実装したものは、carryのためだけの量子ビットを必要とし、それに加え、多くの制御ゲートを使用しています(詳しく調べてないけど、定数加算ならもっと賢い方法でできるのかな?)。

QFT AdderはQFTを行なってから、各量子ビットの位相空間をいじり、IQFTで戻す、というアイデアに基づいています。
これだけでは大きなメリットはなさそうですが(むしろ深さがO(n^2))、いくつもの算術演算が続く時に、ずっと位相空間で作業を行なって、最後にIQFTをすることで、一連の算術演算全体の深さや制御ゲート数の削減を目指すものです。

c.f.) QFT Adderは、Pavlidis and Gizopoulosの一般化手法やRuiz-Perez and Garcia-Escartinのサーベイとかが詳しいです。

定数加算減算についての位相操作は以下のような実装になります。


# Addition by QFT
def subroutine_add_const(num_qubits: int, const: int, to_gate=True) -> Union[Gate, QuantumCircuit]:
    const = const % (2 ** num_qubits)
    qc = QuantumCircuit(num_qubits)
    for i in range(num_qubits):
        for j in range(num_qubits - i):
            if const >> (num_qubits - 1 - (i + j)) & 1:
                qc.p(pi / (2 ** j), i)
                
    return qc.to_gate(label=" [+"+str(const)+"] ") if to_gate else qc
    
# Subtraction by QFT
def subroutine_sub_const(num_qubits: int, const: int, to_gate=True) -> Union[Gate, QuantumCircuit]:
    const = const % (2 ** num_qubits)
    qc = QuantumCircuit(num_qubits)
    for i in range(num_qubits):
        for j in range(num_qubits - i):
            if const >> (num_qubits - 1 - (i + j)) & 1:
                qc.p(-pi / (2 ** j), i)

    return qc.to_gate(label=" [-"+str(const)+"] ") if to_gate else qc

QRAMの実装方法

QRAMへのデータのエンコードにおいて、量子ビット数の多項式時間で完了する例は組み合わせ最適化問題のエンコードによく出てきます。
まず結論を言うと、今回の例は、以下のように実装できます。

list0 = [1,3,5,7]
list1 = [8,6,4,2]

def qram_encoder(index_qubits: int, data_qubits: int, list1: list, list2: list, to_gate = True) -> Union[Gate, QuantumCircuit]:
    
    qr = QuantumRegister(index_qubits, "index")
    qr_data = QuantumRegister(data_qubits, "data")
    # qr_anc = QuantumRegister(1, "ancillary")
    qc = QuantumCircuit(qr, qr_data)
    
    qc.append( QFT(len(qr_data), approximation_degree=0, do_swaps=False, inverse=False, name='QFT').to_gate(), qr_data[::-1] )
    for i, (val1, val2) in enumerate(zip(list1, list2)):
        qc.append(subroutine_add_const(len(qr_data), val1).control(1), qr[i:i+1] + qr_data[:])
        qc.x(qr[i])
        qc.append(subroutine_add_const(len(qr_data), val2).control(1), qr[i:i+1] + qr_data[:])
        qc.x(qr[i])
    qc.append( QFT(len(qr_data), approximation_degree=0, do_swaps=False, inverse=True, name='IQFT').to_gate(), qr_data[::-1] )
    
    return qc.to_gate(label=" QRAM Encoder ") if to_gate else qc

これを量子回路で出力してみると、

qc_qram_encoder = qram_encoder(index_qubits, data_qubits, list0, list1, to_gate = False)
qc_qram_encoder.draw("mpl")

qram_encoder

のようになります。(当然、最初に一様重ね合わせ状態を作ることが前提です。)

上図では、各量子ビットを、新しく作るlistの各indexにおいてlist0, list1のどちらを選択するかという情報に対応させています。

以上の回路が具体的に何をやっているかと言うと、一様重ね合わせ状態にある初期状態から、量子ビットi(つまり、ステップi)において、0を選ぶ場合にそのような選び方全体に対してlist0[i]を足しており、1を選ぶ場合にそのような選び方全体に対してlist1[i]を足しています。
こうすることで、n量子ビットのQRAMが作れます。
これを数式で表現すると、以下のようになります。

\sum_{x\in\{0,1\}^n} |x\rangle|f(x)\rangle\\ \text{ where } f(x) = \sum_{x_i\in x} \text{list}_{x_i}[i]

また、データの量子ビット数は、とりうる和の最大値の桁数をにします。
(以上の実装ではoracleの実装方法の関係で、実際には量子ビット数を1大きくしています。)

Oracle

Grover iterationのoracleとして、今回はQRAMにおいて、ある定数Cより大きい値のデータの位相を反転させるとします。
そのような印の付け方の一例として、carryに注目する方法があります。
具体的には、データの量子ビット数をdとすると、2^d - CをQRAMの全データに一律に加算してあげることで、定数Cより大きい値のみ最上位量子ビットに1が立ちます。
たとえば、今回の場合では、QRAMは次のようにエンコードされますが、

index data
0000 010100
0001 001101
0010 010001
0011 001010
0100 010101
0101 001110
0110 010010
0111 001011
1000 011001
1001 010010
1010 010110
1011 001111
1100 011010
1101 010011
1110 010111
1111 010000

これに2^5 - 22 = 9をたすと、

index data
0000 011101
0001 010110
0010 011010
0011 010011
0100 011110
0101 010111
0110 011011
0111 010100
1000 100010
1001 011011
1010 011111
1011 011000
1100 100011
1101 011100
1110 100000
1111 011001

となり、dataの最上位ビットが1になっているindexが3つあることが確認できます。
以上を実装したのが次の回路です。

def oracle(data_qubits: int, value: int, to_gate = True) -> Union[Gate, QuantumCircuit]:
    
    qr_data = QuantumRegister(data_qubits, "data")
    qr_anc = QuantumRegister(1, "ancillary")
    qc = QuantumCircuit(qr_data, qr_anc)
    
    qc.append( QFT(len(qr_data), approximation_degree=0, do_swaps=False, inverse=False, name='QFT').to_gate(), qr_data[::-1] )
    qc.append(subroutine_add_const(len(qr_data), value), qr_data[:])
    qc.append( QFT(len(qr_data), approximation_degree=0, do_swaps=False, inverse=True, name='IQFT').to_gate(), qr_data[::-1] )
    
    qc.cx(qr_data[0], qr_anc)
    
    # もし、ある定数「以下」のものに印をつけるなら、上のcxを消して次のようにすればいい。
    # qc.x(qr_data[0])
    # qc.cx(qr_data[0], qr_anc)
    # qc.x(qr_data[0])
    
    qc.append( QFT(len(qr_data), approximation_degree=0, do_swaps=False, inverse=False, name='QFT').to_gate(), qr_data[::-1] )
    qc.append(subroutine_sub_const(len(qr_data), value), qr_data[:])
    qc.append( QFT(len(qr_data), approximation_degree=0, do_swaps=False, inverse=True, name='IQFT').to_gate(), qr_data[::-1] )
    
    return qc.to_gate(label=" Oracle ") if to_gate else qc

この量子回路を出力すると、

qc_oracle = oracle(index_qubits, data_qubits, C_complement, to_gate = False)
qc_oracle.draw("mpl")

oracle

のようになります。
ancillaryの部分はQRAM付きGroverでお馴染み、phase kickback用の|-\rangle状態です。
oracleの外でこの状態に初期化してある設定なので、ここでは最初から|-\rangle状態だと思ってください。

実際にoracleとして使うには、エンコードしたQRAMを全て0状態に戻さなければいけないので、完全なoracleは次のようになります。
(実はqram_encoderoracleを一つの関数にすると、使用するQFTをさらに減らせます!)

def oracle_with_encoder(index_qubits: int, data_qubits: int, list1: list, list2: list, value: int, to_gate = True) -> Union[Gate, QuantumCircuit]:
    
    qr = QuantumRegister(index_qubits, "index")
    qr_data = QuantumRegister(data_qubits, "data")
    qr_anc = QuantumRegister(1, "ancillary")
    qc = QuantumCircuit(qr, qr_data, qr_anc)
    
    qc.append(qram_encoder(index_qubits, data_qubits, list1, list2), qr[:] + qr_data[:])
    qc.append(oracle(data_qubits, value), qr_data[:] + qr_anc[:])
    qc.append(qram_encoder(index_qubits, data_qubits, list1, list2).inverse(), qr[:] + qr_data[:])
    
    return qc.to_gate(label=" Oracle with Encoder ") if to_gate else qc

操作を反転するには、.reverse_ops()メソッドか、inverse()メソッドを使います。

qc_oracle_with_encoder = oracle_with_encoder(index_qubits, data_qubits, list0, list1, C_complement, to_gate = False)
qc_oracle_with_encoder.draw("mpl")

oracle_full

Diffusion

Grover iterationのdiffusionの部分は、みなさんおなじみの回路を使用します。

def diffusion(index_qubits: int, to_gate = True) -> Union[Gate, QuantumCircuit]:
    
    qr = QuantumRegister(index_qubits, "index")
    qc = QuantumCircuit(qr)
    
    qc.h(qr)
    qc.x(qr)

    qc.append(ZGate().control(index_qubits - 1), qr)

    qc.x(range(0, index_qubits))
    qc.h(range(0, index_qubits))
    
    return qc.to_gate(label=" Diffusion ") if to_gate else qc

図にすると、

qc_diffusion = diffusion(index_qubits, to_gate = False)
qc_diffusion.draw("mpl")

diffusion

Grover Iterator

Grover itarationは、先述のoracleとdiffusionを組み合わせ、
次のように実装できます。

def grover_iteration(index_qubits: int, data_qubits: int, 
                     list1: list, list2: list, value: int, 
                     to_gate = True) -> Union[Gate, QuantumCircuit]:
    
    qr = QuantumRegister(index_qubits, "index")
    qr_data = QuantumRegister(data_qubits, "data")
    qr_anc = QuantumRegister(1, "ancillary")
    qc = QuantumCircuit(qr, qr_data, qr_anc)
    
    qc.append(oracle_with_encoder(index_qubits, data_qubits, list1, list2, value), qr[:] + qr_data[:] + qr_anc[:])
    qc.append(diffusion(index_qubits), qr)
    
    return qc.to_gate(label=" Grover Iteration ") if to_gate else qc

全体の初期化

関数にするほどでもないですが、一応、以下のように実装できます。

def initialization(num_qubits: int, to_gate = True) -> Union[Gate, QuantumCircuit]:
    
    qc = QuantumCircuit(num_qubits + 1)
    
    qc.x(num_qubits)
    qc.h(range(num_qubits + 1))
    
    return qc.to_gate(label=" Initialization ") if to_gate else qc

Quantum Countingの大枠

Quantum Countingは一言で言うと、Grover iteratorを制御ユニたりにした位相推定です。

残念ながら、ここは指数個のGrover iterationを使用してしまいます。
上手い縮約方法あるのかな...あれば教えてください m(__)m

def quantum_counting(qft_qubits: int, index_qubits: int, data_qubits: int, 
                     list1: list, list2: list, value: int, 
                     to_gate = True) -> Union[Gate, QuantumCircuit]:
    
    qr_qft = QuantumRegister(qft_qubits, "qft")
    qr_index = QuantumRegister(index_qubits, "index")
    qr_data = QuantumRegister(data_qubits, "data")
    qr_anc = QuantumRegister(1, "ancillary")
    
    qc = QuantumCircuit(qr_qft, qr_index, qr_data, qr_anc)

    qc.append(initialization(qft_qubits + index_qubits), qr_qft[:] + qr_index[:] + qr_anc[:])
    
    iterations = 1
    for i_qft in range(qft_qubits):
        for i_iter in range(iterations):
            qc.append(grover_iteration(index_qubits, data_qubits, list1, list2, value).control(1), qr_qft[i_qft:i_qft + 1] + qr_index[:] + qr_data[:] + qr_anc[:])
        iterations *= 2

    qc.append( QFT(len(qr_qft), approximation_degree=0, do_swaps=True, inverse=True, name='IQFT').to_gate(), qr_qft)

    return qc.to_gate(label=" Quantum Counting ") if to_gate else qc

得られた確率分布は、以下の関数を使用して、位相から数(float)に変換します。
以下のコードはQiskit Textbookの解説の実装を多く参考にしました(ほぼ同じ)。

今回はQRAMの実装方法の解説なので、この部分の仕組みの詳細はQiskit Textbookの解説に譲ります。

# OK
def calculate_markers(measured_int, t, n):
    """
    measured_int: 
    t : qft_qubits (number of qft qubits)
    n : index_qubits (number of index qubits)
    """
    # Calculate Theta
    theta = (measured_int / (2 ** t)) * np.pi * 2
    print("Theta = %.5f" % theta)
    # Calculate No. of Solutions
    N = 2 ** n
    M = N * (np.sin(theta / 2) ** 2)
    print("Number of Marked Indices = %.2f" % (N - M))
    # Calculate Upper Error Bound
    m = t - 1 # Will be less than this (out of scope) 
    err = (np.sqrt(2 * M * N) + N / (2 ** (m - 1))) * (2 ** (-m))
    print("Error < %.3f" % err)

    return N - M

回路の実行

全体の回路の実行において、パラメタは次のように置きました。

list0 = [1,3,5,7]
list1 = [8,6,4,2]
C = 22 # 閾値
index_qubits = len(list0) # QRAMのインデックスの量子ビット数
max_l = sum([max(l0, l1) for l0, l1 in zip(list0, list1)]) # QRAMのデータがとりうる最大の値
data_qubits: int = np.ceil(mnp.log2(max_l)) + 1 # QRAMのデータの量子ビット数
C_complement: int = 2 ** (data_qubits - 1) - 1 - C # Oracleで足される定数
qft_qubits = 6 # 位相推定のQFTで使用する量子ビット数(推定精度を決定)

全体の実装についてお見せすると、以下のようになります。

qr_qft = QuantumRegister(qft_qubits, "qft")
qr_index = QuantumRegister(index_qubits, "index")
qr_data = QuantumRegister(data_qubits, "data")
qr_anc = QuantumRegister(1, "ancillary")
cr_qft = ClassicalRegister(qft_qubits, "c_qft")
qc = QuantumCircuit(qr_qft, qr_index, qr_data, qr_anc, cr_qft)
qc.append(quantum_counting(qft_qubits, index_qubits, data_qubits, list0, list1, C_complement), qr_qft[:] + qr_index[:] + qr_data[:] + qr_anc[:])
qc.measure(qr_qft, cr_qft)

(下図は上のコードの出力ではないです。)
quantum_counting

こちらを実行してみると、次のような結果になります。

t1 = time.time()
transpiled_qc = transpile(qc, simulator)
t2 = time.time()
print(t2 - t1, "s")

t1 = time.time()
qobj = assemble(transpiled_qc)
t2 = time.time()
print(t2 - t1, "s")

t1 = time.time()
hist = simulator.run(qobj).result().get_counts()
t2 = time.time()
print(t2 - t1, "s")

plot_histogram(hist)

77.32940173149109 s
1.2563889026641846 s
4.316287994384766 s
hist

(投稿するとなぜかラベルが全部黒塗りになってしまう...)

やはり重いのはtranspileなんですねぇ。
ちなみにtranspile後の回路の深さや基本ゲートの数はだいぶ多いです。

transpiled_qc.num_qubits, transpiled_qc.depth(), transpiled_qc.count_ops()

(17,
99572,
OrderedDict([('mcphase', 48510),
('ccx', 32319),
('cx', 8442),
('u', 7246),
('u1', 3717),
('p', 3529),
('cp', 708),
('u3', 188),
('h', 16),
('measure', 6),
('swap', 3),
('u2', 1)]))

以上の確率分布histから、最終的にC=22より大きい選び方の数は、次のように推定されます。

measured_int = int(max(hist, key=hist.get), 2)
num_markers = calculate_markers(measured_int, qft_qubits, index_qubits)
print(measured_int, num_markers)

Theta = 4.02517
Number of Marked Indices = 2.92
Error < 0.670
41 2.9248537266908325

このように、推定値はおおよそ3だと出力されました。

最後に: 感想や間違いの指摘、改善点の指摘などをいただけると泣いて喜びます!

Discussion

BOBOBOBO

これ、qram_encoderの部分が逆ですね(それに釣られて色々な部分が逆になる)
正しくは、val1とval2を入れ替えて、


def qram_encoder(index_qubits: int, data_qubits: int, list1: list, list2: list, to_gate = True) -> Union[Gate, QuantumCircuit]:
    
    qr = QuantumRegister(index_qubits, "index")
    qr_data = QuantumRegister(data_qubits, "data")
    # qr_anc = QuantumRegister(1, "ancillary")
    qc = QuantumCircuit(qr, qr_data)
    
    qc.append( QFT(len(qr_data), approximation_degree=0, do_swaps=False, inverse=False, name='QFT').to_gate(), qr_data[::-1] )
    for i, (val1, val2) in enumerate(zip(list1, list2)):
        qc.append(subroutine_add_const(len(qr_data), val2).control(1), qr[i:i+1] + qr_data[:])
        qc.x(qr[i])
        qc.append(subroutine_add_const(len(qr_data), val1).control(1), qr[i:i+1] + qr_data[:])
        qc.x(qr[i])
    qc.append( QFT(len(qr_data), approximation_degree=0, do_swaps=False, inverse=True, name='IQFT').to_gate(), qr_data[::-1] )
    
    return qc.to_gate(label=" QRAM Encoder ") if to_gate else qc

なぜなら、最初の設定では、val1が0のデータ、val2が1のデータを表すと決めたから。