位数発見アルゴリズム ～Quantum Zooやっていく【特別編】～

Jan. 27, 2023, 2:50 p.m. edited Feb. 1, 2023, 2:54 p.m.

#Quantum Zoo  #量子情報  #量子力学 

$$ \def\bra#1{\mathinner{\left\langle{#1}\right|}} \def\ket#1{\mathinner{\left|{#1}\right\rangle}} \def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}} $$

Shor の素因数分解アルゴリズムや離散対数問題アルゴリズム、 Hallgren のペル方程式アルゴリズムでの根幹となっている位数発見アルゴリズムをサブルーチンとして用意しておくと便利になるので、一度実装しておくというもの。

位数発見アルゴリズムとは？

ある整数 \(x, N\in\mathbb{Z}\) (\(x<N\)) について

$$x^r\equiv 1\pmod N$$

を満たす最小の正の整数 \(r\) を位数という。 \(N\) を表現するビット数を \(L=\lceil \log N\rceil\) とすると、 \(\text{poly}(L)\) の時間計算量で位数を発見することのできる古典アルゴリズムは存在しない。一方、以下の量子アルゴリズムを用いると \(O(L^3)\) の時間計算量で求めることができる。

この位数発見アルゴリズムをベースとして、有名な Shor の素因数分解アルゴリズムなどを構築することができる。逆にいえば、 Shor の素因数分解アルゴリズムの中で古典と比べて唯一高速なのがこの位数発見アルゴリズムであり、ゆえに非常に重要なサブルーチンであるといえる。

量子アルゴリズム

詳しくは量子コンピュータと量子通信 II を読むのが一番手っ取り早い。変にググって時間を浪費するよりもよっぽど効率よく理解できる。つまり、このブログ読んでないでこういう専門書をまず読んで。このブログはここよりも Qiskit で実装することを目的としているので。

ということで、この本を読んだこと前提で簡潔に書くと、

$$\ket{u_s}=\frac{1}{\sqrt{r}}\sum_{k=0}^{r-1}e^{-2\pi isk/r}\ket{x^k\pmod N}$$

としたときに

$$\frac{1}{\sqrt{r}}\sum_{s=0}^{r-1}e^{2\pi isk/r}\ket{u_s}=\ket{x^k\pmod N}\tag{1}$$

が成り立つこと、および量子位相推定

$$\frac{1}{\sqrt{2^t}}\sum_{j=0}^{2^t-1}e^{2\pi ij\varphi}\ket{j}$$

↓量子フーリエ逆変換
$$\ket{\tilde{\varphi}}$$

（ただし、ここでは \(t,n\in\mathbb{Z}\), \(t>n\), \(\tilde{\varphi}\) は \(\varphi\) の \(n\) ビット近似）

を利用することで、以下のように表される：

1. \(t>2L+1\) を満たす \(t\) 個の量子ビットからなる第一レジスタ、 \(L\) 個の量子ビットからなる第二レジスタを用意： \(\ket{0^t}\ket{0^L}\)
2. アダマール演算により重ね合わせを作る： \(\frac{1}{\sqrt{2^t}}\sum_{j=0}^{2^t-1}\ket{j}\ket{0}\)
3. 冪剰余計算： \(\frac{1}{\sqrt{2^t}}\sum_{j=0}^{2^t-1}\ket{j}\ket{x^j\pmod N}\)
4. これは式 (1) を用いるとおおよそ \(\frac{1}{\sqrt{r2^t}}\sum_{s=0}^{r-1}\sum_{j=0}^{2^t-1}e^{2\pi isj/r}\ket{j}\ket{u_s}\) と式変形できる
5. 第一レジスタに量子フーリエ逆変換： \(\frac{1}{\sqrt{r}}\sum_{s=0}^{r-1}\ket{\widetilde{s/r}}\ket{u_s}\)
6. 第一レジスタを測定すると、 \(\widetilde{s/r}\) を得る

連分数展開により \(r\) を得る

前項で得られたものは \(\widetilde{s/r}\) であり、求めている \(r\) そのものではない。ここから \(r\) を得るために連分数展開を用いる。つまり、

$$\widetilde{s/r}=a_0+\frac{1}{a_1+\frac{1}{a_2+\frac{1}{\dots+\frac{1}{a_m+\frac{1}{\dots+\frac{1}{a_M}}}}}}$$

として、 \(0\leq m\leq M\) それぞれについて、

$$a_0+\frac{1}{a_1+\frac{1}{\dots+\frac{1}{a_m}}}=\frac{s'}{r'}$$

を計算して、 \(r'\) が正しい位数かどうかを \(x^{r'}\pmod N\) を計算して確かめていくというものである。時間がかかりそうに見えるが、個数としては高々 \(\text{poly}(L)\) なので問題ない。

また、測定して得られた \(\widetilde{s/r}\) 由来の \(r'\) に正しい位数が存在しない場合（\(\Omega(1/L)\) の確率で正しい位数は手に入るので単に繰り返すだけでもいける。これはそれすら厭うときのためのもの）、もう一度位数発見アルゴリズムを最初からやって別の \(r'\) を得る。そして最初の \(s'\) と今回の \(s'\) が互いに素なら 2 つの \(r'\) の最小公倍数が求める \(r\) となる。

（上述した本読んでないけどここまで読んで、なにいってんだ、となったら、やはり本読むべし。ちゃんと書いてある）

冪剰余計算

冪剰余の計算はあらかじめオラクルとして与えられているとしているが、実際にこれを実装するのが最も骨が折れる。世の中には、この回路での実装方法を論文としてまとめてくれている Quantum Networks for Elementary Arithmetic Operations が存在し、さらにはこれを日本語で解説してくださっている記事もある。これらを参考に Qiskit で実装した elementary.py を今回の実装例でも用いている。

Qiskit による実装

以下に Qiskit による実装を載せる。コードの全体は GitHub に置いてある。

from typing import Optional
from math import gcd

import numpy as np
import matplotlib.pyplot as plt
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, Aer, transpile, assemble
from qiskit.visualization import plot_histogram
from sympy import Rational
from sympy.ntheory.continued_fraction import continued_fraction, continued_fraction_convergents

from qft import qft
from elementary import ax_modM
from classical_utils import lcm

def order_finding(x: int, N: int, epsilon: Optional[float] = 0.2, show_hist: Optional[bool] = False) -> int:
    """Order-finding algorithm: it finds $r$ of $x^r\equiv 1\pmod N$. It requires 
    Args:
        x (int): $x$
        N (int): $N$
    Returns:
        order $r$
    """

    L = int(np.ceil(np.log2(N)))
    t = 2 * L# + 1 + int(np.ceil(np.log2(3 + 1 / (2 * epsilon)))) # 量子ビット数が多くてつらみなので epsilon 無視してる（けど動いてるのでヨシッ！）

    first_register = QuantumRegister(t)
    second_register = QuantumRegister(2 * t)
    auxiliary_register_mid = QuantumRegister(t)
    auxiliary_register_end = QuantumRegister(6 * t - 2)
    classical_register = ClassicalRegister(len(first_register))

    qc = QuantumCircuit(first_register, auxiliary_register_mid, second_register, auxiliary_register_end, classical_register)

    qc.h(first_register)

    #qc.append(ax_modM(a=x, M=N, N_len=len(first_register)), [first_register, auxiliary_register_mid, second_register, auxiliary_register_end]) # こんな感じで本当は QuantumRegister をつなげたいけど方法がわからぬ
    qc.append(ax_modM(a=x, M=N, N_len=len(first_register)), qc.qubits[:10 * t - 2])

    qc.append(qft(n=len(first_register)).inverse(), first_register)

    qc.measure(first_register, classical_register)

    backend = Aer.get_backend('aer_simulator_matrix_product_state')#('aer_simulator')
    qc = transpile(qc, backend)
    job = backend.run(qc, shots=10000)
    hist = job.result().get_counts()

    if show_hist:
        plot_histogram(hist)
        plt.savefig(f'img/order_finding_x{x}_N{N}_eps{epsilon}.png', bbox_inches='tight')

    all_fractions = []
    for measured_key, _ in sorted(hist.items(), key=lambda x: x[1], reverse=True):
        tilde_s_per_r = Rational(int(measured_key[-t:], 2), 2 ** t)
        if tilde_s_per_r == 0:
            continue

        fractions = []
        for fraction in continued_fraction_convergents(continued_fraction(tilde_s_per_r)):
            if pow(x, fraction.denominator, N) == 1:
                return fraction.denominator
            fractions.append(fraction)

            for other_fraction in all_fractions:
                if math.gcd(fraction.numerator, other_fraction.numerator) == 1:
                    r_candidate = lcm(fraction.denominator, other_fraction.denominator)
                    if pow(x, r_candidate, N) == 1:
                        return r_candidate

    raise Exception('r is NOT found!')


if __name__ == '__main__':
    #print(order_finding(x=5, N=21, show_hist=True)) # 大きくて時間がかかりすぎ
    print(order_finding(x=3, N=5, show_hist=True))

このコードを実行して \(x=3\), \(N=5\) の場合の \(r\) （答えは \(r=4\)）を求めようとすると、次のようなプロットが得られる：

4 本のビンが見えるが、左から順に \(0\), \(2^4=16\), \(2^5=32\), \(2^4+2^5=48\) が第一レジスタの測定により等確率で得られたということを表している。つまり、これらが表す \(\widetilde{s/r}\) はこれらを \(2^t=2^{2\times 3}=64\) で割った、 \(0\), \(16/64=1/4\), \(32/64=1/2\), \(48/64=3/4\) となる。これらのうち、正しく \(r=4\) を得られているのは 2 つのみである。 \(0\) は無視するとして、 \(2^5\) で失敗しているのは \(\tilde{s}\) と \(\tilde{r}\) が互いに素でなかったために約分され、 \(r\) の因数となってしまったためである。それでも、他の \(\tilde{r}\) との最小公倍数を計算すれば正しく \(r\) は求まるので問題ない。

よって、本アルゴリズムにより正しく \(r=4\) を得ることができた。以上。

素因数分解への応用

位数発見アルゴリズムを利用して素因数分解を古典アルゴリズムよりも指数関数的に高速に解くことができる。詳細はFactoring（ショアの素因数分解） 〜Quantum Zooやっていく〜に記述している。