本記事では量子振幅増幅と、それを利用したAbramsらの提唱した方法による量子数値積分法[1]を詳細に説明します。グローバーの検索アルゴリズムから発展した量子振幅推定アルゴリズムは、様々な分野で応用されています。その一つが数値積分です。まずは量子振幅増幅アルゴリズムについて簡単に説明します。

前提知識:

blueqatで使用可能なH,X,回転ゲート,制御NOTゲートなど初歩的な量子ゲートと量子回路、ブラケット記法。

量子振幅増幅アルゴリズム
量子振幅増幅アルゴリズムは、特定の量子状態における存在確率を繰り返しオラクルゲートと増幅ゲートを作用させることで増幅し、それと直交する全状態の中からそれだけを選ぶアルゴリズムです。これは、代表的な量子検索アルゴリズムであるグローバーの検索アルゴリズムを方法の大枠を変えずに拡張したもので、その流れは次のようになります。

1, 作業用量子ビットn個、アンシラ量子ビット1個を用意する。
2, 解となる状態を$\mid \Psi_0 \rangle$、そうでない状態を$\mid \Psi_1 \rangle$とした場合、$A\mid 0 \rangle_{n+1} = sin\theta \mid \Psi_0 \rangle + cos\theta\mid \Psi_1 \rangle$とする操作$A$をn+1量子ビットからなる初期状態$\mid 0 \rangle_{n+1}$に作用させる。

3, 解となる状態のみの位相を逆転させるオラクル演算子$G_{\Psi_0}=A^{\dagger}OA$、$G_A=A^{\dagger}(2\mid 0 \rangle_n \langle 0 \mid_n - I)A$の積$G_{\Psi_0}G_A$を$(\pi - 2\theta)/4\theta$回作用させます。

4, $\mid \Psi_0 \rangle$の基底で観測します。

上の流れを見るに、こちらの記事に示しましたグローバーのアルゴリズムと大枠は同じであることがわかります。量子回路は図1のようになり、こちらもほぼ共通です。すなわち、量子振幅増幅はグローバーのアルゴリズムの一般化でしかないことがわかります。しかしながら、こちらの方は数値積分へ応用が可能です。

図1 量子振幅増幅の量子回路。

計算機上の数値積分

古典計算機は、離散データしか扱うことは出来ません。そのため、計算機上では数値積分も区分求積法で扱われることになります。区分求積法とは、定積分を行う範囲内を各々の積分変数で等間隔に切り分け、図2のように切り分けたそれぞれの面積を近似的に求める方法です。一次元変数xでf(x)を積分する場合、式で表すと、


\begin{align}
F(x) &=& \int f(x)dx \\

& \approx & S(f) \\

S(f)&=&\sum_{x=0}^{2^n-1}p(x)f(x) \\
&=&\sum_{x=0}^{2^n-1}\frac{1}{2^n}f(\frac{x}{2^n-1}) \\
\end{align}

です。また、この式は量子計算機で計算することも可能です。

図2 区分求積法の模式図。

量子数値積分

この節では、いよいよ量子計算機を用いて式(1)を計算します。その方法は3つほど提唱されており、そのうちもっとも古いものは量子位相推定を用いるものです。しかし、これではNISQデバイスにおける実装と忠実度の高い計算を行うのは困難なため、位相推定を用いない方法が考案されました。今回はその1つであるAbramsらのグループが提唱したものを使用して、量子数値積分を行う方法を説明します。この方法は、図3のように量子振幅増幅を行い、その計算結果の誤差範囲を評価し、その誤差範囲で再度量子振幅増幅を行い、その結果を$S(f)$に加算していき、その推定値の精度を高めることで、目的の数値積分を計算する方法です。

図3 Abramsらのグループが提唱した量子数値積分の方法。

その誤差範囲はk回目の計算においては$S_{k-1}(f)-\epsilon_{k-1}<S(f)<S_{k-1}(f)-\epsilon_{k-1}$となります。まず、先ほどの$A$に似通った演算子、


\begin{equation}
B=(P^\dagger \otimes I)R(P\otimes I)
P\mid 0 \rangle = \sum_{x=0}^{2^n - 1}\sqrt{p(x)}\mid x \rangle
R\mid x \rangle \mid 0 \rangle = \mid x \rangle (\sqrt{f(x)}\mid 0 \rangle + \sqrt{1 - f^2(x)}\mid 1 \rangle)
\end{equation}

を定義します。ここでf(x)は式(1)のものと同じです。これをn+1量子ビットの状態$\mid 0 \rangle_{n+1}$に作用させることで、状態は、

 
\begin{equation}
B\mid 0 \rangle_{n+1}=\sum_{x=0}^{2^n-1}p(x)f(x)\mid 0 \rangle + \cdots
\end{equation}

となります。このうち$\mid 0 \rangle$の係数を観測によって求め、その値を$S_1(f)$とします。それを基に$f_2(x)$、$B_2(f)$を、


\begin{equation}
f_k(x)=f(x)-(S_{k-1}(f)-\epsilon_{k-1})
D_k(f)=\sum_{x=0}^{2^n-1}p(x)f_k(x)
\end{equation}

に従って計算し、そこから式(2)と同様に新たな$B_k$演算子を作ります。ここで、$\epsilon_{k-1}$はk回目の計算において許容する誤差です。それを用いて$A=B_k$とした量子振幅増幅アルゴリズムを行います。但し、それにおける反復回数$j_k$は、


\begin{equation}
(2j_k + 1)sin^{-1}\epsilon_k = \pi / 2
\end{equation}

に従って決めます。こうして$k>2$におき得る$S_k(f)$は、


\begin{equation}
S_k(f)=S_{k-1}(f)-\epsilon_k(\frac{1}{2}-2D_k(f))
\end{equation}

となります。ここからは、$S(f)=sin^2\theta$としてこの手法を用いて計算を行う方法とその結果を説明します。まず、$S(f)$の値から、$f(x)=sin2/(2^n-1)x\theta$です。$R$及び$R_k$演算子によって作られる状態において、$f(x)$は単項式で表される状態しか許されないため、$f_k(x)$は、


\begin{equation}
2\delta_k = sin^{-1}(\epsilon_{k-1} - S_{k-1}(f))
f_k(x)=2sin(\frac{1}{2^n-1}x\theta + \delta_k)cos(\frac{1}{2^n-1}x\theta - \delta_k)
\end{equation}

とします。上式は積和の公式を用いて簡単に導出できます。ここで、$\epsilon_k=1/(2(2j_k+1)\sqrt{N_1})$としました。これを図4の回路を用いて実装し、量子振幅増幅を用いて積分値が計算出来ます。$\theta=\pi/2$の場合における4量子ビットを用いて$S(f)$を求めた結果を図5に示します。

図4 $B$を実装するための量子回路。標的ビットが0のときのみsinを掛けるために$\pi/2$位相をずらしています。

図5 $\pi/2$におけるkに対する$S$の推移。$\epsilon_d$、$\epsilon_u$はそれぞれ誤差範囲の下限と上限、exactは区分求積法における厳密解です。

今回、ショット数は25としました。k=2 厳密解は比較的高精度で求められますが、k=3 で厳密解が誤差範囲から逸れてしまう為そのまま誤差範囲が収束します。誤差範囲における初期値はk=1 で厳密解を含むように、また、誤差範囲を拡大した後も厳密解がその範囲に入るように工夫する必要があります。また、量子ビット数が少ないせいで積分範囲の分割が粗いことも関係している可能性があります。

このアルゴリズムは、量子ビットを増やすほどに積分範囲の分割数を$2^n$オーダーで増やすことが可能です。量子ビットで実装する以上、被積分関数は周期性のある三角関数に直す必要がありますが、量子ビット数次第では古典計算機を凌ぐ精度で数値積分が可能になります。今日様々な分野で熱心に研究されている関数の多くは周期関数です。なので、より多くの量子ビットでこの手法が実装された暁には、古典計算機では現実的な時間での計算が不可能な積分が量子計算機によって解かれるでしょう。

[1] 宇野隼平、量子コンピュータを用いた高速数値積分、 https://www.mizuho-ir.co.jp/publication/giho/pdf/010_01.pdf