量子アルゴリズムの基本：Shorのアルゴリズムの確認（その１）

$$
\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のアルゴリズム」です。アルゴリズムを説明した後、自作の量子計算シミュレータqlazyで、動作確認しようと思います。が、前段の説明がとても長くなったので、記事を２つに分けます。本記事（その１）では、アルゴリズムの説明をして、次回（その２）、シミュレータで動作確認してみます。

今回（および次回では）、以下の記事を参考にさせていただきました。実は、Shorのアルゴリズムについて解説した記事はここに書ききれないくらい沢山ありましたが、特に以下に掲載したものは良記事で、とても勉強になりました。＞著者のみなさま、ありがとうございました。

• Shorのアルゴリズム(MDR:加藤さん、湊さん、理研:中田さん)
• 量子位相推定とショアのアルゴリズム(@seiya-sugo)
• 数学好きから見た量子コンピュータ:57を因数分解した話(snuffkin)
• 量子コンピュータ（シミュレータ）で素因数分解する(@converghub)

アルゴリズムの説明

大まかな流れ

以下の[STEP.0]から[STEP.5]までの流れで、与えられた整数$N$の因数が求まります。よく「Shorのアルゴリズム」は「素因数分解アルゴリズム」と紹介されたりしますが、正確には、因数を求めるアルゴリズムです(因数分解アルゴリズム)。なので、素因数分解をするには、これを繰り返し適用する必要があります(RSA暗号解読という文脈では2つの素数の合成数を対象にするので、因数がわかったらその時点で素因数分解ができたということになりますが)。また、以下に説明する[STEP]をすべて量子コンピュータで実行するのではありません。大きな整数になると計算が非常に困難になる、[STEP.4]の位数を発見する計算だけを量子コンピュータでやらせます。それ以外の[STEP.0,1,2,3,5]は古典コンピュータでやる想定です。

[STEP.0] 素数判定

このステップは、Shorのアルゴリズムには通常含まれないと思いますが、与えられた整数が巨大な場合、そもそも因数が存在するのかどうか(つまり素数なのかどうか)わからないことが想定されるので、高コストな因数分解を実行する前にその判定をしておくべきです。一番簡単、かつ、おバカなのは、$2$から$N-1$までの整数で$N$を実際に割り算してみるという手法ですが、いろいろ高速なアルゴリズムが提案されています。

[STEP.1] 因数が自明な場合、計算してアルゴリズム終了

因数が存在するとわかったなら、それが簡単な因数である場合、出力してアルゴリズムを終了します。例えば、$N$が偶数であるかどうかはすぐわかるので(2で割ってあまりが0)、その場合「因数は2です」と言って、終了します。また、$N$がべき乗数、つまり$N = a^b$の形に書けるかどうかの判定および$a$と$b$を求めることも割とやさしくできるので、その場合「因数は$a$です」と言って、終了します。

[STEP.2] aをランダムに選ぶ

ここから、Shorのアルゴリズムに本格突入です。まず、$N$より小さい2以上の整数$a$をランダムに選びます。

[STEP.3] a,Nの最大公約数gcd(a,N)を計算する

$a$と$N$に共通した1ではない約数があるかどうか調べます。それには$a$と$N$の最大公約数$gcd(a,N)$を計算してみるのが一番簡単です(ユークリッド互除法という強力なアルゴリズムがあるので)。もし最大公約数が1以外であれば、「それが因数です」と言って、終了します。そうでない場合(つまり、$a$と$N$が互いに素である場合)、次のステップに行きます。

[STEP.4] a,Nに関する位数を求める(量子位数発見)

ここがShorのアルゴリズムの肝になる、量子コンピュータで実現するアルゴリズムです。何をやるかというと、

f(x) = a^x \space mod \space N

という関数の周期$r$を求めます。この$r$は初等整数論で言うところの「位数」に等しいので、この処理のことを「位数発見」と呼んでいます。位数は$a^x \space mod \space N = 1$を満たす最小の$x$と定義されています。ここで、「そもそもこの$f(x)$って周期関数なの？」という疑問の声が聞こえてきそうですが、ご安心ください。周期関数です。$a^y \space mod \space N = 1$が成り立つとき、$y$は位数$r$の整数倍であるということが証明できますし(証明略。背理法で簡単に証明できます)、また、$a^{x+r} \space mod \space N = a^x \space mod \space N$が成り立ちます。なので、$f(x)$は周期$r$で1が現れる周期関数になります。

このような位数$r$を「量子位数発見」で求めるのですが、その詳細は大事なところなので改めて後述します。

さて、この位数$r$がわかったとして、もしそれが奇数だとすると(以下のステップで都合が悪いので)、再度$a$を選び直す[STEP.2]に戻ります。偶数であれば、次のステップに行きます。実は、このとき偶数になる確率の方が大きいということがわかっているようで、なかなか次のステップに行けないということはないとのことです。

[STEP.5] 位数から因数を求める

位数$r$(偶数)に対して、

a^r \space mod \space N = 1

が成り立つので、適当な整数$m$をもってきて、

\begin{align}
&a^r = mN + 1 \\
&a^r - 1 = mN \\
&(a^{r/2}+1)(a^{r/2}-1) = mN \\
&(a^{r/2}+1)(a^{r/2}-1) \space mod \space N = 0
\end{align}

が成り立ちます。ということは、$(a^{r/2}+1)(a^{r/2}-1)$を$N$で割って割り切れるということなので、$(a^{r/2}+1)$と$(a^{r/2}-1)$の少なくともどっちか一方に$N$の因数が含まれるということになります。どうやったらその因数がわかるかというと、再び最大公約数です。$gcd(a^{r/2}+1,N)$もしくは
$gcd(a^{r/2}-1,N)$のどちらかが、$1$でも$N$でもない値をとったらば、「それが因数です」といって終了します。たまに失敗しますが、その場合は、[STEP.2]に戻ってもう一回$a$をランダムに選び直します。

以上が、「Shorのアルゴリズム」です。

量子位数発見

それでは、先程説明を省略した「量子位数発見」について見ていきます。

ユニタリ演算子の導入

量子位相推定の考え方を導入することで、位数を発見することを行っていくのですが、そのために、まず唐突なのですが、以下のように定義されたユニタリ演算子$U_a$を考えます($0 \leq a,y \leq N-1$、$a$と$N$は互いに素)。

U_a \ket{y} = \ket{ay \space mod \space N}

これがユニタリであることは、すぐにわかります。どんな$a$をとっても右辺のケットの中身は、$0$から$N-1$の値をとり、かつ、違う$y$に対して必ず違う値になります。なぜなら、もし、違う値の$y,y^{\prime}(0 \leq y < y^{\prime} \leq N-1)$に対して、右辺のケットの中身が同じになったとすると、

\begin{align}
&ay \space mod \space N = ay^{\prime} \space mod \space N \\
&a(y^{\prime}-y) \space mod \space N = 0
\end{align}

ということになり、これが成り立つためには、$y-y^{\prime}=0$か、$y-y^{\prime}=mN$($m$は整数)でなければならないため前提と矛盾しますので(背理法)。そうすると、上の変換$U_a$は入力された$0$から$N-1$の整数値を$0$から$N-1$の整数値のどれかに置換する演算(行列表現をイメージすると置換行列)となりますので、$U_a$は、ユニタリ演算子ということが言えます。

このユニタリ演算子$U_a$の固有状態$\ket{u_s}$は、

\ket{u_s} = \frac{1}{\sqrt{r}} \sum_{k=0}^{r-1} exp(-\frac{2\pi isk}{r}) \ket{a^k \space mod \space N}  \tag{0}

のように書けます($0 \leq s \leq r-1$)(ここになんとも不思議なことに位数$r$が登場します)。これが本当に固有状態になっているかどうかは、以下のように計算できることからわかります。

\begin{align}
U_a \ket{u_s} &= \frac{1}{\sqrt{r}} \sum_{k=0}^{r-1} exp(-\frac{2\pi isk}{r}) \ket{a^{k+1} \space mod \space N} \\
&= \frac{1}{\sqrt{r}} \sum_{k^{\prime}=1}^{r} exp(-\frac{2\pi is(k^{\prime}-1)}{r}) \ket{a^{k^{\prime}} \space mod \space N} \\
&= \frac{1}{\sqrt{r}} \sum_{k^{\prime}=1}^{r} exp(\frac{2\pi is}{r}) exp(-\frac{2\pi isk^{\prime}}{r}) \ket{a^{k^{\prime}} \space mod \space N} \\
&= \frac{1}{\sqrt{r}} \sum_{k^{\prime}=0}^{r-1} exp(\frac{2\pi is}{r}) exp(-\frac{2\pi isk^{\prime}}{r}) \ket{a^{k^{\prime}} \space mod \space N} \\
&= exp(\frac{2\pi is}{r}) \ket{u_s}
\end{align}

ここで、1行目から2行目へは、変数変換$k^{\prime}=k+1$を使いました。3行目から4行目へは、$a^0 \space mod \space N = a^r \space mod \space N$を使いました。最後の行からわかる通り、この固有値は、

exp(\frac{2\pi is}{r})

です。

位相推定のための量子状態準備

ユニタリ演算子の固有値に含まれる位相の中に位数$r$が入っているということで、これは以前の記事で説明した「量子位相推定」を召喚！、以上。と一瞬言いたくなるのですが、そう簡単には行きません。この場合、推定するのは$r$そのものではなく、$s/r$です。また、位相推定するためには、固有状態$\ket{u_s}$を準備する必要があるのですが、量子回路でどう実現すれば良いのかわかりません。というか、いま求めたい位数$r$がわからなければ固有状態をそもそも準備できないという罠ですね。

というわけで、一つの固有状態を用意するのは諦めます。代わりに、それらを重ね合わせたものを用意することにします。

\frac{1}{r} \sum_{s=0}^{r-1} \ket{u_s} = \ket{1}  \tag{1}

固有状態を均等に重ね合わせると、実は$\ket{1}$になるという面白い性質があります(ただし、$r$が偶数の場合です)。量子回路的に$\ket{1}$を用意するのは簡単なので、これを代わりに使おうというわけです。ちなみに上の等式が成り立つことは以下のように証明できます。

\frac{1}{r} \sum_{s=0}^{r-1} \ket{u_s} = \frac{1}{r} \sum_{s=0}^{r-1} \sum_{j=0}^{r-1} exp(-\frac{2\pi isj}{r}) \ket{a^j \space mod \space N} \tag{2}

ここで、右辺の$s$についての和、

\sum_{s=0}^{r-1} exp(-\frac{2\pi isj}{r})

に注目します。$j=0$の場合、和の中身は$1$なので、

\sum_{s=0}^{r-1} exp(-\frac{2\pi isj}{r}) = r \tag{3}

です。$j \ne 0$の場合は、

\begin{align}
&\sum_{s=0}^{r-1} exp(-\frac{2\pi isj}{r}) \\
&= \sum_{s=0}^{r/2-1} exp(-\frac{2\pi isj}{r}) + \sum_{s=r/2}^{r-1} exp(-\frac{2\pi isj}{r}) \\
&= \sum_{s=0}^{r/2-1} exp(-\frac{2\pi isj}{r}) + \sum_{s^{\prime}=0}^{r/2-1} exp(-\frac{2\pi i(s^{\prime}-r/2)j}{r}) \\
&= \sum_{s=0}^{r/2-1} exp(-\frac{2\pi isj}{r}) + \sum_{s=0}^{r/2-1} exp(-\frac{2\pi isj}{r}) exp(\frac{2\pi r/2}{r}) \\
&= \sum_{s=0}^{r/2-1} exp(-\frac{2\pi isj}{r}) + \sum_{s=0}^{r/2-1} exp(-\frac{2\pi isj}{r}) exp(i\pi) \\
&= \sum_{s=0}^{r/2-1} exp(-\frac{2\pi isj}{r}) - \sum_{s=0}^{r/2-1} exp(-\frac{2\pi isj}{r}) \\
&=0 \tag{4}
\end{align}

ということで$0$になります。ここで、2行目から3行目へは、$s^{\prime}=s+r/2$という変数変換を行いました。ちょっと面倒な式展開をしてしまいましたが、上の和は複素空間上で考えると、単位円に内接する正$r$角形の頂点の和に相当するので、$0$であることは直感的にわかると思います。これら(3)式と(4)式を合わせると、

\sum_{s=0}^{r-1} exp(-\frac{2\pi isj}{r}) = r \delta_{j,0} \tag{4}

となり、(2)式に代入すると、

\begin{align}
&\frac{1}{r} \sum_{s=0}^{r-1} \ket{u_s} = \frac{1}{r} \sum_{s=0}^{r-1} \sum_{j=0}^{r-1} exp(-\frac{2\pi isj}{r}) \ket{a^j \space mod \space N} \\
&= \frac{1}{r} \sum_{j=0}^{r-1} r \delta_{j,0} \ket{a^j \space mod \space N} \\
&= \ket{1}
\end{align}

となります。これで(1)式が証明できました。

というわけで、話を戻します。量子位相推定を行うに際し、固有状態を一つ用意する代わりに、それらを重ね合わせたものを用意します、という話をしていました。

\frac{1}{r} \sum_{s=0}^{r-1} \ket{u_s} = \ket{1}  \tag{1}

ですね。以前の記事で説明した量子位相推定では、$\ket{\psi}$を固有状態として、補助ビットと合わせ、$\ket{0} \otimes \ket{\psi}$というテンソル積にしておいた上で、$\ket{0}$にアダマールをかけて、

\frac{1}{\sqrt{2^n}} \sum_{k=0}^{2^{n}-1} \ket{k} \otimes \ket{\psi} \tag{5}

という状態にしておき、それに、いま考えているユニタリ演算子を制御ユニタリ演算子にして、それをズラズラ適用して、最後に補助量子ビットの方に逆量子フーリエ変換を施すのでした。今回はちょっと違っていて、(5)式の代わりに、

\begin{align}
&\frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^{m}-1} \ket{k} \otimes \ket{1} \\
&= \frac{1}{\sqrt{2^m}} \sum_{k_0=0,1} \ket{k_0} \cdots \sum_{k_{m-1}=0,1} \ket{k_{m-1}} \otimes \ket{1} \tag{6}
\end{align}

を用意します。これに対して、今回の場合、

\Lambda_{0}(U_a^{2^0}) \Lambda_{1}(U_a^{2^1}) \cdots \Lambda_{m-1}(U_a^{2^m-1}) \tag{7}

という演算を適用します。ここで、$\Lambda_{i}(U)$は$i$番目のビットを制御ビットにした制御$U$演算子を表すものとします。$m$は補助量子ビット数です。$m$を大きくとれば、その分推定精度はあがります。

さて、それでは、(6)式に(7)式の演算子を適用してみます。

\begin{align}
&\frac{1}{\sqrt{2^m}} \sum_{k_0=0,1} \Lambda_{0}(U_a^{2^0}) \ket{k_0} \cdots \Lambda_{m-1}(U_a^{2^{m-1}}) \sum_{k_{m-1}=0,1} \ket{k_{m-1}} \otimes \ket{1} \\
&= \frac{1}{\sqrt{2^m}} (\ket{0}+\ket{1}U_a^2) \cdots (\ket{0}+\ket{1}U_a^2) \otimes \ket{1} \\
&= \frac{1}{\sqrt{r}} \sum_{s=0}^{r-1} \frac{1}{\sqrt{2^m}} (\ket{0}+\ket{1}U_a^2) \cdots (\ket{0}+\ket{1}U_a^2) \ket{u_s} \\
&= \frac{1}{\sqrt{r}} \sum_{s=0}^{r-1} \frac{1}{\sqrt{2^m}} (\ket{0}+exp(\frac{2\pi is}{r}\cdot 2)\ket{1}) \cdots (\ket{0}+exp(\frac{2\pi is}{r}\cdot 2)\ket{1}) \ket{u_s} \\
&= \frac{1}{\sqrt{r}} \sum_{s=0}^{r-1} \frac{1}{\sqrt{2^m}} \sum_{k_0=0,1} \cdots \sum_{k_{m-1}=0,1} \ket{k_{m-1} \cdots k_0} exp(\frac{2\pi is}{r} \sum_{l=0}^{m-1} k_l 2^l) \ket{u_s} \\
&= \frac{1}{\sqrt{r}} \sum_{s=0}^{r-1} \frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^m-1} exp(\frac{2\pi is}{r} \cdot k) \ket{k} \ket{u_s} \\
&= \frac{1}{\sqrt{r}} \sum_{s=0}^{r-1} \frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^m-1} exp(\frac{2\pi ik (\tilde{s}/r)}{2^m}) \ket{k} \ket{u_s}  \tag{8}
\end{align}

となります(ここで、$\tilde{s}=2^m s$とおきました)。これに対して逆量子フーリエ変換をかけると、

\frac{1}{r} \sum_{s=0}^{r-1} \ket{\tilde{s}/r} \ket{u_s}  \tag{9}

となり、この状態を測定することで、$\tilde{s}/r$のうちのどれかの値を取得することができます。

測定と位相推定

はじめにちゃんとした固有状態を用意しておく通常の位相推定とは違い、今回、固有状態の重ね合わせからスタートしたため、逆量子フーリエ変換しても残念ながら重ね合わせ状態のままです。ですが、$\ket{\tilde{s}/r}$のうちのどれかが観測されることになるので、ここから$r$を推定することはできます。観測されるのは何らかのビット列になるので、これを小数表現に直します。その小数に対して「連分数展開」(いろいろ解説されているネット記事は容易に見つかると思うので、本記事では説明省略します)を行い、近似する分数を計算し、その分母を位数$r$だとするのです。当たっている保証はないですが、何度か測定を繰り返せば、いずれ当たます。当たっているかどうか(つまりそれが本当の位数$r$なのかどうか)は、実際に$a^r \space mod \space N$を計算して$1$かどうかをみれば良いです。

べき剰余との関係

理屈の説明はだいたい以上なのですが、一つ大事なことを言っておかないといけません。注目していただきたいのは、(8)式の最後から2行目です。式(0)に示した$\ket{u_s}$をこの行に代入すると、どうなるかというと、

\frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^m-1} \ket{k} \ket{a^k \space mod \space N} \tag{10}

となります。見てわかるように、これは、前回の記事で説明した「べき剰余」の量子計算そのものです。つまり、どういうことかと言うと、ユニタリ演算子の位相推定のために、制御ユニタリ演算子をズラズラとかけるようなことは一切必要なく、代わりに「べき剰余」をやって逆フーリエ変換をやれば良いということになります。「Shorのアルゴリズム」では「べき剰余」が本質的な役割を果たすというお話を聞いていましたが、最後にようやくつながりました。¹

量子位数発見の処理ステップ

いままで説明してきた流れを量子計算で行うステップを以下に示します。

[STEP.4-1] 第1レジスタを準備

m量子ビットのレジスタを準備して$\ket{0}$に初期化しておいて、各々にアダマールをかけます。

\ket{0} \rightarrow \frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^{m}-1} \ket{k}

mをどの程度用意するかで位相推定(つまり位数推定)の精度は変わってくるのですが、因数分解したい$N$のビット数の2倍よりも大きくするのが良いとされているようです。

[STEP.4-2] 第2レジスタを準備

$N$のビット数と同じ量子ビット数のレジスタを準備して$\ket{1}$に初期化します($\ket{0}$としておいて、0番目の量子ビットにXゲートをかけて反転させれば良いです)。

[STEP.4-3] べき剰余の計算

第1レジスタと第2レジスタのテンソル積状態を対象にべき剰余計算を実行します。

\frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^{m}-1} \ket{k} \otimes \ket{1} \rightarrow \frac{1}{2^{m-1}} \sum_{k=0}^{2^{m}-1} \ket{k} \otimes \ket{a^k \space mod \space N}

[STEP.4-4] 第1レジスタを逆量子フーリエ変換

\frac{1}{\sqrt{r}} \sum_{s=0}^{r-1} \ket{\tilde{s}/r} \ket{u_s}

が得られます($\tilde{s}=2^m s$)。

[STEP.4-5] 第1レジスタを測定して位数推定

$\ket{\tilde{s}/r}$のうちのどれかが観測されるので、ビット列を小数に直して連分数展開で分母$r$を得ます。本当の位数かどうかは$a^r \space mod \space N$を計算してチェックします。1であればそのときの$r$を返します。何度測定してもダメなら諦めて[STEP.2]からやり直します。

おわりに

実は、数論にはあまり馴染みがなかったので「そもそも位数って何だっけ？」というあたりから勉強しないといけませんでした。幸いネットには良質な解説記事が多数ありました。特に、以下の資料はわかりやすかったです(ありがとうございました)。ご参考まで。

• 整数論と代数の初歩

「はじめに」で述べたように、次回は実際にシミュレータで動作確認をしてみます。「べき剰余」の部分はやはりシミュレータでは量子ビット数が足りないので工夫しました。量子状態を(人為的にですが)作って、量子計算シミュレートしてみました。という意味で「なんちゃってべき剰余」です。詳細は次回ということで。

以上

1. と言ってみましたが、今回対象にしているユニタリ演算子を制御ユニタリ化したものは、制御-剰余乗算演算子そのものです。なので、位相推定をしたとしても、代わりにべき剰余で計算したとしても(前回説明したようなべき剰余回路とすると)、量子回路的には、だいたい同じになります(厳密に同等かどうかはちゃんと見てないですが)。ただ、べき剰余はいろいろ研究されていて効率的な量子回路の組み方があるので、Shorのアルゴリズムはべき剰余を使ってやるんだと思っておけばいんじゃないかな、と思います。 ↩