量子位相推定による固有値・固有ベクトル問題

はじめに

これまで研究での活用がメインであった量子情報や量子コンピュータがいよいよ実用化に向けて民間利用が進んできました。その中で量子コンピュータの機能を活用してできることを模索するのはとても大事です。

量子コンピュータが期待される理由として、

１、既存コンピュータが性能限界にきている
２、既存コンピュータの問題点を解決するためにサイズを小さくすると量子力学を考慮する必要があり、量子力学に基づいて計算原理を組み直すとデータの使い方が変わる。
３、量子力学の原理を使うとこれまでのコンピュータと異なることができる。

などです。今回はその中でも特に特徴的な波動関数や状態ベクトルと呼ばれる既存計算機では持っていない仕組みを簡単に理解して、計算機科学としてどのように活用すれば良いかを確認してみたいと思います。やってみて思うのは波動関数とか、状態ベクトルとか名前が物理っぽくなっていますが、実際には物理記号を使わなくてもほぼ数学の行列演算だけに落とし込むことができるので、先入観をなくせば理解しやすいと思います。

基本的演算

量子コンピュータのかなめとなる要素は、「状態ベクトル」と「ゲート演算」になります。この２つを把握して早速使いこなしに進みましょう。

状態ベクトル

量子コンピュータではデータの持ち方がこれまでのコンピュータと異なり、「状態ベクトル」と呼ばれるデータの01の組み合わせの配列を持ちます。これまでのコンピュータは01110101011001...010101のように、操作した結果としてのッビット配列を持ちますが、量子コンピュータではこれを拡張し、すべてのビットの組み合わせの配列を情報として持ちます。それを状態ベクトル（まはた波動関数）と呼びます。

全てのビットの組み合わせ配列は二進数で表され、それぞれの組みあわせは確率振幅と呼ばれるそれぞれの状態の出現確率の平方根で情報が格納されます。この状態ベクトルを操作し、把握することで計算を行います。状態ベクトルはビットの組み合わせに対応するので、スペックとしてある量子ビット数 $N$ に対して、 $2^N$ の組み合わせを状態ベクトルとして持ちます。これはNが増えるごとに扱えるデータ量が指数増大することを意味し、この機能をうまく使うことで既存計算機では原理的に計算できないことを実行できるのではと期待されているのが量子コンピュータです。

その中でももっとも有名なのが素因数分解を実行するshorのアルゴリズムで、現在のスパコンで何万年もかかるような計算を数時間で解いてしまうことができると言われています。現在はまだ量子コンピュータのサイズが小さいため実質的に解くことができません。

状態ベクトルは、

\mid \psi \rangle = \begin{bmatrix} 0.01\\ 0.2\\ 0.003\\ \vdots\\ 0.1 \end{bmatrix}

のように縦ベクトルで表現され、上から順番にビットの配列の組み合わせの出現確率の平方根が格納されます。例えば上記が6量子ビットの場合、

\mid \psi \rangle = \begin{bmatrix} \mid 000000 \rangle\\ \mid 000001 \rangle\\ \mid 000010 \rangle\\ \vdots\\ \mid 111111 \rangle \end{bmatrix}

のように、

上から順番に計算結果が最終的に000000になる確率が $0.01^2=0.0001$ とか、000001になる確率が $0.2^2=0.04$ のように読み取ることができます。

最終的な計算結果は1回の演算で1つしか取り出せませんので、なんども計算をする必要が出ることがあります。

量子ゲート

上記の状態ベクトルを操作するために必要なのが量子ゲートです。量子ゲートは、

X = \begin{bmatrix} 0&1\\ 1&0\\ \end{bmatrix}

のように行列演算の形で表現され、１量子ビットの演算と２量子ビットの演算との２種類があります。基本的に１量子ビットを操作するだけで、複数の状態ベクトルの要素に影響を与えるところが特徴です。量子ゲート演算を使いことなすことで、ビットそのものを操作するだけでなく、ビットの組み合わせそのものを操作できます。これらを使ってできることを確認したいと思います。

基本的には行列の演算ですので、状態ベクトルにゲート演算をかけることによって状態ベクトルを操作することができます。

\mid \psi \rangle = X\mid 0 \rangle = \begin{bmatrix} 0&1 \\ 1&0 \end{bmatrix} \begin{bmatrix} 1 \\ 0 \end{bmatrix} = \begin{bmatrix} 0 \\ 1 \end{bmatrix}

量子位相推定による固有値・固有ベクトル問題

量子コンピュータは現在のコンピュータのデジタル01計算の上位互換と呼ばれています。特徴的な機能として、現在のコンピュータでは性能的に難しいと言われているサイズの行列の固有値・固有ベクトルの問題を期待値を元に解くというのが期待されています。

具体的には $N$ 量子ビットの量子コンピュータで $2^N$ の状態ベクトルを持ち、それを操作するのに、 $2^N * 2^N$ の巨大な量子ゲート演算を行います。現在のスパコンでもメモリの性能限界で $2^{48} * 2^{48}$ の行列サイズの計算が限界と言われていますので、それ以上のサイズの量子コンピュータによって、今まででは原理的に計算のできない計算ができることが期待されています。

今回、量子ゲートで行列を構成し、その固有ベクトルを状態ベクトルで表現すると、固有値を求めることができるアルゴリズムの量子位相推定に関してみてみたいと思います。

行列の固有値を求める

wikipediaによると、

有限次元線形空間V上の線形変換Aに対して、次の方程式

Ax = \lambda x

を満たす零でないベクトル x とスカラー λ が存在するとき、x を A の固有ベクトル（右固有ベクトル）、λ を A の固有値と呼ぶ。

https://ja.wikipedia.org/wiki/固有値

ということです。将来的にAのサイズが大きくなっても効率的に固有値がもとまると様々な問題が解けるようになります。

行列の固有値を状態ベクトルとゲート演算で求める

行列の固有値は、線形変換Aとベクトルxでスカラーλを求めるのでした。

Ax = \lambda x

ですので、Aに「ゲート演算」、xに「状態ベクトル」を利用することでAの固有値であるλを求める手法があれば良いことになります。状態ベクトルは $\mid \psi \rangle$ の形で書かれますので、

H \mid \psi \rangle = \lambda \mid \psi \rangle

とできればOKです。Hはハミルトニアンと呼ばれる演算子で、ゲートの塊＝行列と考えれば大丈夫です。

Hの行列は求めたい行列で固定になりますので、あとは $\mid \psi \rangle$ が作れれば良くなります。

量子状態の時間発展

シュレディンガー方程式の式変形から下記の時間発展が導き出されます。

\mid \psi (t)\rangle = exp(-i/ \hbar H(t-t_0))\mid \psi(t_0)\rangle

こちらは時間t=0の状態にあるハミルトニアンを作用されると量子状態を変化させることができるというものです。これを断続的に行い状態ベクトルの $\mid \psi \rangle$ を作っていきます。もっと簡単に書くと、

\mid \psi\rangle = e^{-iHt}\mid \psi_0\rangle

となります。初期状態から少しずつ状態を変化させて $e^{-iHt}$ に対応する操作を量子ゲートで実現することで状態ベクトル $\mid \psi \rangle$ を作り出すことができます。

行列指数関数と量子ゲートによる対角化、トロッタ展開による近似

時間発展は、指数の計算を量子ゲートのユニタリ演算に落とし込むことができます。基本的には対角行列のZゲートを使って、それ以外はZゲートを使って対角化を行います。それらを組み合わせて計算をしますが、指数の肩のハミルトニアンHがユニタリでない場合には、トロッター展開を使ってステップ数nを利用して近似計算で処理をします。

e^{-iHt} = e^{-i(A+B)t} = (e^{-iAt/n}e^{-iBt/n})^n

上記のような展開方法を使って時間発展のシミュレーションをします。時間発展のパラメータとしてトロッタ展開のnと時間発展の時間tがありますが、状態ベクトルは最終的に固有ベクトルとして固有値の期待値を求める際に重要になりますが、簡単にもとまるわけではなく調整に苦労しそうです。

量子位相推定アルゴリズム（Quantum Phase Estimation）

（注意）最初に注意しておくのは、この位相推定アルゴリズムは現在の量子コンピュータではエラーが多すぎて実際の実行が難しいことがあります。VQEというアルゴリズムの方が現在のエラーの多いマシンでは実用的なのでその点を理解した上で読んでいただくことをお勧めします。

量子位相推定アルゴリズムは固有値 $\lambda$ を角度の形で求める方法です。上記の方法で量子状態 $\mid \psi \rangle$ が準備できたら、２段階の方法で固有値を求めます。

\lambda = e^{-i\alpha}

固有値 $\lambda$ をやはり波動の形で書いて、 $\alpha$ を求めます。状態ベクトルは直接は人間は見ることができません。まずは状態ベクトル $\mid \psi \rangle$ を位相を取り出せる形に変形します。

位相推定アルゴリズムの全体を見る前にこの量子位相推定アルゴリズムの要である「量子フーリエ変換」を先に見ます。量子フーリエ変換は入力を01ビットの配列から、その配列に対応する位相（角度）をもつ量子状態の形に変形をすることができます。

{QFT:\mid x \rangle \mapsto \frac{1}{\sqrt{N}}\sum_{k=0}^{N-1} \omega_n^{xk}\mid k\rangle }

$\omega_n = e^{\frac{2\pi i}{N}}$ とすると、

{F_N= \frac{1}{\sqrt{N}} \left[ \begin{array}{rrrr} 1 & 1 & 1 & \cdots &1\\ 1 & \omega_n&\omega_n^2&\cdots&\omega_n^{N-1}\\ 1 & \omega_n^2&\omega_n^4&\cdots&\omega_n^{2(N-1)}\\ 1 & \omega_n^3&\omega_n^6&\cdots&\omega_n^{3(N-1)}\\ \vdots&\vdots&\vdots&&\vdots\\ 1 & \omega_n^{N-1}&\omega_n^{2(N-1)}&\cdots&\omega_n^{(N-1)(N-1)} \end{array} \right] }

のようになります。既存の計算機の高速フーリエ変換とほぼ同じなので、詳しく勉強したい方は高速フーリエ変換から学ぶと良いと思います。

量子フーリエ変換は数学的な処理なので、みて覚えてしまっていいと思いますが、 $x_1$ から $x_n$ のビットを入力として入れると、それに対応した位相が量子状態の形で出力されます。

{QFT(\mid x_1,x_2,…,x_n \rangle) = \frac{1}{\sqrt{N}}(\mid 0 \rangle + e^{2\pi i [0.x_n]} \mid 1 \rangle) \otimes … \otimes (\mid 0 \rangle + e^{2\pi i [0.x_1x_2…x_n]} \mid 1 \rangle) }

出力された量子状態はそれぞれの量子ビットは位相の角度の精度に対応された形に出力されますが、それぞれの位相に対応する係数は量子ビットから見ると絶対値はすべて1なので、測定しただけでは0と1が完全に50%ずつ出てしまうのが特徴です。

{[0.x_1x_2…] = \frac{x_1}{2}+\frac{x_2}{2^2}+… }

位相は小数点の形で取り出されます。

量子フーリエ変換はフーリエ変換できますが、結果を測定することができません。今回はこの量子フーリエ変換を逆から利用し、逆量子フーリエ変換を使って、量子状態から逆に位相のビット表記を取り出します。

量子コンピュータは可逆回路といって回路を逆向きに実行すると計算結果を戻すことができます。量子フーリエ変換の逆である逆量子フーリエ変換は、フーリエ変換の回路を逆向きに実行するだけで実現できます。

つまり上記量子位相推定では、作った状態ベクトル $\mid \psi \rangle$ の位相である $\lambda = e^{-i\alpha}$ を求めたい場合、状態ベクトルを上記の量子フーリエ変換の形に対応させた形で変形させ、それの状態ベクトルを再度逆量子フーリエ変換にかけることで、任意の制度で位相を小数点表記で取り出すことができます。

量子位相推定アルゴリズムの実装回路をみてみる

量子位相推定アルゴリズムの全体像は下記の通りです。

まず①の重ね合わせですが、Hゲートを最終的に位相を取り出すビットにかけます。仕組みはあとで確認します。次に②の量子状態の作成を行います。こちらは求めたいHに対応する量子状態を作成します。

（すみません、①と③をセットにした方がよかったです。。。今後そうします。）

まず量子状態を作成するには時間発展シミュレーションをすればいいのでした。これ上記でざっくりと概要が書いてありますので、それにそってあとでもう少し詳しくみてみたいと思います。

位相情報を移す

量子状態が時間発展で準備できて、最終的に位相をビット配列の形で逆量子フーリエ変換で取り出す下準備ができたところで、量子状態から位相を取り出す操作である位相推定のアルゴリズムを見てみたいと思います。

こちらが大変参考になりました。
「量子位相推定アルゴリズム (Quantum Phase Estimation) とは」

上記の量子フーリエ変換に関して、任意の精度でのビットの出力に関して式をまとめると、

\frac{1}{\sqrt{2^n}}\sum_{k=0}^{2^n-1}e^{i 2\pi k \frac{j}{2^n}}\mid k \rangle

となります。こちらを量子フーリエ変換すればいいのでした。まずは最初に準備した $\mid 0 \rangle$ にアダマールゲートを準備することで、

\frac{1}{\sqrt{2^n}}\sum_{k=0}^{2^n-1}\mid k \rangle

が準備できます。量子フーリエ変換する前と比較すると、係数 $e^{i 2\pi k \frac{j}{2^n}}$ が異なっており、これを実現できれば回路は実現できそうです。ここで、 $\frac{j}{2^n}$ は定数なので、 $\phi$ と書き直すと、

\frac{1}{\sqrt{2^n}}\sum_{k=0}^{2^n-1} e^{i2\pi k\phi}\mid k \rangle

となりますが、 $e^{i\theta}$ はゲートの回転に相当し、係数kによって対応するビットだけ回転をすればいいということになります。不思議ですね。。。

位相キックバック

基本的には固有値固有ベクトル問題を利用し、ハミルトニアンHを導入したユニタリ回路を制御付きにすることにより、固有ベクトル $\mid \psi \rangle$ に対応する固有値 $e^{2\pi i \phi}$ をコントロールビットの方の $\mid 1 \rangle$ 状態にかけることができます。

対応するハミルトニアンに対する固有値固有ベクトルは、

U\mid \psi \rangle = e^{2\pi i \phi} \mid \psi \rangle

で書くことができます。量子状態 $\mid \psi \rangle$ は準備されてますので、ハミルトニアンを導入したユニタリ回路Uを準備することで、位相を移すことができます。

まず初期状態からアダマールゲートを使って結果を取り出したい方の量子ビットを重ね合わせ状態にします。

\mid 0 \rangle \mid \psi \rangle \rightarrow \frac{\mid 0\rangle + \mid 1 \rangle}{\sqrt{2}} \mid \psi \rangle = \frac{\mid 0\rangle \mid \psi \rangle + \mid 1 \rangle \mid \psi \rangle}{\sqrt{2}}

次に制御付きのユニタリ回路（ハミルトニアンに対応した形）を導入することで、制御付きなので量子ビットが $\mid 1 \rangle$ の方にだけユニタリがかかります。

\frac{\mid 0\rangle \mid \psi \rangle + \mid 1 \rangle U \mid \psi \rangle}{\sqrt{2}}

そして、これは固有値に対応しているので、

\frac{\mid 0\rangle \mid \psi \rangle + \mid 1 \rangle e^{2\pi i \phi} \mid \psi \rangle}{\sqrt{2}} = \frac{\mid 0\rangle \mid \psi \rangle + e^{2\pi i \phi} \mid 1 \rangle \mid \psi \rangle}{\sqrt{2}} = \frac{\mid 0\rangle + e^{2\pi i \phi} \mid 1 \rangle}{\sqrt{2}} \mid \psi \rangle

このようになりました。これによって、固有値をかけることができます。あとは、

\frac{1}{\sqrt{2^n}}\sum_{k=0}^{2^n-1} e^{i2\pi k\phi}\mid k \rangle

のように、取り出したい桁に対応するkを導入すれば良いですが、これは回転角に対応しており、固有値をk回かけるということをすると、自然と対応することができます。つまり $U^k$ のように同じユニタリ操作をk回実行すれば良いことになります。ということで、k回Controlled-Unitary操作を行うことで求めることができます。

\frac{\mid 0\rangle + U^k \mid 1 \rangle}{\sqrt{2}} \mid \psi \rangle = \frac{\mid 0\rangle + e^{2\pi i k \phi} \mid 1 \rangle}{\sqrt{2}} \mid \psi \rangle

これにより、対応する桁kに対してk回の制御付きユニタリゲートをかけることで実行が終了します。

素材は揃ったのであとは実行

これですべての回路が揃いました。

①位相キックバックのための重ね合わせの準備
②時間発展による量子状態の準備
③位相キックバックによる固有値固有ベクトルによる位相の反映
④逆量子フーリエ変換により準備された量子状態の位相をビットの配列で取り出す

これによって初期にときたい行列をハミルトニアンとして準備し、時間発展によって量子状態を準備。量子状態が準備できれば、今度はハミルトニアンをユニタリ操作として制御付きで反映し、固有値固有ベクトルで量子状態を作成。逆量子フーリエ変換によって位相情報のビット配列を取り出す。

これによって量子ビットが増えれば将来的に大きな行列の固有値を効率的（？）に求めることができそうです。shorのアルゴリズムでも、この量子位相推定を利用してうまく位相を取り出すことで素因数分解に役立てることができています。

blueqatで実装

量子コンピュータを使った固有値固有ベクトル問題を解くためのアルゴリズムである位相推定アルゴリズムを実装してみます。ときたい行列をハミルトニアンの形にして、量子状態を準備し、それに対する固有値を推定します。早速実装してみましょう。

簡単な位相推定

まずは少ない量子ビットでやってみます。回路は主に、

１、固有ベクトルとしての量子状態の準備（時間発展・今回省略）
２、固有値を量子状態として取り出す（位相キックバック）
３、固有値の量子状態をビットに変換（量子フーリエ変換）

の３つからなります。

Zゲートの位相を推定する

ひたすら例題をやってみます。まずはハミルトニアンとしてZゲートを用意してみます。

Z = \begin{bmatrix} 1&0\\ 0&-1 \end{bmatrix}

まずは手計算で答えを確認してみると、固有値は、

det\begin{bmatrix} 1-\lambda&0\\ 0&-1-\lambda \end{bmatrix}

を計算して、 $\lambda = 1,-1$ で、固有ベクトルは、

\begin{bmatrix} 1\\ 0 \end{bmatrix}, \begin{bmatrix} 0\\ 1 \end{bmatrix}

こうなるはずです。早速確かめてみます。

0、回路の全体概要

今回の回路の全体概要です。

|0> ----H--*--iQFT--M
           |
|0> -------Z--------

まず２量子ビット準備します。0番目、1番目と名前をつけます。どちらも量子ビットは0からスタートします。

１、今回は最初っから量子状態を準備します。
２、次に0番目の量子ビットを重ね合わせ状態にして、CZゲートで制御付きのハミルトニアンを実行して位相を0番目の量子ビットにキックバックします。
３、最後に出来上がった量子状態から位相をビット表記で取り出すための逆量子フーリエ変換を実行し、測定して終わりです。

１、量子状態の準備

固有ベクトルを|0>にして計算しますので、何もしないで大丈夫そうです。

２、位相キックバック

次に量子状態から位相をキックバックします。

|0> --H--*--iQFT--M
         |
|0> -----Z--------

量子状態が準備できました。0番目の量子ビットを重ね合わせにして、量子状態を固有ベクトルとして、固有値を0番目の方に情報としてうつします。

from blueqat import Circuit
Circuit().h[0].cz[0,1].run()
#array([ 0.70710678+0.j,  0.70710678+0.j,  0.        +0.j, -0.        +0.j])

array([ 0.70710678+0.j,  0.70710678+0.j,  0.        +0.j, -0.        +0.j])

となりました。

３、最後量子フーリエ変換

1量子ビットのフーリエ変換はアダマールゲートHと同じですので、Hを適用すると、

Circuit().h[0].cz[0,1].h[0].run()
#array([1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j])

array([1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j])

\phi = 0.0

となりますので、求める固有値は、

e^{2\pi i \phi} = e^{2\pi i *0} = e^0 = 1

となり、固有ベクトルと固有値は、

\mid \psi \rangle = \begin{bmatrix} 1\\ 0 \end{bmatrix}, \lambda = 1

となりました。

続きまして固有状態を|1>に

最初の量子状態を|1>にしてみます。

|0> --H--*--iQFT--M
         |
|0> --X--Z--------

Circuit().x[1].h[0].cz[0,1].h[0].run()
#array([0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j])

array([0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j])

こちらは11となりました。0番目の量子ビットが1となるので、 $\phi = 1/2 = 0.5$ となりました。これを計算して、

e^{2\pi i *0.5} = -1

となりましたので、固有値 $\lambda = -1$ となり対応が取れています。

続きましてXゲート

Xゲートは、

X = \begin{bmatrix} 0&1\\ 1&0 \end{bmatrix}

となります。固有ベクトルとしてまずは、

\mid \psi \rangle = \begin{bmatrix} 1\\ 1 \end{bmatrix}

を考えると、

|0> --H--*--H--M
         |
|0> --H--X-----

を実行してみます。

Circuit(2).h[:].cx[0,1].h[0].m[0].run(shots=100)
#Counter({'00': 100})

Counter({'00': 100})

こちらは、 $\phi = 0$ となり、

\lambda = e^0=1

となりました。続いて固有ベクトルが、

\mid \psi \rangle = \begin{bmatrix} 1\\ -1 \end{bmatrix}

を考えます。

|0> --H---*--H--M
          |
|0> --HZ--X-----

HZゲートで量子状態を準備します。

Circuit(2).h[:].z[1].cx[0,1].h[0].m[0].run(shots=100)
#Counter({'10': 100})

Counter({'10': 100})

これより、 $\phi = 0.5$ となり、

\lambda = e^{2\pi i *0.5}=-1

となりました。

考察

今回の量子位相推定は汎用量子コンピュータの実現を待つ必要があり、数十年かかると言われています。現在使われているVQEと言われる既存計算機とのハイブリッド方式も今後使い方を見てみたいと思います。

量子位相推定による固有値・固有ベクトル問題

Yuichiro Minato

はじめに

基本的演算

状態ベクトル

量子ゲート

量子位相推定による固有値・固有ベクトル問題

行列の固有値を求める

行列の固有値を状態ベクトルとゲート演算で求める

量子状態の時間発展

行列指数関数と量子ゲートによる対角化、トロッタ展開による近似

量子位相推定アルゴリズム（Quantum Phase Estimation）

量子位相推定アルゴリズムの実装回路をみてみる

位相情報を移す

位相キックバック

素材は揃ったのであとは実行

blueqatで実装

簡単な位相推定

Zゲートの位相を推定する

0、回路の全体概要

１、量子状態の準備

２、位相キックバック

３、最後量子フーリエ変換

続きまして固有状態を|1>に

続きましてXゲート

考察