戻る

量子コンピュータの原理

作成: 2022/09/25

目次

量子力学と量子コンピュータ、量子ビット

量子というのは、粒子と波の両方の性質を両方もつ原子、電子、中性子、陽子、光を粒子とみんなした光子などが量子の代表的なものです。量子力学的な現象を用いているわけですから、量子力学とは切っても切れない関係にあるわけです。ここでは、まず、量子コンピュータで様々な問題を解く際に押さえておくべき基本的な概念をみてみます。

量子の振る舞いは、次のシュレディンガーの波動方程式を解くことによってもとまります。

この式の両辺はエネルギーであって、粒子が持つエネルギーは、波動関数にハミルトニアンという運動エレルギーと位置エネルギーの和からなる演算子$H$$H$を作用させることによって計算できます。波動関数の時間微分がエネルギーに等しいということです。つまり、量子状態が時間的に変化するには、エネルギーが必要で、時間変化の原因となるものがエネルギーだということができます。

波動関数に作用する演算子であるハミルトニアンの固有値と固有関数を求めれば、波動関数が解け、量子力学的な状態がわかります。ハミルニアンという名称は、１８３０年代に活躍したウィリアム・ローワン・ハミルトンにちなんでいます。ハミルトニアンを、波動関数に作用させるとエネルギーが計算できます。ハミルトニアンは行列ですから、エネルギー変換行列とでも呼んだほうが分かりやすいかも知れません。

ハミルトニアンの固有値と固有関数を求め、最低エネルギー状態とそれより一つエレルギーが高い状態が求まったとします。

２つの状態だけを考えて、量子状態は、この２つの状態に対応する波を重ね合わせたものとして記述できます。複素数 $\alpha$$\alpha$と $\beta$$\beta$ を用いて、以下のように書けます。

${\psi }_0(x)$$\psi_{0}(x)$ や ${\psi }_1(x)$$\psi_{1}(x)$ は、定数ですから、それを明示的に書く必要はなく、係数だけで表す ことができます。これが量子ビットです。

ハミルトニアンとエネルギー

量子状態にハミルトニアンをかけるとエネルギーになるということは重要です。たとえば、量子古典ハイブリッドアルゴリズム（NISQアルゴリズム）で解くことの出来る問題の典型は、問題がハミルトニアン、つまり状態をエネルギーに変換する行列が与えられ、安定なエネルギー状態を求めることに置き換えることができます。つまり、固有値と固有ベクトルを求めることになります。または、言い換えると安定なエネルギー状態のうち、一番エネルギーが小さいものを求めるということになります。

量子コンピュータの動作原理

量子コンピュータは、「重ね合わせ」や「量子もつれ」と言った量子力学的な現象を用いて従来のコンピュータでは現実的な時間や規模で解けなかった問題を解くことが期待されるコンピュータです。量子コンピュータの世界では、従来のコンピュータは、量子力学的な現象に基づいていないという意味で、古典コンピュータと呼ばれます。

量子コンピュータは、量子力学的な重ね合わせによって、$n$$n$個の量子ビットを用いて$2^n$$2^n$の状態を同時に処理できるため、古典コンピュータでは、時間がかかりすぎ解くことが著しく困難な問題が解けることが期待されています。

古典ビットと量子ビットの違い

図１にしめすように、古典コンピュータでは、ビットは、０と１かのどちらか値をとります。量子ビットでは、基底ベクトル |0⟩ |1⟩ に係数を掛けて足し合わせた状態ベクトル $|\psi ⟩=\alpha |0⟩+\beta |1⟩$$|\psi\rangle= \alpha|0⟩+\beta|1⟩$ を用いて計算が行われます。量子ビットの状態ベクトルは、直接知ることはできなくて、測定という操作をすると、状態ベクトルの要素、$\alpha$$\alpha$、$\beta$$\beta$ に応じて、結果が０になる確率と１になる確率が決まります。測定という操作を複数回繰り返すことによって、量子ビットの状態が得られます。測定を行うと、状態ベクトルが変化してしまうのも忘れてはいけません。

図１　古典ビットと量子ビット

１量子ビット

式(1)、(2)で示すように、一つの量子ビットは、２つの基底ベクトル、 $|0⟩,|1⟩$$|0\rangle,|1\rangle$ を用いて量子状態を表します。

これは状態ベクトルや量子ビットと呼ばれます。式（３）でしめす足し合わせた状態のことを重ね合わせ状態と呼びます。$\alpha \text{や}\beta$$\alphaや\beta$ は、複素確率振幅と呼ばれます。もともと、$\alpha \text{や}\beta$$\alpha や\beta$ は、２つの波の重ね合わせの状態を表すものなので、状態ベクトル$|\psi ⟩$$|\psi\rangle$ も波の性質をもち、干渉するということは、量子プログラミングを理解するにあたって重要なことです。

測定

残念ながら、量子力学では、この複素確率振幅は、直接知ることができません。測定という操作をしたときに、０か１かが確率的に決まります。また、測定を行うと、状態ベクトルは、測定結果が $0$$0$ の場合は$|0⟩$$|0\rangle$、$1$$1$ の場合は$|1⟩$$|1\rangle$に変化してしまいます。したがって、量子ビットを初期状態に設定し、操作し、測定を行うということを、繰り返すことによって、$0$$0$ と $1$$1$ の分布を求めます。

状態ベクトル$|\psi ⟩$$|\psi\rangle$を、基底ベクトル$|0⟩,|1⟩$$|0\rangle,|1\rangle$で測定したときに、測定結果が0になる確率 $p_0$$p_{0}$と１になる確率 $p_1$$p_{1}$ は、式（４），（５）のように、複素確率振幅の絶対値の2乗になります。すべての確率の和は1ですので、複素確率振幅は、式（６）の関係を満たします。

このような表現の方法は、ディラック記法と呼ばれ、$|\psi ⟩$$|\psi \rangle$ は、ケットベクトル、$\alpha ,\beta$$\alpha ,\ \beta$ の複素共役（複素数の虚数部の符号を反転したもの）を要素とした行ベクトル $⟨\psi |=({\alpha }^{\ast }{\beta }^{\ast })$$\langle \psi | = (\alpha^* \ \beta^*)$ を、ブラベクトルと呼びます。

$⟨\psi |$$\langle \psi|$ と $|\psi ⟩$$|\psi\rangle$の行列積をとると、内積になります。$⟨\psi ||\psi ⟩$$\langle\psi || \psi\rangle$は、$⟨\psi |\psi ⟩$$\langle\psi | \psi\rangle$とも、表記されます。ブラベクトルとケットベクトルの内積をもちいると、上の条件は、式（７）のように表記できます。

２量子ビット

２個の量子ビットの状態ベクトル $|\psi ⟩$$|\psi\rangle$ は、基底ベクトル $|00⟩,|01⟩,|10⟩,|11⟩$$|00\rangle, |01\rangle,|10\rangle, |11\rangle$ と、４個の複素数を用いて、式（８）、式（９）のように表せます。

ここでも、状態ベクトルの絶対値が量子の存在確率を表しますので、$c_{ij}$$c_{ij}$ は、式（１０）のような関係を満たす複素数とします。

１量子ビットの場合と同様、$c_{ij}$$c_{ij}$ を複素確率振幅と呼び、基底ベクトルで測定すると、確率を得ることができます。

２つの１量子ビット

2 のの1量子ビット $|\psi ⟩,|\varphi ⟩$$|\psi\rangle,|\phi\rangle$ に対して、テンソル積 $|\psi ⟩\otimes |\varphi ⟩$$|\psi\rangle \otimes|\phi\rangle$を、式（１１）、式（１２）のように定義します。

基底ベクトル$|0⟩,|1⟩$$|0\rangle,|1\rangle$に対しては、式（１２）がなりたちます。

つまり、上の2量子ビットの基底ベクトル $|00⟩,|01⟩,|10⟩,|11⟩$$|00\rangle,|01\rangle,|10\rangle,|11\rangle$ は $|0⟩\otimes |0⟩$$|0\rangle \otimes|0\rangle$, $|0⟩\otimes |1⟩,|1⟩\otimes |0⟩,|1⟩\otimes |1⟩$$|0\rangle \otimes|1\rangle,|1\rangle \otimes|0\rangle,|1\rangle \otimes|1\rangle$ と、テンソル積を用いて表すことができるということです。

エンタングルメント（量子もつれ）

重ね合わせ状態に加えて、重要な量子ビット間の関係に、エンタングルメント（量子もつれ）があります。２つの１量子ビットが、上のテンソル積で記述できない状態が存在します。例えば、式（１３）のような状態です。この状態のことをエンタングルメント（量子もつれ）といい、量子コンピュータにおいては、非常に重要な役割を果たします。

量子もつれとは、複数の量子ビットがつくられるときや、相互作用するとき、空間的に近接した場合に起こる物理現象です。量子ビットが大きく離れている場合も含めて、量子ビットの量子状態が他の量子ビットの状態に依存するという不思議な現象です。

状態ベクトルを視覚化するブロッホ球

量子ビットを単位球面上の点として視覚化するのが、Bloch（ブロッホ）球です。状態ベクトル $|\psi ⟩$$|\psi\rangle$ を式（１４）のように変形して、長さ１のベクトルとして表示したものです。ここでは、下の式の右辺全体にかかるグローバル位相と呼ばれる項（$e^{i\lambda }$$e^{i\lambda}$ )を、省略しています。量子ビットを複数もちいた量子回路における量子ビットの回転操作を視覚化するのに用いられます。

図２　ブロッホ球

量子ゲート

古典コンピュータでは、論理ゲートは、AND, OR, NOTなどですが、量子コンピュータでは、量子ビットの回転操作を行う量子ゲートによって、計算が行われます。量子ゲートの回転操作は行列によて表わされます。式（１５）のように、状態ベクトル$|\psi ⟩$$|\psi\rangle$ に行列をかける演算を繰り返すことによって計算が行われます。1量子ビットの操作は、 ２行×２列の 行列 $U$$U$ で表されます。 この行列$U$$U$を量子ゲートと呼び、量子ゲートを状態ベクトルにかけることによって、状態ベクトルに回転操作をすることができます。

この行列 $U$$U$ は、状態ベクトルに$U$$U$をかけた後のベクトルも量子ビットでなければなりませんので、任意の$\alpha ,\beta$$\alpha, \beta$ について$⟨\varphi |U^{†}U|\varphi ⟩=1$$\langle \phi|U^\dagger U|\phi \rangle =1$ がなりたつ必要があり、そのためには、$U^{†}U=U{U}^{†}=I$$U^\dagger U = U U^\dagger = I$ である必要があります。この$U^{†}$$U^\dagger$ は、$U$$U$の随伴行列で、行列$U$$U$ に対して、$U$$U$ の転置と成分の複素共役（実部はそのままで虚部の符号を反転する）をとったものです。このような行列 $U$$U$ は、ユニタリ行列と呼ばれます。以下で説明する各量子ビットゲートは、ユニタリ行列になっており、特に意識する必要はありません。

１量子ビットゲート

１量子ビットゲートには、パウリゲート、アダマールゲート、Sゲート、Tゲート、位相ゲート、回転ゲートなどがあります。

パウリゲート

ブロッホ球のX軸、Y軸、Z軸で、180°回転させるゲートで、Xゲート、Yゲート、Zゲートと呼ばれます。

図３　Xゲート、Yゲート、Zゲート

アダマールゲート

Z軸とX軸にたいして45度傾く軸で180°させるゲートで、Hゲートと呼ばれます。アダマールゲートを使うことにより、量子ビット $|0⟩,|1⟩$$|0\rangle,|1\rangle$ の重ね合わせ状態を作ることができます。

$|0⟩$$|0\rangle$ に、アダマールゲートを作用させたときは、以下のよう操作になります。

図４　アダマールゲート

S, Tゲート

位相を回転させるゲートです。

$S$$S$ ゲートは $\frac{\pi }{2}$$\frac{\pi}{2}$ 、 $T$$T$ ゲートは $\frac{\pi }{4}$$\frac{\pi}{4}$ 位相を回転します。 $S^{†},T^{†}$$S^{\dagger}, T^{\dagger}$ は位相を逆に回転します。

図５　ｓゲート、Tゲート

位相ゲート

量子ビットを任意の位相だけ回転させるゲートです。

図６　位相ゲート

回転ゲート

各軸にたいして量子ビットを回転させるゲートです。

図７　回転ゲート

２量子ビットゲート

２量子ビットゲートには、CNOTゲート、制御ユニタリゲート、SWAPゲートなどがあります。SWAPゲートを除く、2量子ビットゲートは制御タイプです。制御２量子ビットゲートは　最初の量子ビットが |1⟩ の時に、２番目の量子ビットに対して単独の量子ビットゲートを適用します。このとき最初の量子ビットを制御量子ビットと呼び、２番目の量子ビットをターゲット量子ビットと呼びます。

基底ベクトル表記の順番

通常、状態をケットであらわすときに、「1番目」の量子ビット、「2番目」の量子ビットの状態に対応する0と1を左から順番に並べて表記されます。例えば、最初の qubit が |0⟩ 状態で、２番目のものが |1⟩ 状態の場合、その結合状態は |01⟩ です。ところが、今回、プログラミングにもちいるQiskit では、古典コンピューター上でのビット列表現と同じ順序同様にし、測定が実行された後のビット文字列から整数への変換を容易にするという理由から、 最上位ビット (Most Significant Bit - MSB) を左にとり、最下位ビット (Least Significant Bit - LSB) が右にとられます。つまり、「1番目」の量子ビット、「2番目」の量子ビットの状態に対応する0と1を右から順番に並べて表記されます。上の例でいえば、最初の qubit が |0⟩ 状態で、２番目のものが |1⟩ 状態の場合、その結合状態は |10⟩ となります。この点は、注意が必要です。 今回は、Qiskitをもちいたプログラミングを紹介しますので、以下では、Qiskitの基底ベクトル表記を用います。

CNOTゲート

制御Xゲートとも呼ばれます。制御量子ビットが、$|1⟩$$|1\rangle$ のとき、対象となる量子ビットを反転、つまりXゲートを作用させます。MSBを制御ビットとしたときは、以下のようになります。Qiskitでは、 cx(q[1],q[0])。

LSBを制御にした場合は、以下のようになります。Qiskitではcx(q[0],q[1])。

制御ユニタリゲート

CNOTでは、ターゲット量子ビットに、Xゲートを作用させましたが、その他の量子ゲート $U$$U$ を作用させることができます。

制御Zゲート

制御量子ビットが |1⟩ である場合には、Z ゲートがターゲット量子ビットの位相を反転させます。 この行列は、制御量子ビットが、MSBでも、LSBでも、同じになります。

制御アダマールゲート

制御量子ビットが |1⟩ のとき、ターゲット量子ビットに H ゲートを作用させます。制御がLSB量子ビットの時は、以下のようになります。

SWAPゲート

SWAPゲートは2量子ビットの値を交換します。

主な量子ゲートまとめ

ゲート	回路	行列
パウリXゲート		$\begin{array}{c}\left(\begin{array}{ll}0& 1\\ 1& 0\end{array}\right)\end{array}$$\left(\begin{array}{ll}0 & 1 \\ 1 & 0\end{array}\right)$
パウリZゲート		$\begin{array}{c}\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right)\end{array}$$\left(\begin{array}{cc}1 & 0 \\ 0 & -1\end{array}\right)$
アダマールゲート		$\begin{array}{c}\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& 1\\ 1& -1\end{array}\right)\end{array}$$\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1 & 1 \\ 1 & -1\end{array}\right)$
Tゲート		$\begin{array}{c}\left(\begin{array}{cc}1& 0\\ 0& e^{i\frac{\pi }{4}}\end{array}\right)\end{array}$$\left(\begin{array}{cc}1 & 0 \\ 0 & e^{i \frac{\pi}{4}}\end{array}\right)$
CNOTゲート		$\begin{array}{c}\left(\begin{array}{llll}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\\ 0& 0& 1& 0\end{array}\right)\end{array}$$\left(\begin{array}{llll}1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0\end{array}\right)$

図８　主な量子ゲートのまとめ

量子コンピュータはどのように計算をするのか

これまで、量子ビットとはなにか、量子ゲートにはどのようなものがあるかを見てきました。では、量子コンピュータは、量子ビット、量子ゲートを用いて、どのように計算をするのでしょうか。

量子コンピュータでは、まず、すべての状態の重ね合わせを準備します。この状態に、量子ビット間を干渉させるための操作つまり量子ゲートによる操作を順次施します。状態ベクトルつまり量子ビットは、波の性格をもつため干渉するというのがポイントです。量子ビットの間の干渉の結果を測定することによって計算が行われます。複数の量子ビットの状態を同時に操作できるのも量子コンピュータの大きな特徴の一つです。

図９　量子コンピュータの基本

量子ビットの操作は具体的にどうするのか

これまで、量子ビットや量子ゲートを、ベクトルや行列で説明してきましたが、量子コンピュータでは、具体的に、どのように量子ビットを操作しているのかでしょうか。代表的な量子ビットである超電導量子ビットの場合で、見てみたいとおもいます。

超電導量子ビットは、図１０に示すように、２つの超電導電極を薄い絶縁膜でつないだものです。

図１０　超電導量子ビットの構造

これは、非線形インダクタとコンデンサによるLC共振回路とみなせます。量子状態 $|0⟩$$|0\rangle$ と $|1⟩$$|1\rangle$ を持ち、２つの状態の間のエネルギー差のマイクロ波（典型的には４～６GHｚ）を印加することによって、量子状態を、 $|0⟩$$|0\rangle$ から $|1⟩$$|1\rangle$ に変化させることができます。

図１１　超電導量子ビットの操作 (12)

重ね合わせ状態をつくるには、$|0⟩$$|0\rangle$ から$|1⟩$$|1\rangle$への操作を途中でやめます。様々な量子ゲート操作を行うには、それにしたがったパルス操作を行います。ゲート操作を順次繰り返し、最後に測定をすることによって、計算が行われます。また、量子もつれ操作を行うため、量子ビットの間が共振器で結ばれます。

図１２　超電導量子チップの構造 (13)

出所： https://www.ibm.com/blogs/research/2017/11/qubit-explained/

ゲート操作の例

図１３の上の縦軸は電磁波の強度を表しています。実線が１量子ビットの操作で、破線が２量子ビットの操作です。まず、初期状態を準備し、量子回路を通した後、計測を行います。このように、マイクロ波や交流パルスを、順次、印加することによって、量子ビットを操作し、計算を行っていきます。図１３の下に対応するゲートからなる量子回路をしめしています。

図１３ 量子ビットゲート操作の例 (14)

出所: https://journals.aps.org/prx/pdf/10.1103/PhysRevX.6.031007

量子コンピュータは万能

古典コンピュータが、NANDゲートをくみあわせれば任意の計算が可能なように、CNOT、アダマールゲート、Tゲートがあれば任意の計算が可能なことがわかっています。ただし、ノイズの影響により正しい計算ができるかどうか、その計算が早いかどうかは別問題だということに注意が必要です。

図１４　古典コンピュータと量子コンピュータ

量子コンピュータの制約事項・課題

理想的な量子コンピュータは万能であり、古典コンピュータを凌ぐ能力があることはわかっています。ただし、まだまだ、多くの制約事項や課題をかかえています。

• コヒーレンス時間 量子ビットの状態は、古典ビットと比べると、ノイズに弱く、変化してしまいやすいものです。|1⟩という状態にしたはずの量子ビットが、時間が経つと |0⟩ という状態になってしまいます。それにくわえて|0⟩、|1⟩ の重ね合わせ状態がくずれてしまうことも起こります。一定時間経つと、量子ビットの重ね合わせの状態は全てくずれてしまいます。こような量子性が保たれる間の時間をコヒーレンス時間といいます。コヒーレンス時間内に計算を終えないといけません。
• エラー 古典コンピュータにおいても、メモリなどでビットエラーが生じるように、量子コンピュータでも、温度、磁場のゆらぎ、隣接量子ビットの動作などのノイズの影響を受けて、ゲート操作時や測定時にエラーが起こることがあります。

このように、現時点のハードウェアでは、演算中に発生するエラーで正しい計算を行うことが困難なことに加えて、コヒーレンス時間の制約から、深い量子回路の計算も困難です。誤り訂正ありの量子コンピュータが実現されて初めて、実用的な計算を行うことができると言われています。

量子コンピュータの現在

誤り訂正ありの量子コンピュータの実現は、まだまだ先のことです。現在実現されているハードウェアは、エラーの原因となるノイズがある中規模量子デバイス、Noisy Intermediate-Scale Quantumです。今後は、NISQと呼びます。現時点では、性能・機能に制約が多いNISQのうまい使い方の研究開発が進められていたり、誤り訂正ありの量子コンピュータの応用研究がすすめられている段階です。また、現在のNISQで解ける問題の規模は限られており、古典コンピュータで十分高速に解くことができます。

規模が大きくなると古典コンピュータで解くことが困難な問題の一つに、巡回セールスマン問題などの組み合わせ最適化問題があります。セールスマン巡回問題とは、セールスマンが所定の複数の都市を1回だけ巡回する最短経路を計算するというものです。都市数が増加すると計算量が急速に増加しますので、都市数が多くなると現実的な時間で計算が終わりません。このような問題を現実的な時間で解くことが、量子コンピュータに期待されていますが、残念ながら、現状では、量子コンピュータで扱える都市数が少ないとともに、古典コンピュータに勝る計算速度が実現されているわけではありません。今後の量子コンピュータの発展が期待されています。

量子アルゴリズムの種類

量子コンピュータでは、すべての状態の重ね合わせを準備し、順次、量子ゲートの操作を加えていき、量子ビット間の干渉を起こさせ、回答が高確率で得られるような特別なアルゴリズムが必要です。これは古典コンピュータで、様々なアルゴリズムが作られているのと同じです。量子アルゴリズムは、大きく２つに分類されます。

• FTQCアルゴリズム、long-termアルゴリズム 誤り訂正ありの量子コンピュータ（FTQC、Fault Tolerant Quantum Computer）でしか実行が困難なアルゴリズム
• NISQアルゴリズム、量子古典ハイブリッドアルゴリズム 実現しつつあるノイズありの中規模量子デバイス(NISQ,Noisy Intermediate-Scale Quantum)で実行可能なアルゴリズム 一定のエラーがある前提で、量子コンピュータと古典コンピュータを組み合わせたアルゴリズム

FTQCアルゴリズム

FTQCアルゴリズムの代表的なものには、以下のようなものがあります。

• ショア(Shor)のアルゴリズム 素因数分解を解くアルゴリズムです。
• グローバー(Grover)のアルゴリズム 非構造化データを高速に検索するアルゴリズムです。
• 量子フーリエ変換 離散フーリエ変換を行うアルゴリズムです。他の量子計算の要素として用いられることが多いアルゴリズムです。
• 量子位相推定（Quantum Phase Estimation, QPE） 行列の固有値を、量子状態ベクトルの位相として、求めるアルゴリズムです。ハミルトニアン（つまりエネルギー変換行列）の固有値を求めることは、他の多くのアルゴリズムで用いられます。

NISQアルゴリズム

NISQアルゴリズムには、以下のものがあります。

• 量子変分アルゴリズム（Variational Quantum Eigen solver, VQE ) 量子変分原理にもとずいて、行列の最小の固有値を発見的に求めるアルゴリズムです。 量子変分原理は、行列の期待値を最小化することによって、最小の固有値がえられるという原理です。
• QAOA（Quantum Approximate Optimization Algorithm） : 量子近似最適化アルゴリズム VQEの派生系で、量子断熱計算を行うことによって、最適化問題に応用できるアルゴリズムです。
• 量子機械学習アルゴリズム SVM（サポートベクトルマシン）のカーネルや、ニューラルネットワークを量子回路に変えて、機械学習を行うアルゴリズムです。

Python と IBM Quantum, Qiskit でプログラミング

現在利用可能なNISQを用いる前提で、どのような問題が解けるのか、アルゴリズムは実際どのようなものなのかを知るために、現在利用できるフレームワークを用いたプログラムを見てみましょう。

量子コンピューティングのフレームワークとしては、以下のように多くのものがあります。この中で、資料が充実していて、実機も無料で使えることなどを総合的に判断して、今回は、IBM Qiskit とPythonをもちいて量子コンピュータのプログラミングを解説していきます。

• IBM Qiskit (Quantum Information Software Kit)

  • オープンソースの量子ソフトウェア開発フレームワークです。
  • ハードウェアレベルの命令言語 QASM、量子ゲートを操作するローレベルAPI Terra、量子化学、最適化問題、量子機械学習などのアルゴリズムをサポートするハイレベルAPI Aqua、ノイズモデルを含むシミュレーター Aer から成り立っています。
  • 豊富な資料と包括的な教科書が準備されています。ただし、教科書は、一見簡単そうに見えますが、量子物理に精通した人向けにかかれており、初心者にとっては、難易度が高いものとなっています。
  • Qiskitのコミュニティ活動も盛んです。
  • IBMクラウドで、５量子ビットの実機を無料で使うことができます。ただし計算に要する時間は長く、結果がでるまでに数時間かかります。５量子ビット程度の量子回路であれば、シミュレーターの方が圧倒的に早く計算が終了します。

• Google Cirq

  • Qiskit同様、オープンソースのNISQアルゴリズム用のフレームワークです。公式ドキュメントも豊富で、コミュニティも発達してきています。
  • Google Colabでプログラムを実行できます。

• Amazon AWS Braket

  • ２０２０年に発表された比較的新らしいフレームワークです。教育資料もありますが、限られたものです。従量課金方式がとられています。

• Microsoft Q# and Azure Quantum

  • マイクロソフトは、プログラム言語として、Q#を採用しています。シミュレーターをクラウドで提供しています。教育資料は良いものですが、チュートリアルの数が限られています。

• Rigetti Forest

  • PyQuil というNISQアルゴリズムに最適化した独自のフレームワークを提供しています。RigettiやAmazon Braketのクラウドサービスでプログラムを実行できます。教育資料は限られています。

• Xanadu Pennylane

  • Strawberry Fields というフレームワークを提供しています。また、PennyLaneといういろんな量子コンピュータがあつかえるクロス・プラットフォームのPythonライブラリも提供されています。

• blueqat

  • blueqat株式会社の国産の量子コンピューティングpythonフレームワークです。クラウド経由で実機を使うことができます。チュートリアルも充実しています。

• Qulacs

  • 大阪大学の藤井研究室で開発され株式会社QunaSysによって新機能の開発とメンテナンスが行われています。
  • Cirq, PennyLaneなどのサービスのバックエンドとして利用可能なオープンソース量子回路シミュレーターです。大きな量子回路やノイズがある量子回路にも対応しています。

IBM Quantum実行環境

IBM Quantumのアカウント作成

まず、IBM Quantum URL: https://quantum-computing.ibm.com/ にアクセスして、Create an IBMid account で、IBMidを取得してください。

IBM Quantum サインイン画面

取得したIBMidで、サインインするとIBM Quantum Welcome 画面が表示されます。API tokenは、IBM QuantumのAPIにアクセスするときに、必要になります。

IBM Quantum Welcome 画面

ここで、Launch Composer をクリックすると、下のようにGUIで量子回路の作成ができます。ここでは、Pythonを用いたプログラミングを扱いますので、詳細は割愛します。

Welcome画面で、Launch Lab をクリックすると、Jupyter Lab が起動します。標準のJupyter Lab とは、デザインが異なりますが、基本的な使い方は、同じです。Jupyter Labの使い方のは省略します。

Notebookの下のアイコンをクリックすると、Jupyter Lab が立ち上がります。その後、プログラミングを行うという流れです。

IBM Quantum ローカル実行環境

以下のようにローカルにもインストールできますが、各種モジュールのバージョン依存関係などによりエラーが発生する可能性少なからずあります。IBM Quantum lab を使うことを、推奨します。

Anacondaのインストールと環境設定

Qiskitは、Python3.6以降がサポートされており、データサイエンス用のプラットフォームAnacondaが推奨されています。Anacondaのインストールについては、ネット上にも様々な情報がありますので、ここでは省略します。

Anacondaのconda コマンドを使ってQiskit専用の仮想環境を作成することをお勧めします。Anaconda プロンプトを使って以下のように環境を作成します。環境名である ENV_NAME は、各自選んだ適当な名前に変更してください。

xxxxxxxxxx
conda create -n ENV_NAME python=3

新しい環境をアクティブにします。

xxxxxxxxxx
conda activate ENV_NAME

pip をアップグレードします。

xxxxxxxxxx
pip install --upgrade pip

Qiskitパッケージをインストールします。

xxxxxxxxxx
pip install qiskit

可視化機能やJupyter notebookを使うために、 visualization 機能を追加してインストールします。

xxxxxxxxxx
pip install qiskit[visualization]

Qiskit フレームワークの全体像

Terra

量子プログラムを実行する基盤です。いろいろなモジュールから成り立っています。以下、代表的な機能です。

• 量子回路の量子ゲート操作の演算を行い、測定結果を出力します。
• コンパイラ（トランスパイラーとも呼ばれます） 量子ゲート数と量子回路の実行時間を減らすために、量子回路の最適化を行います。
• プロバイダ 量子回路をバックエンドで実行します。
• 様々な形で可視化を行います。

Aer

Aer は、Terraでコンパイルされた回路を実行するためのC++言語で記述された 量子回路のシミュレーターです。Aerは実機で発生するエラーのノイズ・シミュレーションも含みます。３つのシミュレーターがあります。

1. Qasm シミュレーター 量子回路を、ノイズがなし、またはノイズありで複数回実行し、結果を返します。
2. 状態ベクトル（state vector）シミュレーター 量子回路をノイズなしで一回実行し状態ベクトルを返します。実機では、状態ベクトルは直接観測できず計測をしないといけませんが、状態ベクトルシミュレーターでは実行結果の状態ベクトルを知ることができます。あくまでシミュレーションであることに注意が必要です。
3. ユニタリー・シミュレーター。 量子回路をノイズなしで、一回実行し、結果のユニタリー行列を返します。測定やリセット操作を含めることはできません。

Ignis

研究者を対象とした、エラー特定、ゲートの改良、ノイズ存在下での計算などを行うものです。

Aqua

化学や最適化、金融、AI などのアプリケーション用アルゴリズムを実装した量子計算のハイレベルAPIです。量子ゲートの操作などを意識することなく量子化学計算、量子機械学習、最適化、金融への応用などの量子コンピューティングのアプリケーションを作ることができます。

Qiskit とPython を用いた量子コンピュータのプログラミング

Qiskitのプログラムの基本的なステップ

ビルド：問題を解くための量子回路を作成します。

コンパイル：量子コンピュータの実機やシミュレーターで動作させるために回路をコンパイルします。

実行：コンパイルされた回路の演算を実行します。

分析：実行結果の統計処理を行い、実行結果を可視化します。

ビルドとコンパイルは、クラウドのIBM Quantum lab でも、ローカルにQiskit をインストールした環境でも、かまいません。ローカル環境から、実機を用いる場合は、APIを経由して、実機で計算が実行します。

プログラミングの実例

IBM Quantum にサインインして、Launch Lab をクリックし、Qiskit ノートブックを立ち上げます。セルに順次プログラムを入力し、セルを選択して、シフトキーとリターンキーを押すと、セルに入力したプログラムが実行されます。

ステップ 1 : パッケージをインポートする。

必要なパッケージをインポートし、使用するシュミレータを指定します。ここでは、Qasmシュミレータを用います。

xxxxxxxxxx
import numpy as np
from qiskit import QuantumCircuit, transpile # 量子回路のクラスとコンパイラをインポートします。
from qiskit.providers.aer import QasmSimulator #Aerの回路シミュレータです。
from qiskit.visualization import plot_histogram # ヒストグラム作成モジュールです。

simulator = QasmSimulator() # シュミレータとして、qasm_simulatorを設定します。

ステップ 2 : 量子回路を定義する。

xxxxxxxxxx
circuit = QuantumCircuit(2, 2)

２つの量子ビットと２つの古典ビットからなる回路を作成します。おのおのビットはゼロに初期化されます。

ステップ 3 : ゲートを追加する。

ステップ４で可視化される回路を参照するとこのプログラムの意味が、よく理解できると思います。

xxxxxxxxxx
circuit.h(0)                  # 量子ビット０にアダマールゲートで重ね合わせ状態をつくる。
circuit.cx(0, 1)              # 量子ビット０と１の間に、エンタングル状態をつくる 
circuit.measure([0,1], [0,1]) # 計測する。

circuit.h(0) : 初期化された量子ビット０を、アダマールゲートで重ね合わせ状態をつくります。

circuit.cx(0, 1) : 量子ビット０を制御ゲート、１をターゲットゲートとしたCNOTゲートで、量子ビット０と量子ビット１を、エンタングル状態にします。$|\psi ⟩=(|00⟩+|11⟩)/\sqrt{2}$$|\psi\rangle=(|00\rangle+|11\rangle) / \sqrt{2}$

circuit.measure([0,1], [0,1])　: 測定をします。 measure(qubit, cbit)　qubit: 計測する量子ビット、cbit: 計測結果を保存する古典ビット

ステップ ４: 回路を可視化する

draw()を用いて、回路を可視化します。

xxxxxxxxxx
circuit.draw()

ステップ 5: 演算のシミュレーション

シュミレータ Qiskit Aer を用いてシミュレーションを行います。シミュレーションには、qasm_simulator を使います。

xxxxxxxxxx
simulator = QasmSimulator()                       # シュミレータとして、qasm_simulator を使う。
compiled_circuit = transpile(circuit, simulator)  # 上で定義した回路をコンパイルする。
job = simulator.run(compiled_circuit, shots=1000) # シミュレーションを、1000回実行する。
result = job.result()                             # シミュレーション結果を取得する。
counts = result.get_counts(circuit)               # シミュレーション結果から、回数を取得する。
print("\nTotal count for 00 and 11 are:",counts)

シミュレーション結果

xxxxxxxxxx
Total count for 00 and 11 are: {'11': 483, '00': 517}

ノイズを含むqusm_simulotorを使って、１０００回シミュレーションを繰り返して、11となるのが483回、00となるのが、517回という結果になりました。状態ベクトルは、$|\psi ⟩=(|00⟩+|11⟩)/\sqrt{2}$$|\psi\rangle=(|00\rangle+|11\rangle) / \sqrt{2}$ で、演算結果が０になる確率 $p_0=|\alpha {|}^2=1/2$$p_{0}=|\alpha|^{2}=1/2$ 、１になる確率は、 $p_1=|\beta {|}^2=1/2$$p_{1}=|\beta|^{2} =1/2$ ですので、０になる確率が、$517/1000=0.517$$517/1000=0.517$ と期待どおりの結果となりました。０．５にならないのは、ノイズの影響です。

ステップ６: 演算結果の可視化

Qiskit には、多くの可視化のやり方が含まれていますが、ここでは、plot_histogramという関数を用いて可視化を行います。

xxxxxxxxxx
plot_histogram(counts)

実機で実行

上では計算の実行にシュミレータを用いましたが、実機を用いて計算を実行させるには、ここまでのプログラムに、以下のプログラムを追加します。

まず、アカウント設定します。 MY_API_TOKEN には、ログイン時に表示されるAPI token をコビー＆ペーストします。

xxxxxxxxxx
from qiskit import IBMQ
IBMQ.save_account('MY_API_TOKEN')

Welcome画面

次に、設定したアカウントで使用可能な実機のリストを取得します。量子ビット数が１ビットのものがあり、１ビットの実機では、このプログラムは動作しませんので、稼働している実機のうち、量子ビット数が５以上のものをリストアップするようにしています。backendというのは、プログラムを実行する環境のことで、実機とシミュレーターがあります。

xxxxxxxxxx
provider = IBMQ.load_account()
provider.backends(filters=lambda x: x.configuration().n_qubits >= 5
                                    and not x.configuration().simulator
                                    and x.status().operational==True)

稼働している実機のリストが表示されます。実機の稼働状況は変わりますので、リストは取得時点により変わります。

xxxxxxxxxx
[<IBMQBackend('ibmq_bogota') from IBMQ(hub='ibm-q', group='open', project='main')>,
 <IBMQBackend('ibmq_lima') from IBMQ(hub='ibm-q', group='open', project='main')>,
 <IBMQBackend('ibmq_belem') from IBMQ(hub='ibm-q', group='open', project='main')>,
 <IBMQBackend('ibmq_quito') from IBMQ(hub='ibm-q', group='open', project='main')>,
 <IBMQBackend('ibmq_manila') from IBMQ(hub='ibm-q', group='open', project='main')>]

次に、一番空いている実機を選択します。

xxxxxxxxxx
from qiskit.providers.ibmq import least_busy
backend_lb = least_busy(provider.backends(filters=lambda x: x.configuration().n_qubits >= 5
                                    and not x.configuration().simulator
                                    and x.status().operational==True))
print("Least busy backend: ", backend_lb)

実機でプログラムを実行します。

xxxxxxxxxx
from qiskit import execute
backend = backend_lb # 最も空いている実機を使います。

result = execute(compiled_circuit, backend, shots=1000).result() #実機(backend)で1000回実行します。
print(result.get_counts(circuit)) #結果の出力

plot_histogram(result.get_counts(circuit)) #結果のヒストグラムを表示します。

結果

xxxxxxxxxx
{'00': 504, '01': 24, '10': 24, '11': 448}

この結果を見ると、理論上は得られない$|01⟩$$|01⟩$と$|10⟩$$|10⟩$が表れています。これは、実機の演算プロセスで発生したエラーによるものです。量子コンピュータの実機では、このようなエラーがあることが分かります。

汎用量子計算

現在、稼働していて我々が使うことのできる量子コンピュータは、量子ビット数が限られ、量子ビット間の接続にも制限があり、エラー率も高いので、NISQでの実行は難しいと考えられています。したがって、ここでは、汎用量子計算の代表的なアルゴリズムの紹介にとどめます。

量子フーリエ変換（Quantum Fourier Transform, QFT）

離散フーリエ変換を量子コンピュータで実行するアルゴリズムです。離散フーリエ変換は、式（２７）に従って、ベクトル $(x_0,\dots ,x_{N-1})$$\left(x_{0}, \ldots, x_{N-1}\right)$ をベクトル $(y_0,\dots ,y_{N-1})$$\left(y_{0}, \ldots, y_{N-1}\right)$ に変換します。$x_{N-1}$$x_{N-1}$ をAD変換された時間信号とすると、$y_{N-1}$$y_{N-1}$ は周波数スペクトルになります。古典コテンコピュータで、この変換を高速に行うのが、高速フーリエ変換 (fast Fourier transform, FFT) です。古典コンピュータのFFTをさらに高速に実行しようとするのが、量子フーリエ変換です。量子フーリエ変換は、いろいろな量子アルゴリズムの構成要素として用いられます。回路が複雑で、入力状態をどう準備するのかの問題があります。

量子フーリ工変換は、式（２８）に従って、量子状態 $\sum _{i=0}^{N-1}{x}_i|i⟩$$\sum_{i=0}^{N-1} x_{i}|i\rangle$ を、 量子状態 $\sum _{i=0}^{N-1}{y}_i|i⟩$$\sum_{i=0}^{N-1} y_{i}|i\rangle$ に変換します。 $|i⟩$$|i\rangle$ は、整数 $i$$i$ の進数での表示 $i_1\cdots i_n(i_m=0,1)$$i_{1} \cdots i_{n}\left(i_{m}=0,1\right)$ に対応する量子状態 $|i_1\cdots i_n⟩$$\left|i_{1} \cdots i_{n}\right\rangle$ を表します。（例えば、$|2⟩=|0\cdots 010⟩,|5⟩=|0\cdots 0101⟩$$|2\rangle=|0 \cdots 010\rangle,|5\rangle=|0 \cdots 0101\rangle$ などです。)

ケットベクトルを用いて、式（２９）のように表すこともできます。

ユニタリ行列として表現すると、式（３０）のようになります。QFTは、まとめてユニタリ行列として、扱われることがあります。

回路図は以下のようになります。

量子フーリエ変換回路図

$R_l$$R_l$ は、式 (31) の角度 $2\pi /2^l$$2 \pi / 2^{l}$ の一般位相ゲートです。

量子位相推定（Quantum Phase Estimation, QPE）

ユニタリ行列の固有値を求めるアルゴリズムです。多くの量子アルゴリズムの構成要素として用いられます。ユ二タリ行列 $U$$U$が、$U|\psi ⟩=e^{2\pi i\theta }|\psi ⟩$$U|\psi\rangle=e^{2 \pi i \theta}|\psi\rangle$ であたえられた場合に、$\theta$$\theta$ を推定するアルゴリズムです。ここで、 $|\psi ⟩$$|\psi\rangle$ は、固有ベクトルで、$e^{2\pi i\theta }$$e^{2 \pi i \theta}$ は、固有ベクトルに対応する固有値です。固有ベクトルが重ね合わせられた$|\psi ⟩$$|ψ⟩$から固有値を推定するものだと言うこともできます。

量子位相推定、QPEの応用としては、以下のようなものがあります。

1. 素因数分解問題が解ける。 位数という数を推定し、素因数分解をすることができる。
2. 量子系のエネルギー計算ができる。 エネルギーはハミルトニアンと呼ばれる行列の固有値によって与えらます。基底状態のエネルギーは最小の固有値、状態ベクトルは固有ベクトルで与えらます。ハミルトニアン自体はユニタリー行列ではありませんが、工夫を加えることによって、量子位相推定でエネルギーを計算することができます。
3. 連立一次方程式 $Ax=b$$Ax=b$ を解くことができる。 固有値・固有ベクトルが求められるということは、連立一次方程式が解けるということです。また、機械学習への応用も期待されています。

Shorのアルゴリズム

素因数分解を高速にとく量子アルゴリズムです。素因数分解は、自然数 $N(1,2,3,\cdots )$$N (1,2,3,\cdots)$ に対して、$N=pq$$N=pq$ となる素数、$p,q$$p,q$ を求めるというものです。この問題が、古典コンピュータでは高速に解けないことが、RSAなど暗号の基礎となっています。

この素因数問題の解き方は、$N$$N$ を偶数として、以下のようになります。

1. $1$$1$ から $N-1$$N-1$ までの数 $x$$x$ を１つ選びます。
2. $x$$x$ と $N$$N$ の最大公約数を計算します。 最大公約数が１でなければ、その最大公約数が素因数です。 最大公約数が１の場合、次に進みます。
3. $x^r－１$$x^r－１$ を、$N$$N$で割った余りが１になる $r$$r$ を求めます。数式で書くと $x^{r}modN=1$$x^r \ mod \ N = 1$ です。 $r$$r$ が偶数で、 $(x^{r/2}+1)modN\ne 0$$(x^{r/2} +1) \ mod \ N \not = 0$ の条件を満たせば、$x^{r/2}-1$$x^{r/2}-1$と$N$$N$の最大公約数が素因数となります。条件が満たされないなら、１．に戻ります。

この中で、$x^{r}modN=1$$x^r \ mod \ N = 1$ を求める問題を量子アルゴリズムを用いて高速に解くというのがShorのアルゴリズムです。

回路は、量子位相推定と似たものですが、ここでは割愛します。

Shorのアルゴリズムを用いても、現在の量子コンピュータでは、エラーがあることと、量子ビットの数が限られていることから、すぐさま暗号が解読されるレベルの素因数分解の解を求めることはできません。

量子古典ハイブリッドアルゴリズム

量子アルゴリズムだけのプログラムは、現時点では、エラーが多く発生し、また量子ビット数が限らることから、正しく動作しません。そこで、現在も利用可能で、量子ビット数が着実に増えているNISQ（Noisy Intermediate-Scale Quantum technology、小・中規模でノイズを含む量子コンピュータ）の活用方法の開発が盛んに行われています。その多くが、量子コンピュータでは、ビット数の少ないパラメータ付き量子回路を動作させ、古典コンピュータで、そのパラメータを探索・最適化する量子古典ハイブリッドアルゴリズムです。

量子古典ハイブリッドアルゴリズムの基本

回転角度をパラメータ $\theta$$\theta$ にする回転ゲートなどを用いてパラメータ付き量子回路 $U_n$$U_n$ をつくり連結した回路を準備します。その回路に、初期状態を入力し、計測結果からコストを計算した後、コストが最低になるようにパラメータ付き量子回路のパラメータを最適化します。この量子回路は、深層学習におけるモデルと似ています。機械学習、深層学習においては、モデルのパラメータをコストを下げるように、繰り返し、学習を進行させます。深層学習のレイヤーが、各パラメータ付き量子回路 $U_n$$U_n$ に相当します。同様の手法を用いるのが、量子古典ハイブリッドアルゴリズムです。コストが小さくなるようにパラメータを調節するのが、オプティマイザの役割です。コスト関数は、以下にの具体例で述べるように、問題を定式化（値が小さくなるほど望ましい数式で表すこと）したものを用います。最低のコストが求まれば、問題が解けたことになります。

古典コンピュータの深層学習モデルより、量子回路は少ない量子ビット数でも、高い表現力を持つところから、NISQでも実行可能な実用的な計算ができることが期待されています。量子化学計算、組み合わせ最適化問題、機械学習などへの応用が進められています。

ハミルトニアンと期待値

量子状態にハミルトニアンを掛けるとエネルギーになりますので、ハミルトニアンは、量子状態をエネルギーに変換する行列だということができます。量子状態 $|\psi ⟩$$| \psi\rangle$ に対するエネルギーの期待値は、$⟨\psi |H|\psi ⟩$$\langle\psi|H| \psi\rangle$ と書けます。量子は、通常、エネルギーが一番低い状態となっていることが多く、それは、どのような状態なのかが問題になります。また、最適化問題は、問題をエネルギーの形で書き表して、エネルギーの期待値を最小化するという問題と理解することができます。

量子変分アルゴリズム VQE

Variational Quantum Eigen solver（変分量子固有値ソルバー、VQEと略されます）は、量子状態を用いてエネルギーが最小になる基底状態を探索するアルゴリズムです。

分子の性質は、電子の状態によって決まり、その電子の状態は、シュレディンガー方程式 $H|\psi ⟩=E|\psi ⟩$$H|ψ⟩=E|ψ⟩$ を解くことによって、明らかになります。ハミルトニアンは、分子の構造などで決まります。 シュレディンガー方程式を解くということは、ハミルトニアン $H$$H$ の固有値 $E_i$$E_i$ と対応する固有ベクトル $|\varphi i⟩$$|ϕi⟩$を求めることと同じです。このとき固有値 $E_i$$E_i$ は、固有状態 $|\varphi i⟩$$|ϕi⟩$のエネルギーとなります。この問題をとくことは、量子化学計算と呼ばれています。電子の数が増加すると指数的に難しくなるため、厳密に解くことは実質不可能になります。そこで、いかに解を近似するかが課題になります。

式（３２）に示すように、VQEは、様々な量子状態 $|\psi ⟩$$|ψ⟩$ についてそのエネルギーを計算すると、その期待値が基底エネルギー $E_0$$E_0$よりも高くなるという原理に基づいています。

網羅的に様々な状態にたいするエネルギーを計算するのは、効率が悪いため、パラメータ付きの量子状態 $|\psi (\theta )⟩$$|ψ(θ)⟩$を用いて、$⟨\psi (\theta )|H|\psi (\theta )⟩$$⟨ψ(θ)|H|ψ(θ)⟩$ を最小化する $\theta$$\theta$ を発見的に見つけようというのがVQEです。

1. 量子コンピュータ上で量子状態 $|\psi (\theta )⟩$$|ψ(θ)⟩$ を生成します。
2. $⟨H(\theta )⟩=⟨\psi (\theta )|H|\psi (\theta )⟩$$⟨H(θ)⟩=⟨ψ(θ)|H|ψ(θ)⟩$ を測定します。
3. 測定結果をもとに、古典コンピュータによって $⟨\psi (\theta )|H|\psi (\theta )⟩$$⟨ψ(θ)|H|ψ(θ)⟩$ が小さくなるように $\theta$$θ$ を更新し、１．に戻ります。

これを $⟨\psi (\theta )|H|\psi (\theta )⟩$$⟨ψ(θ)|H|ψ(θ)⟩$ が収束するまで繰り返します。

Qiskit プログラム例

ここではH-He+(水素化ヘリウムイオン)の基底エネルギーを求めてみましょう。量子コンピュータの実機では、測定を行うためにパラメータ量子状態に、ハミルトニアンを対角化するためのゲートを追加する必要がありますが、ここでは、簡単のため計算途中の量子状態を得ることができる量子状態シミュレータのプログラムを示します。 以下はdojo.qulacs.orgに掲載されたコードをQiskitで書き直したものです。

まず、必要なモジュールをインポートします。

xxxxxxxxxx
import numpy as np
from qiskit import QuantumCircuit, Aer
from qiskit.opflow import I, X, Z
from qiskit.circuit import Parameter
from qiskit.quantum_info import Statevector
from scipy.optimize import minimize

ハミルトニアンを準備します。H-He間の距離が0.9オングストロームの時のハミルトニアンです。このハミルトニアンの固有状態を計算することによって、$H-H{e}^{+}$$H-He^+$分子の性質がわかります。

xxxxxxxxxx
hamiltonian = ((-3.8505 * I^I) + (-0.2288 * X^I) + (-0.2288 * I^X) + (-1.0466 * Z^I) + (-1.0466 * I^Z) + (0.2613 * X^X) + (0.2288 * X^Z) + (0.2288 * Z^X) + (0.2356 * Z^Z))/2

パラメータ付き量子回路を作成する関数を準備します。

xxxxxxxxxx
def PQC(phi):
    qc = QuantumCircuit(2)
    qc.rx(phi[0], 0)    # 量子ビット０に角度 phi[0]のX軸回転ゲートを追加する。
    qc.rz(phi[1], 0)    # 量子ビット０に角度 phi[1]のZ軸回転ゲートを追加する。
    qc.rx(phi[2], 1)    # 量子ビット1に角度 phi[2]のX軸回転ゲートを追加する。
    qc.rz(phi[3], 1)    # 量子ビット1に角度 phi[3]のZ軸回転ゲートを追加する。
    qc.cx(1, 0)         # 量子ビット1をコントロールビットとして、量子ビット０にCNOTゲートを追加する。
    qc.rz(phi[4], 1)    # 量子ビット1に角度 phi[4]のZ軸回転ゲートを追加する。
    qc.rx(phi[5], 1)    # 量子ビット1に角度 phi[5]のx軸回転ゲートを追加する。
    return qc           # 作成した量子回路を返します。

出来上がった量子回路を、パラメータを、$\theta 0\text{か}\text{ら}、\theta 5$$\theta0 から、\theta5$ として、次のプログラムで表示すると以下のようになります。

xxxxxxxxxx
theta = [Parameter('θ{:}'.format(i)) for i in range(6)] 
PQC(theta).draw()

パラメータ $\theta$$\theta$ に対する期待値 $⟨\psi (\theta )|H|\psi (\theta )⟩$$⟨ψ(θ)|H|ψ(θ)⟩$ を計算する関数を定義します。

xxxxxxxxxx
def cost_func(theta):
    backend = Aer.get_backend('statevector_simulator') # 量子状態シミュレータをつかいます。
    job = backend.run(PQC(theta), shots = 1024)        # パラメータ付き量子回路を、パラメータを θ として計算します。
    result = job.result()                              # 計算結果を取得します。
    q_state = Statevector(result.get_statevector())    # 計算結果の状態ベクトルを取得します。　
    return np.real(q_state.expectation_value(hamiltonian))　# 期待値を返します。

scipy の古典最適化関数 minimize で、基底エネルギをあたえる$\theta$$\theta$ をもとめます。ランダムな初期値から始めて、上で定義したcost_func を呼び出し、損失関数の最小値となるまでパラメータ $\theta$$\theta$ を更新して計算を繰り返します。

xxxxxxxxxx
init = np.random.rand(6)    # ランダムな初期値からスタートします。
res = minimize(cost_func, init, method='Powell')    # 損失関数が最小となる θ を求めます。
cost_func(res.x)    # その時の損失関数の値を求めます。

計算結果です。これで、H-He間の距離が0.9オングストロームの場合の基底エネルギーが求まりました。この問題を厳密（ミルトニアンを対角化）に解いた解、-2.8626207640766816 と小数点以下第三位まで一致しています。

xxxxxxxxxx
-2.8623982023201857

• 参考資料: quantum dojo 第5章 変分量子回路に基づくアルゴリズム

QAOA

QAOA は、Quantum Approximate Optimization Algorithm の略で、量子近似最適化アルゴリズムと呼ばれるものです。組み合わせ最適化問題を解くためのアルゴリズムです。組み合わせ最適化問題には、巡回セールスマン問題、配車ルートの問題、集合被覆問題、ナップサック問題、スケジューリング問題などがあります。これらの問題には、意味のある時間内で正確な解を見つけることが非現実的なものが含まれるため、近似アルゴリズムが必要となります。このアルゴリズムを理解するには、量子断熱計算と呼ばれる知識が必要ですが、ここではQAOAのアルゴリズムの紹介と使い方の説明をすることにとどめます。

• 組み合わせ最適化問題とは、いろんな制約条件のもとで、数多くの選択肢の中から、目的となる指標（コスト関数）を最も小さくする変数の組合せを求めることです。
• 問題をこコスト関数と制約条件からなるモデルで表すことを、定式化するといいます。

まず、0か1 の値をとる$n$$n$個の変数 $z_1,z_2,\cdots ,z_n$$z_1,z_2,⋯,z_n$　について、式（３３）のコスト関数が最小になるような $z_n$$z_n$ の組み合わせを探索する問題として、最適化問題を定式化します。

初期状態をつくるハミルトニアンを $B$$B$、 最適解に相当するハミルトニアンを $C(Z)$$C(Z)$ とし、 $\beta =({\beta }_1,\cdots {\beta }_p),\gamma =({\gamma }_1,\cdots {\gamma }_p)$$\beta=(\beta_1, \cdots \beta_p), \gamma=\left(\gamma_1, \cdots \gamma_p\right)$ をパラメータとして、式（３４）のような量子状態を考えます。

ここで $|+⟩=\frac{1}{\sqrt{2}}(|0⟩+|1⟩)$$|+\rangle=\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$ で、 $U_C（\gamma ),U_B(\beta )$$U_C（\gamma), U_B(\beta)$ は 式（３５）のとおりです。

初期状態 $B$$B$ の代表的なものは、$X$$X$ （パウリ演算子）です。$X$$X$ 以外にも、様々なハミルトニアンが工夫されていて、ミキサ・ハミルトニアンと呼ばれます。

そして、 $F(\beta ,\gamma )=⟨\beta ,\gamma |C(Z)|\beta ,\gamma ⟩$$F(\beta, \gamma)=\langle\beta, \gamma|C(Z)| \beta, \gamma\rangle$ を最小にするような $\beta ,\gamma$$\beta, \gamma$ を探索することで、元々の 最適化問題の答えを求めようとするのが、QAOAアルゴリズムです。

Qiskit プログラム例 Maxcut

手順は、以下の通りです。

1. 問題の定式化
2. コスト関数に変換
3. 量子回路の作成
4. 量子回路と古典オプティマイザをもちいて最適なパラメータを求める
5. 得られたパラメータでQAOA回路をつくり結果を分析する。

まず、必要なパッケージをインポートします。

xxxxxxxxxx
import networkx as nx
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister
from qiskit import Aer, execute
from qiskit.circuit import Parameter

4個の頂点 $V=0,1,2,3$$V={0,1,2,3}$ と、4つの辺 $E=(0,1),(1,2),(2,3),(3,0)$$E={(0,1),(1,2),(2,3),(3,0)}$ を持つグラフを考えます。このグラフのを２つに分割するときに、分割される辺の数の最大値を求める問題を考えます。これは、Maxcut問題と呼ばれています。

Python のnetworkxというモジュールでグラフを作成し描画するとこのようになります。

xxxxxxxxxx
G = nx.Graph()　　　　　　　　　　　　　　　　　　　　　　　    # 空のグラフを用意する。
G.add_nodes_from([0, 1, 2, 3])　　　　　　　　　　　　　　   # 頂点を追加する。
G.add_edges_from([(0, 1), (1, 2), (2, 3), (3, 0)])      # 辺を追加する。
nx.draw(G, with_labels=True, alpha=0.8, node_size=500)  # グラフを描画する。

このグラフ $G(V,E)$$G(V,E)$ において、頂点を２つに分割して、一つのグループにたいして、$x_i=0$$x_i=0$ 、他のグループに対して、$x_i=1$$x_i=1$ を与えると、分割した各々のグループの頂点を結ぶ辺の数をコスト関数とすれば、式（３６）のように、問題を定式化できます。

このような形式は、QUBO (Quadratic Unconstrained Binary Optimization) と呼ばれます。

これをハミルトニアンに、変換するには、$x_i=(1-Z_i)/2$$x_i = (1-Z_i)/2$ 　($Z_i$$Z_i$ は、パウリ$Z$$Z$演算子）とすればよく、式（３７）のようになります。

この問題のグラフでは、式（３８）のようになります。

$2$$2$ は、定数ですので、コスト関数としては、$2$$2$ を考慮する必要はありません。したがって、コスト関数は式（３９）のようになります。

初期状態としては、式（４０）の通り、$X$$X$ （パウリ演算子）をもちいます。

xxxxxxxxxx
def maxcut_obj(x, G):
    # ２つのグループに分割した結果xに対して、辺の数を返す関数です。
    # 辺の数は負の数として返され、それを最小化します。
    # 頂点が同じグループに属さないとき、その頂点間の辺がカットされます。
    # 下のcompute_expectationで使われます。
    obj = 0
    for i, j in G.edges():
        if x[i] != x[j]:
            obj -= 1            
    return obj

def compute_expectation(counts, G):
    # 計測結果countsから、辺の数の期待値を計算する関数です。
    # bitstring は、どのようにグループ分けするかの結果、counts は、そのグループ分けがされた回数です。
    # 繰り返し回路を実行した結果の辺の数の期待値（平均値）を返します。
    avg = 0
    sum_count = 0
    for bitstring, count in counts.items():
        obj = maxcut_obj(bitstring, G)
        avg += obj * count
        sum_count += count
    return avg/sum_count

def create_qaoa_circ(G, theta):
    # グラフGに対応し、パラメータθ(β,γ)をもつQAOA回路をつくる関数です。
    nqubits = len(G.nodes())     # ノードの数
    p = len(theta)//2            # ユニタリ行列の繰り返し回数
    
    qc = QuantumCircuit(nqubits) # nqubitsの量子回路を準備する。
    beta = theta[:p]             # パラメータβ
    gamma = theta[p:]            # パラメータγ
    
    # 初期状態を設定するためアダマールゲートを追加して重ね合わせ状態をつくります。
    for i in range(0, nqubits):
        qc.h(i)                  # アダマールゲートを追加する。
    
    for irep in range(0, p):
        # 問題のハミルトニアン UC(γ) 
        for pair in list(G.edges()):
            qc.rzz(2 * gamma[irep], pair[0], pair[1])

        # ミキサ・ハミルトニアン UB(β)
        for i in range(0, nqubits):
            qc.rx(2 * beta[irep], i)         
    
    # 計測します。
    qc.measure_all()
    return qc

def get_expectation(G, p, shots=512):
    # グラフGをユニタリ行列の繰り返し回数pで計算をおこないます。
    
    backend = Aer.get_backend('qasm_simulator')
    backend.shots = shots
    
    def execute_circ(theta):
        # 上のcreate_qaoa_circ関数をもちいて量子回路をつくります。
        qc = create_qaoa_circ(G, theta)  
        # 作成した量子回路を動作させ指定したショット数だけ計算を実行し、各量子ビットのカウント値を求めます。
        counts = backend.run(qc, seed_simulator=10, nshots=512).result().get_counts()
        # カウント値から、辺の数期待値を計算します。
        # 辺の数はマイナスで数えられているので、この期待値を最小化することで問題を解きます。
        return compute_expectation(counts, G)
    return execute_circ

作成される回路は次の通りです。

Maxcut 回路図

次に、scipyの古典最適化モジュールで、最適値を求めます。

xxxxxxxxxx
from scipy.optimize import minimize

expectation = get_expectation(G, p=1) # ユニタリ行列の繰り返し回数 1 での期待値計算を対象とします。
res = minimize(expectation, [1.0, 1.0], method='COBYLA') 
# β、γの初期値を、1.0, 1.0 として、辺の数の期待値を最小化する最適値を求めます。
# ここでは、COBYLA という最適化手法を用います。

res # 結果を表示します。

fun: -2.994140625 ​ maxcv: 0.0 message: 'Optimization terminated successfully.' ​ nfev: 30 ​ status: 1 success: True ​ x: array([1.9793337 , 1.16663483])

最適化が終了し、x すなわち $\beta \text{と}\gamma$$\betaと\gamma$ の最適値 [1.9793337 , 1.16663483] が求まりました。

次に、最適化されたパラメータ $\beta \text{と}\gamma$$\betaと\gamma$ で、QAOA回路をつくり、結果の分析をします。

xxxxxxxxxx
from qiskit.visualization import plot_histogram

backend = Aer.get_backend('aer_simulator')      #バックエンドとしてaer_simulator を用います。
backend.shots = 512

qc_res = create_qaoa_circ(G, res.x)             # 最適化されたパラメータ β、γ でQAOA回路を作ります。
counts = backend.run(qc_res, seed_simulator=10).result().get_counts()　# 量子回路計算を実行します。

plot_histogram(counts)

ビット列「0101」と「1010」（すなわち頂点０と頂点２、頂点３と頂点１でカットする）の確率が最も高く、図２５に示すように、最大の辺の数が、4つである分割となっています。

Maxcut の計算結果

このプログラムを見て、複雑だと思われた読者も多いと思いますが、心配はいりません。IBM Qiskit には、QAOAのゲート回路などをプログラミングすることなく問題を解くことができるハイレベルAPIが準備されています。

[上記のプログラムは、IBM qiskit テキストブックの Solving combinatorial optimization problems using QAOA https://qiskit.org/textbook/ch-applications/qaoa.html の一部で、コメントを日本語に修正しています。]

Qiskitを用いてVQE,QAOA で実問題を解く

ここでは、IBM quantum challenge fall 2021(https://github.com/qiskit-community/ibm-quantum-challenge-fall-2021) から、興味深い実問題の解法を紹介します。Python と Qiskit を用いたプログラムです。Qiskit のハイレベルAPIを用いることによって、実問題を簡単に解くことができます。

ポートフォリオ最適化

ポートフォリオとは、株式などの金融資産の組み合わせのことです。ポートフォリオ最適化のゴールは、リスクを最小化し、リターンを最大化することですが、リスクとリターンは通常トレードオフの関係にあります。積極的に大きなリターンを得ようとするなら、リスクが大きい金融資産を多く所有し、手堅く運用したいなら、リスクが小さい金融資産に多くを配分することになります。最適な配分を決めることは、そう簡単ではありません。

最適な配分を決める方法として、現代ポートフォリオ理論というものがあります。リスクを避けながら、いかにリターンを大きくするかという方針にもとづいて最適な配分を行うというものです。

そのためには、以下のことを行います。

• 過去のデータに基づいて、ポートフォリオの銘柄の平均、分散と相関を計算する。
• リスクを抑えつつリターンを大きくする銘柄の重みと期待リターンを計算する。

四銘柄のポートフォリオ最適化問題

こでは、株式４銘柄の、リスクとリターンの間のトレードオフを最小化する配分を見つけることを問題とします。

問題の定式化

$x_i$$x_i$: 銘柄の配分、保有するときは１、保有しない場合は０

$n$$n$: 保有対象の銘柄数

${\sigma }_{i,j}$$\sigma_{i, j}$:２つの銘柄の値動きが、どう相関しているかの指標である共分散です。２つの銘柄が同時に上がる関係にあるなら正の値、一つが上がればもう一つが下がる関係にあるなら、負の値をとり、値の絶対値が、その度合を示します。リスクは、同じ値動きをする銘柄を保有するときに高く、逆の銘柄を保有するときは、低くなるような指標です。行列の形にしたものを、共分散行列 $\Sigma$$Σ$ と表記します。

$q$$q$: リスクファクタ（リスク許容度）、０から１の間の値をとり、どの程度のリスクを取るか、取れるかを決める指標です。

${\mu }_i$$\mu_i$: 各銘柄の期待リターン

$B$$B$: 保有する銘柄の数

として、（４１）のように定式化できます。第一項が、ポートフォリオのリスク（保有する銘柄間の共分散を足し合わせるたものに保有する銘柄数を足したもの）にリスクファクタを掛けたもの、第二項がポートフォリオのリターンを表しています。つまり、リスクを抑えつつ、リターンを大きくしようとすることを数式化したものです。。

ゴール

どの銘柄を選ぶか、選ばないかの $x_i$$x_i$ を見つけることです。

仮定

問題を簡単にするため、以下の仮定を置きます。

• すべての銘柄の価格は同じ。
• 予算 $B$$B$ を使い切る。

ライブラリのインポート

xxxxxxxxxx
from qiskit import Aer
from qiskit.algorithms import VQE, QAOA, NumPyMinimumEigensolver
from qiskit.algorithms.optimizers import *
from qiskit.circuit.library import TwoLocal
from qiskit.utils import QuantumInstance
from qiskit.utils import algorithm_globals
from qiskit_finance import QiskitFinanceError
from qiskit_finance.applications.optimization import PortfolioOptimization
from qiskit_finance.data_providers import *
from qiskit_optimization.algorithms import MinimumEigenOptimizer 
from qiskit_optimization.applications import OptimizationApplication
from qiskit_optimization.converters import QuadraticProgramToQubo
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import datetime
import warnings
from sympy.utilities.exceptions import SymPyDeprecationWarning
warnings.simplefilter("ignore", SymPyDeprecationWarning)

時系列データ（金融データ）の生成

RandomDataProviderクラスを使って、４銘柄のランダムな値動きデータを作成します。バジェットを２として、４銘柄から２銘柄を選択することとし、リスクファクタを0.5にします。

xxxxxxxxxx
# 銘柄とリスクファクターのパラメータを設定します。
num_assets = 4     # 銘柄数　４
q = 0.5            # リスクファクター 0.5
budget = 2         # バジェットは、２
seed = 132     

# ４銘柄のランダムな値動き時系列データを作成します。
stocks = [("STOCK%s" % i) for i in range(num_assets)]
data = RandomDataProvider(tickers=stocks,
                 start=datetime.datetime(1955,11,5),   
                 end=datetime.datetime(1985,10,26),  
                 seed=seed)
data.run()

xxxxxxxxxx
#４データをプロットします。
for (cnt, s) in enumerate(data._tickers):
    plt.plot(data._data[cnt], label=s)
plt.legend()
plt.xticks(rotation=90)
plt.xlabel('days')
plt.ylabel('stock value')
plt.show()

期待リターンと共分散行列の計算

get_period_return_mean_vector() メソッドで、期待リータンが計算できます。

xxxxxxxxxx
mu = data.get_period_return_mean_vector()   # おのおのの銘柄の期待されるリターンを返します。
print(mu)

[1.59702144e-04 4.76518943e-04 2.39123234e-04 9.85029012e-05]

４銘柄に対しての期待するリターンが、求まりました。次に、共分散行列の計算を行い、表示します。

xxxxxxxxxx
# 共分散行列 Σ のプロット
sigma = data.get_period_return_covariance_matrix() # ４つの銘柄の共分散行列を返します。
fig, ax = plt.subplots(1,1)
im = plt.imshow(sigma, extent=[-1,1,-1,1])
x_label_list = ['stock3', 'stock2', 'stock1', 'stock0']
y_label_list = ['stock3', 'stock2', 'stock1', 'stock0']
ax.set_xticks([-0.75,-0.25,0.25,0.75])
ax.set_yticks([0.75,0.25,-0.25,-0.75])
ax.set_xticklabels(x_label_list)
ax.set_yticklabels(y_label_list)
plt.colorbar()
plt.clim(-0.000002, 0.00001)
plt.show()

共分散行列 $\Sigma$$Σ$ は、２つの銘柄の期待リターンがどのように関連するかを表します。対角成分（図の黄色で示された値）は、ある銘柄自身との関係を示しています。 共分散行列の見方は次の通りです。

• 2つの銘柄が同時に増加・減少すれば、共分散の値は正の値になります。
• 一方が上昇し、他方が下降した場合、共分散の値は負になります。

違う動きをする銘柄に投資することはリスクを減らす効果があります。

Qiskit Finance のポートフォリオ最適化クラス

PortfolioOptimization クラスを使うと、$\mu ,\mathrm{\Sigma },q,B$$\mu , \Sigma, q, B$ から自動的に、QUBOで問題が定式化されます。ポートフォリオの資産選択の変数が、バイナリ変数（０，１）なので、「bounds = None」と設定します。

xxxxxxxxxx
portfolio = PortfolioOptimization(expected_returns=mu, covariances=sigma, risk_factor=q, budget=budget, bounds=None)
qp = portfolio.to_quadratic_program()    # μ,Σ、qからQUBOを作成し返します。
print(qp)

定式化した結果が出力されます。

Minimize （最小化するコスト関数、x_0, x_1,x_2,x_3 が変数） obj: - 0.000159702144 x_0 - 0.000476518943 x_1 - 0.000239123234 x_2 - 0.000098502901 x_3 + [ 0.000048831990 x_0^2 - 0.000002157372 x_0x_1 - 0.000004259230 x_0x_2 + 0.000001413200 x_0x_3 + 0.000997360142 x_1^2 + 0.000007031887 x_1x_2 + 0.000000737432 x_1x_3 + 0.000287365468 x_2^2 + 0.000006416382 x_2x_3 + 0.000192316728 x_3^2 ]/2 Subject To　（制約条件） c0: x_0 + x_1 + x_2 + x_3 = 2

Bounds　（値の範囲） 0 <= x_0 <= 1 0 <= x_1 <= 1 0 <= x_2 <= 1 0 <= x_3 <= 1

Binaries　（バイナリ変数） x_0 x_1 x_2 x_3

Minimum Eigen Optimizer（最小固有値オプティマイザー）

このポートフォリオ最適化問題は、ハミルトニアンの基底状態を求める問題として解くことができます。ハミルトニアンを考えると、シミュレーションしたいシステムのエネルギー関数と考えることができます。そのためには、図に示すように、QUBOで定式化したものを、イジングモデルと呼ばれる数学モデルで表現し、さらにハミルトニアンに変換するという複雑な計算を行わないといけません。Qiskit は、これらの変換を自動的にやってくれます。MinimumEigenOptimizer（最小固有値オプティマイザー）が、QUBOをハミルトニアンに変換し、VQEやQAOAなどを呼び出して基底状態を計算し、最適化の結果を返します。

最適化問題の解法の流れ

古典コンピュータの最適化アルゴリズムで問題を解く

この規模の問題は、古典コンピュータの最適化アルゴリズムで解くことができます。

xxxxxxxxxx
exact_mes = NumPyMinimumEigensolver()    # Numpy最小固有値ソルバーを使います。
exact_eigensolver = MinimumEigenOptimizer(exact_mes)    #  最小固有値オプティマイザーに、Numpy最小固有値ソルバーを使う設定をします。
result = exact_eigensolver.solve(qp)    # Step 4. で作成したQUBOを解きます。
print(result)