量子コンピュータの応用 - HHLアルゴリズムについて論文を読んでみる

Tetsuro Tabata

2021/06/02 10:12

#HHLアルゴリズム

§ この記事の目的

量子コンピューティング理論の一つであるHHLアルゴリズムについて、その内容を確認します。
本記事では、HHLアルゴリズムとその背景となる量子コンピューティング理論を体系的に説明している以下の論文をもとに説明しています。

■HHLアルゴリズム
Quantum linear systems algorithms: a primer
https://arxiv.org/pdf/1802.08227.pdf

§ イントロダクション

HHLアルゴリズムとは、ホロー(Aram W.Harrow)氏、ハシディム(Avinatan Hassidim)氏、ロイド(Seth Lloyd)氏らによって2009年に発表された理論で、量子位相推定を使用して効率的に逆行列計算を行う手法に関する理論です。

■Quantum Algorithm for Linear Systems of Equations
(Aram W. Harrow, Avinatan Hassidim, and Seth Lloyd)
https://arxiv.org/pdf/0811.3171.pdf

この理論では、下図[図1]の(a)(b)(c)ように、大きく分けて3つの処理からなります。

(a)は位相推定を行う処理、(b)は制御回転ゲートを作用させる処理、そして(c)は逆演算により解を求める処理です。

順を追って詳しく見ていきます。

§ 位相推定の一般論

まずはじめに、位相推定の理論及び数式を確認します。
位相推定は文字のごとく、位相の値(角度)を求めるアルゴリズムです。あるユニタリ行列によるユニタリ変換を行ったときに、この操作による位相の変化がどのくらいになるかを近似的に求めます。また、与えられた固有状態から精度の高い固有値の近似値を求めることでもあります。

例を示します。
あるユニタリ行列を $U$ 、固有ベクトルを $|\psi〉$ 、固有値を $\lambda$ とすると、以下が成立します。

U|\psi〉=\lambda|\psi〉\quad(式1)

$\lambda$ は実数または複素数で、 $\lambda=e^{i2\pi\phi}$ とすると、

U|\psi〉=e^{i2\pi\phi}|\psi〉\quad(式2)

$2\pi\phi$ は位相の角度を表します。角度の範囲を $0\le2\pi\phi\le2\pi$ とすると $\phi$ は、

0\le\phi\le1\quad(式3)

となります。
古典コンピュータにも共通して言えますが、量子コンピュータも2進数であるビット(量子ビット)を単位として計算するため、求める位相(角度)は近似値となります。

$\phi$ の値が下限値 $\phi=0$ 及び上限値 $\phi=1$ を取るとき、 $\lambda$ も1を取ります。
$\phi$ の値が $0\lt\phi\lt1$ を取るとき、 $\phi$ を2進小数で表現すると、

\phi=0.\tilde{\phi_1}\tilde{\phi_2}...\tilde{\phi_n}

となります。
ここで記号 $\tilde{}$ (チルダ)は2進表現(0または1)であることを意味しています。
よって、以下のように書き表すことができます。

\begin{align} \phi&=0.\tilde{\phi_1}\tilde{\phi_2}...\tilde{\phi_n}\\ &=\frac{\tilde{\phi_1}}{2^1}+\frac{\tilde{\phi_2}}{2^2}+\frac{\tilde{\phi_3}}{2^3}+...\frac{\tilde{\phi_n}}{2^n}\\ &=\sum_{k=1}^{n}\tilde{\phi_k}2^{-k}\quad(式4) \end{align}

位相推定の処理手順を回路で表現すると下図[図2]のようになります。

左上から順に第1量子ビット、第2量子ビット、…と番号(インデックス)を付けます。最下部ビットの $\psi$ を除くすべてのビットは $|0〉$ に初期化されます。

図2の各量子ビットについて、上からインデックスをkとして各量子状態を $|u_k〉$ とします。

$|u_k〉$ の内容を確認していきます。ここでは簡単のために、第1量子ビットに着目して確認します。

まず、量子ビットについてHマークのアダマール変換を適用します。
初期状態は $|0〉$ なので、これにアダマール変換をかけると、

|u_1〉=\frac{1}{\sqrt{2}}(|0〉+|1〉)\quad(式5)

次に、この量子ビットの結果(式5)を入力(コントロールビット)として、 $|\psi〉$ に制御ユニタリゲート(Ctrl-U)を適用します。制御ユニタリゲート(Ctrl-U)は以下のように定義されます。

C_{trl}U=|0〉〈0|⊗I + |1〉〈1|⊗U

よって量子ビットの状態は、

\begin{align} |u_1〉&=C_{trl}U\frac{1}{\sqrt{2}}(|0〉+|1〉)|\psi〉\\ &=(|0〉〈0|⊗I + |1〉〈1|⊗U)\frac{1}{\sqrt{2}}(|0〉+|1〉)|\psi〉\\ &=\frac{1}{\sqrt{2}}(|0〉+U|1〉)|\psi 〉\quad(式6) \end{align}

ここで(式2)より、Uは

U_k=e^{i2\pi\phi_k}

これを(式6)に代入すると、

\begin{align} |u_1〉&=\frac{1}{\sqrt{2}}(|0〉+e^{i2\pi\phi_k}|1〉)|\psi 〉\quad(式7) \end{align}

となります。
(式7)では第1量子ビットのみを確認しましたので、これを[図2]に従い、以下のようにnビット分を計算します。

|u_1〉⊗|u_2〉⊗|u_3〉...|u_n〉|\psi〉

試しにn=3として3ビット分の計算をしてみます。

\begin{align} |u_1〉⊗|u_2〉⊗|u_3〉|\psi〉&=\frac{1}{\sqrt{2}}(|0〉+e^{i2\pi\phi_1}|1〉) ⊗\frac{1}{\sqrt{2}}(|0〉+e^{i2\pi\phi_2}|1〉)⊗\frac{1}{\sqrt{2}}(|0〉+e^{i2\pi\phi_3}|1〉)|\psi 〉\\ &=\frac{1}{\sqrt{2^3}}\bigg( |00〉 + e^{i2\pi\phi_2}|01〉+ e^{i2\pi\phi_1}|10〉+ e^{i2\pi(\phi_1 + \phi_2)}|11〉\bigg)⊗\frac{1}{\sqrt{2}}(|0〉+e^{i2\pi\phi_3}|1〉)|\psi 〉\\ &=\frac{1}{\sqrt{2^3}}\bigg( |000〉 + e^{i2\pi\phi_3}|001〉+ e^{i2\pi\phi_2}|010〉+ e^{i2\pi(\phi_2 + \phi_3)}|011〉+ e^{i2\pi\phi_1}|100〉+ e^{i2\pi(\phi_1 + \phi_3)}|101〉+ e^{i2\pi(\phi_1 + \phi_2)}|110〉+ e^{i2\pi(\phi_1 + \phi_2 + \phi_3)}|111〉\bigg)|\psi〉\\ &=\frac{1}{\sqrt{2^3}}\bigg( e^{i2\pi t_1}|000〉+e^{i2\pi t_2}|001〉+e^{i2\pi t_3}|010〉+e^{i2\pi t_4}|011〉+e^{i2\pi t_5}|100〉+e^{i2\pi t_6}|101〉+e^{i2\pi t_7}|110〉+e^{i2\pi t_8}|111〉\bigg)|\psi 〉\quad(式8) \end{align}

ここで、 $t_k$ を以下のように置き換えました。

t_1 = 0,\; t_2=\phi_3,\; t_3=\phi_2,\; t_4=\phi_2 + \phi_3,\\ t_5 = \phi_1,\; t_6=\phi_1 + \phi_3,\; t_7=\phi_1 + \phi_2,\; t_8=\phi_1 + \phi_2 + \phi_3

ここで、第1量子ビットから第n量子ビットについて、慣例に従い、インデックスを0から数えることにします。
(式8)を用いて、nビットの量子状態を一般的に書くと、

\begin{align} |u〉|\psi〉&=|u_0〉⊗|u_1〉...|u_{n-1}〉|\psi〉\\ &=\frac{1}{\sqrt{2^n}} \bigg(\sum_{k=0}^{2^n-1}e^{i2\pi t_k}|\tilde{k}〉\bigg)|\psi〉\quad(式9) \end{align}

ただし、 $\tilde{k}$ は2進数表示(例： $|0010110〉$ )を意味しています。
以降では $k$ を10進数で書き表し、 $|k〉$ とします。

(式9)では一般的な記載をしました。
よって下図[図3]のようにまとめることができます。

[図3]の通り、制御ユニタリゲート(Ctrl-U)を適用した次の処理として、 $|u〉$ に量子フーリエ逆変換を掛けます。
詳細は省略しますが、量子フーリエ逆変換は以下のように定義されます。

|k〉=\frac{1}{\sqrt{2^n}}\sum_{j=0}^{2^n-1}exp \bigg(-i\frac{2\pi jk}{2^n} \bigg)|j〉

(式9)に代入すると、

\begin{align} |u〉|\psi〉&=\frac{1}{\sqrt{2^n}}\sum_{k=0}^{2^n-1}\bigg( e^{i2\pi t_k}\frac{1}{\sqrt{2^n}}\sum_{j=0}^{2^n-1}e^{(-i\frac{2\pi jk}{2^n})}|j〉 \bigg)|\psi〉\\ &=\frac{1}{2^n}\sum_{k=0}^{2^n-1}\sum_{j=0}^{2^n-1}exp\bigg( i2\pi(t_k-\frac{jk}{2^n})\bigg)|k〉|\psi〉\quad(式10) \end{align}

以上が位相推定の一般的な理論となります。

§ HHLアルゴリズムの基本概念

HHLアルゴリズムは、以下のような線形代数問題で、行列Aとベクトルx,bの関係式が成り立つとき、ベクトルxを求めるためのアルゴリズムです。

Ax=b

解析的に解く場合、Aの逆行列が存在すると仮定すると、以下のように解くことになります。

x=A^{-1}b\quad(式11)

以降で、この問題を量子コンピューティングで解くHHLアルゴリズムについて確認します。

§ (a)位相推定処理

まずは、位相推定の理論を応用した処理である[図1]の(a)位相推定の処理について確認します。
あるエルミート行列である行列Aについて、その固有ベクトルを $|u〉$ 、固有値を $\lambda$ とします。

A|u〉=\lambda|u〉\quad(式12)

線形代数の公理として、 $A$ の逆行列 $A^{-1}$ の固有値は以下となります。

\begin{align} A^{-1}|\lambda〉&=\lambda^{-1}|u〉\\ &=\frac{1}{\lambda}|u〉\quad(式13) \end{align}

また、[図1]のregister Iとして入力される量子ビット $|b〉$ は(式11)の右辺のベクトルbと同じものとします。つまり、(式10)において $|\psi〉$ が $|b〉$ となります。

さらに(式10)において、以下のように置き換えを行います。

|\lambda_i'〉=\frac{1}{2^n}\sum_{k=0}^{2^n-1}\sum_{j=0}^{2^n-1}exp\bigg( i2\pi(\lambda_i-\frac{jk}{2^n})\bigg)|k〉\quad(式14)

ここで $t_k$ を $\lambda_i$ に置き換えました。これは、(式2)における固有値の表現を位相角度で表す代わりに固有値 $\lambda$ で表したことが理由です。

また、 $|b〉$ を固有ベクトル $|u_i〉$ で基底変換します。 $\beta_i$ は基底変換のための変換用定数です。

|b〉=\sum_{i=0}^{n-1}\beta_i|u_i〉\quad(式15)

(式14)(式15)より(式10)は

(式10)=\sum_{i=0}^{n-1}\beta_i|\lambda_i'〉|u_i〉\quad(式16)

となります。

以上が(a)位相推定処理です。
次に、(b)制御回転ゲートを作用させる処理に入ります。

§ (b)制御回転ゲート処理

[図1]に従い、(b)制御回転ゲートを作用させる処理を行います。
ここで、レジスタS(register S：メモリ上の位置)にAncillaビット $|0〉$ を追加します。 Ancillaビットは、アンシラビットまたは補助ビットと呼ばれます。

アンシラビット(補助ビット)を追加すると、量子状態(式16)は以下となります。

(式16) = \sum_{i=0}^{n-1}\beta_i|\lambda_i'〉|u_i〉|0〉\quad(式17)

$|\lambda_i'〉$ をコントロールビットとし、アンシラビットに制御回転ユニタリ演算を行います。具体的には、パウリY行列を適用します。

まずは、制御回転ゲートについて定義を確認します。

\begin{align} C_{trl}U&=\sum_{\tilde{\theta}=0,1}^{}|\tilde{\theta}〉〈\tilde{\theta}|⊗exp(-i\theta\sigma_y)\\ &=|0〉〈0|⊗I+|1〉〈1|⊗exp(-i\theta\sigma_y)\quad(式18) \end{align}

さらに、 $exp(-i\theta\sigma_y)$ に関しては、以下となります。

exp(-i\theta\sigma_y) = \begin{pmatrix} \cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{pmatrix}\quad(式19)

(式18)を(式17)に適用し、(式17)の $|\lambda_i'〉$ の値が $|1〉$ であるときのみを取り出すと、

\begin{align} (式17) &= \sum_{i=0}^{n-1}\beta_i|u_i〉|\lambda_i'〉|0〉 ( |\lambda_i'〉〈\lambda_i'| ⊗e^{-i\theta \sigma_y} )\\ &=\sum_{i=0}^{n-1}\beta_i|u_i〉|\lambda_i'〉e^{-i\theta \sigma_y}|0〉\\ &=\sum_{i=0}^{n-1}\beta_i|u_i〉|\lambda_i'〉\begin{pmatrix} \cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{pmatrix}\begin{pmatrix} 1 \\ 0 \end{pmatrix}\\ &=\sum_{i=0}^{n-1}\beta_i|u_i〉|\lambda_i'〉(\cos \theta|0〉+\sin \theta|1〉)\quad(式20) \end{align}

上記計算途中で(式19)も適用しました。

ここで、 $\sin \theta = \frac{C}{\lambda_i'}$ とおきます。 Cは規格化定数です。

(式20)は以下となります。

\begin{align} (式20) &= \sum_{i=0}^{n-1}\beta_i|u_i〉|\lambda_i'〉\bigg(\sqrt{1-\frac{c^2}{\lambda_i'^2}}|0〉+\frac{C}{\lambda_i'}|1〉\bigg)\quad(式21) \end{align}

以上で(b)制御回転ゲート処理を終わります。
最後に、(c)逆演算により解を求める処理を行います。

§ (c)逆演算により解を求める

最後に、逆演算(uncomputation)と呼ばれる処理を行います。
[図1]の(c)の通り、レジスタCである $|\lambda_i'〉$ に対して $|0〉$ に戻す処理を行います。不要な量子ビットは $|0〉$ に戻すことが目的です。

(式21)は以下となります。

\begin{align} (式21) &= |0〉\sum_{i=0}^{n-1}\beta_i|u_i〉\bigg(\sqrt{1-\frac{c^2}{\lambda_i'^2}}|0〉+\frac{C}{\lambda_i'}|1〉\bigg)\quad(式22) \end{align}

上記で(c)逆演算処理が終わりました。
この後は、アンシラビットの測定を複数回行います。そしてアンシラビットの値が $|1〉$ となる場合のみを抽出します。その結果、(式22)から以下が得られることがわかります。

\sum_{i=0}^{n-1}\beta_i|u_i〉\frac{C}{\lambda_i'}

上記は[図1]の右下の出力結果の $|x〉$ に相当することになります。
書き直すと、

|x〉=C\sum_{i=0}^{n-1}\frac{\beta_i}{\lambda_i'}|u_i〉

さらに(式15)により変換すると、

|x〉=C\frac{1}{\lambda_i'}|b〉\quad(式23)

(式23)は(式11)と同じ形となっています。つまり、(式23)を算出することで最終目的の $|x〉$ を求めることができます。

以上がHHLアルゴリズムの内容となります。

§ ご意見など

ご意見、間違い訂正などございましたらお寄せ下さい。