概要
量子計算の最も基本的なアルゴリズムの一つである量子フーリエ変換(QFT)をqiskitのモジュールでやろうとしたらエンディアンでハマったのでメモ。
参考文献
Addition on a Quantum Computer, Draper (1998)
(図もこの論文から転載)
量子フーリエ変換(QFT)
QFTは、抽象的には量子状態$|x\rangle = \sum_i x_i |i\rangle$ ($|i\rangle$は基底状態)に対して次のような状態$|\phi(x)\rangle = \sum_i \phi_i(x) |i\rangle$を返す計算です
$$
\phi_i (x) = \frac{1}{\sqrt N} \sum_{j=1}^N x_j e^{2\pi \sqrt{-1} ij/ N}
$$
通常 $N=2^n$ の形として、以下のような回路で実現されます。
ここで丸つきの数字$i$は、位相$2\pi/2^i$の(制御)位相ゲートです。
これを用いることで様々な量子計算が効率的に実現できるため、非常に基本的で重要なアルゴリズムとなっています。
実装
qiskitから以下のモジュールをインポート、また変数を定義
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, Aer
from qiskit.visualization import plot_state_qsphere
from qiskit.circuit.library import QFT
# データを二進法で表すための量子ビット数
data_qubits = 4
# データを表現する量子ビット
qr_data = QuantumRegister(data_qubits, "data")
# 量子回路
qc = QuantumCircuit(qr_data)
先ほどの図の通りに実装すると、次のようになります。
(番号付けが1~nと0~n-1でずれるので注意)
for i in range(data_qubits):
qc.h(qr_data[-1-i])
for j in range(i+1,data_qubits):
qc.cp(2 * math.pi / 2**(j+1-i), qr_data[-1-j], qr_data[-1-i])
また、これはqiskitのQFTモジュールを使って1 lineで書くこともできます...
qc = qc.compose(QFT(data_qubits))
が、ここでハマったというのが今回の話です
実はこのQFTモジュールの出力は上のべた書き実装とはエンディアンが異なります
そもそもqiskitのエンディアンの約束は若干紛らわしく、慣れるまで間違えることが多いと思います。
まずqiskitで量子ビットを宣言したときの順番とそれを測定したときのビット列の順番について以下の実装を見て先に整理しておきます。
data_qubits = 4
qr_data = QuantumRegister(data_qubits, "data")
qc = QuantumCircuit(qr_data)
qc.x(qr_data[0])
qc.measure_all()
backend = Aer.get_backend('qasm_simulator')
count = execute(qc, backend, shots=128).result().get_counts()
print(count)
qc.x(qr_data[0])によって0ビット目を反転しています。
このとき、測定結果は
{'0001': 128}
となることが確認できます。
つまり0ビット目は測定値の最下位ビットに対応していること(リトルエンディアン)がわかります。
先ほどのべた書き実装では、これを踏まえて上位ビットから順にHと回転ゲートをかけていったのでした。
しかし実際にはこうするとQFTモジュールの出力のエンディアンと逆になります。
これを見るために、QFTを使って足し算をするQFT加算器を考えます。
これは、QFTを行った後適切に位相シフトを行ってから逆QFTを行うことで足し算が実現できるというものです。
この図では $a = \sum_{i=1}^n a_i 2^{i-1}$ に $b = \sum_{i=1}^n b_i 2^{i-1}$ を(周波数空間上で)足す回路を表現しています。
これを使って先ほどの0ビット目の反転に対応するコードを書いてみます。
(上図で $b_1 = 1$ に対応)
data_qubits = 4
qr_data = QuantumRegister(data_qubits, "data")
qc = QuantumCircuit(qr_data)
### QFT ###
for i in range(data_qubits):
qc.h(qr_data[-1-i])
for j in range(i+1,data_qubits):
qc.cp(math.pi / 2**(j-i), qr_data[-1-j], qr_data[-1-i])
# 周波数空間で「1を加算」
for i in range(data_qubits):
qc.p(math.pi / 2**i, qr_data[i])
### 逆QFT ###
for i in range(data_qubits):
for j in range(i):
qc.cp(-math.pi / 2**(i-j), qr_data[i], qr_data[j])
qc.h(qr_data[i])
qc.measure_all()
backend = Aer.get_backend('qasm_simulator')
count = execute(qc, backend, shots=128).result().get_counts()
print(count)
出力：{'0001': 128}
ちゃんと「0に1を足す」ができました。
一方、このQFTをqiskitのQFTモジュールに置き換えると…
data_qubits = 4
qr_data = QuantumRegister(data_qubits, "data")
qc = QuantumCircuit(qr_data)
### QFT ###
qc = qc.compose(QFT(data_qubits))
# 周波数空間で「1を加算」
for i in range(data_qubits):
qc.p(math.pi / 2**i, qr_data[i])
### 逆QFT ###
qc = qc.compose(QFT(data_qubits, inverse=True))
qc.measure_all()
backend = Aer.get_backend('qasm_simulator')
count = execute(qc, backend, shots=128).result().get_counts()
print(count)
出力：{'1100': 26, '0100': 29, '1000': 55, '1001': 3, '1010': 10, '1011': 1, '0010': 3, '0101': 1}
観測値がばらけ、特に'1000'が最も多いのが分かります。
先の結果と同じにならなかったのは、QFT後のエンディアンが逆になっているため。
そこで途中の位相ゲートを
# 周波数空間で「1を加算」
for i in range(data_qubits):
qc.p(math.pi / 2**i, qr_data[-i-1]) # 上位ビットを大きく位相シフトする
と書き換えてあげます。
すると出力がちゃんと{'0001': 128}になりました。
まとめ
qiskitでQFTモジュールを使うときは出力のエンディアンに注意しましょう。
ちなみにwikipediaでQFTを検索したら、qiskitと同じエンディアンでの回路図が出てきたので、参照した論文の方向がマイナーだったのかもしれない…