この記事の目的

量子コンピュータを応用したアルゴリズムの一つである**VQE(Variational Quantum Eigensolver)**について、その基本からblueqatを使用したプログラミングまでの解説を行います。

VQEとは

VQEは Variational Quantum Eigensolver の略で、変分量子固有値ソルバーと言います。ソルバーとは、数式や数値解析を用いて最適な解を求める機能またはソフトウェアのことです。
よってVQEとは、変分原理を用いて量子力学の理論から最適な固有値を求めるもの、という意味になります。

変分原理と量子力学および固有値の関係について順を追って解説します。

VQEアルゴリズムを使用することで、化学分野では材料科学計算や量子化学計算へ、また、組合せ最適化問題の分野への応用などに期待されています。

それでは、ここから変分量子固有値ソルバーという言葉を手綱に、基本からしっかり解説していきましょう。

固有値とは何か?

まず、固有値について解説します。固有値は数学の分野の話になります。
固有値の数学的な定義は以下です(難しく感じる場合はこの定義部分を一旦飛ばして読んで下さい)。

固有値・固有ベクトル・固有空間

Fをn次元ベクトル空間Vの線形変換、Aをn次正方行列とするとき、\\
　1.　F(\mathbf{x})=\lambda\mathbf{x}　かつ　\mathbf{x}\ne\mathbf{0}　なる　\mathbf{x}\in V　が存在するとき、\\
\lambdaをFの固有値、\mathbf{x}を固有値\lambdaに対するFの固有ベクトル、\\
W(\lambda)=\lbrace\mathbf{x}|F(\mathbf{x})=\lambda\mathbf{x}\rbrace\\
を固有値\lambdaに対するFの固有空間という。\\
　2.　A\mathbf{x}=\lambda\mathbf{x}　かつ　\mathbf{x}\ne\mathbf{0}　なる　\mathbf{x}\in \mathbb{C}^{n}　が存在するとき、\\
\lambdaをAの固有値、\mathbf{x}を固有値\lambdaに対するAの固有ベクトル、\\
W(\lambda)=\lbrace\mathbf{x}|A\mathbf{x}=\lambda\mathbf{x}\rbrace\\
を固有値\lambdaに対するAの固有空間という。

平たく説明すると、$A=(a_{jk})$をn x n行列、$\mathbf{x}$を列ベクトルとして、以下の一次方程式

A\mathbf{x}=\lambda\mathbf{x}

について、ベクトル$\mathbf{x}$に行列$A$を作用させても、もとのベクトル$\mathbf{x}$の$\lambda$倍になる特別なベクトル$\mathbf{x}$を固有ベクトル、その根となる$\lambda$を固有値と言います。

固有値は一般に複素数です。また、固有値は重根となる場合があり、その場合を縮退していると言います。

量子力学について次の節で詳しく説明しますが、量子力学においてはしばしば、固有ベクトルを固有関数、または固有状態と言います。
量子力学では、固有値と固有関数を以下のように記述します。

Aを物理量とするとき、\\
A\varphi(x)=\alpha\varphi(x)　　...式(1)\\
ならば、\alphaをAの固有値、\varphi(x)を固有関数(固有状態)という。

量子力学とVQE

量子コンピュータは基本的に量子力学の理論を応用したコンピュータと言えます。

量子力学とは、物質を構成する原子や電子、イオン化されたこれらの粒子、さらにこれらを構成するさらに小さい各種の素粒子(まとめて量子と言います)の挙動(振る舞い)を説明する理論です。
これらの量子の振る舞いを記述する最も基本的な式がシュレーディンガーの方程式です。
以下に時間依存しないシュレーディンガー方程式を紹介します。

-\frac{\hbar^2}{2m}\frac{d^2\varphi(x)}{dx^2}+V(x)\varphi(x) = E\varphi(x)\\
\hbar：プランク定数\\
V：ポテンシャルエネルギー\\
E：粒子のエネルギー

上記の式について、$-\frac{\hbar^2}{2m}\frac{d^2}{dx^2}+V(x)$をハミルトニアン$H$とおくと以下のように簡潔に記述することができます。

H\varphi(x)=E\varphi(x)　　...式(2)

上式を見ると、前節で説明した固有値・固有ベクトルの関係である式(1)と同じ形をしていることがわかります。シュレーディンガー方程式とはハミルトニアン$H$の固有関数$\varphi(x)$、固有値$E$を求める式であると言えます。

シュレーディンガーの方程式【式(2)】を応用して粒子の状態である固有値(エネルギー状態)を求めることがVQEの理論の裏付けとなります。

期待値の公式

変分原理についての説明に入る前に、期待値の求め方について説明します。期待値はVQEでは以下のように書かれますが、なぜそのように書くのかを確認します(この部分の説明は飛ばしても構いません)。

期待値 = 〈\psi|A\psi〉　　...式(3)\\
Aは物理量

一般に物理量とは、位置や運動エネルギーなどの測定可能な量のことです。物理量の測定値は実数となります。

ここで以下の量子力学の公理を確認します。

量子系ではHilbert空間(ヒルベルト空間)上で考察を行います。
物理量はHilbert空間上のエルミート演算子で表現され、測定値はその固有値とします。量子状態$\psi$の下で物理量Aを測定するとき、測定値が固有値$a$となる確率$Pr$は以下の式で与えられます。

Pr(A=a | |\psi〉)=〈\psi|P_a\psi〉　　...式(4)\\
P_aはAの固有値aに対する固有射影演算子

式(4)は確率を求める式です。この式をBornの確率規則と言います。
上式は、エルミート演算子Aに対応する物理量の測定を行った時のAの固有値のうち一つが得られる確率を意味しています。

ここでは有限次元を考えますので、Hilbert空間は一般のベクトル空間と考えて良いです。
固有値$a$は実数値です。また、$P_a$は固有値aに対する固有射影演算子です。

固有射影演算子とは、固有値を与える射影演算子で、さらに射影演算子とは以下を満たす演算子です。

P_a=P_a†=P_a^2　　...式(5)\\
\sum_{a}^{}P_a=\mathbb{I}　　...式(6)

式(5)を用いて、式(4)について以下が言えます。
これは、$Pr$が0以上を取ることを意味します。

〈\psi|P_a\psi〉=〈\psi|P_a^2\psi〉=〈P_a\psi|P_a\psi〉≧0　　...式(7)

同様に、式(6)を用いて式(4)について以下が言えます。
これは、各確率の和を取ると1になることを意味します。

\sum_{a}^{}〈\psi|P_a\psi〉=〈\psi|(\sum_{a}^{}P_a)\psi〉=〈\psi|\mathbb{I}\psi〉=||\psi||^2 = 1　　...式(8)

式(7)式(8)より、式(4)は確率分布となることが分かります。
上記の公理を用いて期待値を算出することができます。

状態$|\psi〉$のもとで物理量Aを測定するときの測定値の期待値Eは固有値$a$と確率$Pr$の積$aP_a$の和を取り、

\begin{align}
E_\psi[A] &= \sum_{a}^{}aPr(A=a||\psi〉)\\
&=\sum_{a}^{}a〈\psi|P_a\psi〉\\
&=〈\psi|(\sum_{a}^{}aP_a)\psi〉\\
&=〈\psi|A\psi〉　　...式(9)
\end{align}

上式において、以下のスペクトル分解定理を使用しました。

任意の正規演算子は以下のように表すことができる。\\
A=\sum_{a∊\sigma(A)}^{}aP_a

よって、式(3)の期待値を求めることができました。
これを期待値の公式と呼びます。

変分原理を理解する

変分原理と言っても古典力学や電磁気学、量子力学などいくつかの分野で定義された異なる理論がありますが、ここでは量子力学での変分原理であるリッツの変分原理を取り上げます。

期待値の公式【式(3)】より、期待値を状態$|\psi〉$と物理量Aで求めことができることがわかりました。
すでに説明した通り、物理量Aは位置やエネルギーなどさまざまな量を表す演算子です。ここで、その物理量の一つでありエネルギー量を表す演算子の$H$(ハミルトニアン)をAの代わりに用います。

ハミルトニアン$H$と量子状態$|\psi〉$にあるときの期待値は、ある固有値(固有エネルギー)$E_0$を用いて、以下のように書き表すことができます。

〈\psi|H|\psi〉≧E_0　　...式(10)

これは、量子力学の観点から以下のように説明できます。

すなわち、さまざまな状態において、電子などの粒子は量子力学の法則に従いいくつかのエネルギー状態を取ります。
最もエネルギーの低い状態を基底状態と呼びます。また、エネルギーが高くなるにつれて、順に第一励起状態、第二励起状態、さらに上位に状態が変化します(下図)。

よって、以下が成り立ちます。
適当な境界条件を持つ任意の状態$|\psi〉$に対するハミルトニアン$H$の期待値は、基底状態のエネルギー$E_0$よりも常に大きいか等しい。
これは式(10)が成立することを意味します。
これが変分原理です。

ψへのパラメータθの導入

上記の変分原理を応用し、エネルギーが最小となる固有値$E_0$を求めるためにはどうすればよいでしょうか。この答えは単純で、様々な異なる状態の$\psi$を与えたときに、期待値が最も小さくなるときの$\psi$を探し出せば良いわけです。
そのときの期待値が固有値$E_0$となります。

様々な異なる$\psi$を与えるために、新たに角度のパラメータ$\theta$を用意します。
そして、$\psi$がパラメータ$\theta$に依存する一種の関数とみなします。
期待値は以下のように書き直すことができます。

〈\psi(\theta)|H|\psi(\theta)〉

VQEアルゴリズムを用いたプログラミングの実践

変分原理を理解した上で、VQEを用いたプログラミングに入りたいと思います。
プログラミングには、量子コンピューターシミュレータとして評価の高い、blueqat株式会社(東京・本郷)が提供するblueqat(ブルーキャット)を使用します。

blueqat (https://blueqat.com/)

また、解く問題としては組合せ最適化問題の一つである頂点被覆問題を取り扱います。
今回解く頂点被覆問題については、こちらをご覧下さい。

Blueqatで簡単な頂点被覆問題を解く
https://qiita.com/ttabata/items/b4ab222ef9c5ee99233c

(1)変分原理によるVQEアルゴリズムの計算手順

変分原理を用いたVQEアルゴリズムは以下の手順で計算を行います。

　①パラメータ$\theta$を用いた量子回路を準備(ansatz：アンザッツ)
　②量子ビットを測定し、期待値を算出する
　③古典最適化計算により、期待値を最小にするパラメータ$\theta$を探索する

VQEでは、量子コンピュータによる期待値の測定と、古典最適化計算(古典コンピュータの計算)によるパラメータ$\theta$の探索処理の両方を用いて解を算出します。

処理①~③を模式図で示すと以下となります。こちらはVQEの論文からの引用です。
引用：(https://arxiv.org/pdf/1304.3061.pdf)

処理①に記載のansatzとは、量子計算に必要な量子状態$|\psi〉$を生成するための量子回路のことで、任意回転ゲートと入力パラメータ$\theta$を用いて量子状態を生成します。この量子状態を用いて期待値を計算します。

(2)blueqatによるVQEアルゴリズム計算のプログラミング

blueqatには上記の計算手順①~③を内部関数で実行してくれる機能が備わっており、プログラマーはパラメータや期待値計算、最小値の探索に手をかけなくても良くなっています。これはblueqatが優れている点の一つと言えます。

プログラマーがすべきことは、パラメータとハミルトニアン$H$を定義することのみです。

では、さっそくパラメータ及びハミルトニアン$H$の定義を行いましょう。
前掲の例題である頂点被覆問題から、ハミルトニアンを以下とします。

H=5-\frac{3}{2}q_0-\frac{3}{2}q_1-\frac{3}{2}q_2-\frac{5}{2}q_3-\frac{1}{2}q_4\\
+q_0q_1+q_0q_3+q_1q_2+q_2q_3+q_3q_4

例題では変数を$x$としていますが、ここでは量子ビット(quantum bit)の意味で、アルファベット$q$を使用しました。また、インデックス番号を0始まりとしています。

基本的にはこのハミルトン式を組み込めば良いです。
プログラムのソースコードは以下となります。

sample_vqe.py

from blueqat import Circuit
from blueqat.pauli import X, Y, Z, I
from blueqat.pauli import qubo_bit as q
from blueqat.vqe import AnsatzBase, Vqe

class QubitAnsatz(AnsatzBase):
    def __init__(self, hamiltonian):
        super().__init__(hamiltonian, 10)
        self.step = 1
    
    def get_circuit(self, params):
        a, b, c, d, e, f, g, h, i, j = params   # パラメータの定義1
        return Circuit().ry(a)[0].rz(b)[0].ry(c)[1].rz(d)[1].ry(e)[2].rz(f)[2].ry(g)[3].rz(h)[3].ry(i)[4].rz(j)[4]   # パラメータの定義2

H = 5 - 3*q(0)/2 - 3*q(1)/2 - 3*q(2)/2 - 5*q(3)/2 - 1*q(4)/2 + q(0)*q(1) + q(0)*q(3) + q(1)*q(2) + q(2)*q(3) + q(3)*q(4)   # ハミルトニアンの定義
H = H.to_expr().simplify()

runner = Vqe(QubitAnsatz(H))
result = runner.run()

print('Result by VQE')
print(result.most_common())

ソースコード中で「パラメータの定義1」と「パラメータの定義2」の記述を行いました。
今回の例題では、ハミルトニアン内に含まれる量子ビット(q)がq0～q4まで5個あります。それぞれの量子ビットに対して回転ゲートry,rzを適用するため、回転ゲートの総数は10個となり、それぞれにパラメータを設定するため、パラメータ数も10個です。
よって、「パラメータの定義1」に10個分のパラメータa～jを宣言しています。

「パラメータの定義2」では上記に定義した10個のパラメータを用いて、回転ゲートを宣言します。量子ビット(q)がq0～q4まで5個分あるので、回転ゲートry,rzもそれぞれ、ry(0)～ry(4)の5個分と、rz(0)～rz(4)の5個分を宣言します。

「ハミルトニアンの定義」にハミルトニアンを定義します。
プログラミングは以上です。

実行結果は以下となりました。結果の見方は例題と同じです。

q0～q4は以下の値を取ります。

q0	q1	q2	q3	q4
0	1	0	1	0

ここでは、0,1はクラス分けを意味しています。
0を取る頂点のクラスと、1を取る頂点のクラスに分かれます。

このクラス分けが今回の問題の解答となります。

関連情報

【ほぼ決定版】量子コンピュータVQE/QAOAセミナーまとめ
https://qiita.com/YuichiroMinato/items/3f68b9fed67ff3a01370

VQE algorithm (Variational Quantum Eigensolver)
https://qiita.com/KeiichiroHiga/items/c900d430dc33d4342b24

VQE/QAOAの期待値測定の部分を学ぶ
https://qiita.com/gyu-don/items/efbd66c691981d7185c5

Blueqatで簡単な頂点被覆問題を解く
https://qiita.com/ttabata/items/b4ab222ef9c5ee99233c

ご意見など

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