量子コンピュータを応用したアルゴリズムの一つである**VQE(Variational Quantum Eigensolver)**について、その基本からblueqatを使用したプログラミングまでの解説を行います。
VQEとは
VQEは Variational Quantum Eigensolver の略で、変分量子固有値ソルバーと言います。ソルバーとは、数式や数値解析を用いて最適な解を求める機能またはソフトウェアのことです。
よってVQEとは、変分原理を用いて量子力学の理論から最適な固有値を求めるもの、という意味になります。
変分原理と量子力学および固有値の関係について順を追って解説します。
VQEアルゴリズムを使用することで、化学分野では材料科学計算や量子化学計算へ、また、組合せ最適化問題の分野への応用などに期待されています。
それでは、ここから変分量子固有値ソルバーという言葉を手綱に、基本からしっかり解説していきましょう。
固有値とは何か?
まず、固有値について解説します。固有値は数学の分野の話になります。
固有値の数学的な定義は以下です(難しく感じる場合はこの定義部分を一旦飛ばして読んで下さい)。
◆固有値・固有ベクトル・固有空間
_F_を_n_次元ベクトル空間_V_の線形変換、_A_を_n_次正方行列とするとき、
1. F(x)=λx かつ x≠0 なる x∈V が存在するとき、1. F(x)=λx かつ x≠0 なる x∈V が存在するとき、1. F(x)=λx かつ x=0 なる x∈V が存在するとき、
λをFの固有値、xを固有値λに対するFの固有ベクトル、λをFの固有値、xを固有値λに対するFの固有ベクトル、λをFの固有値、xを固有値λに対するFの固有ベクトル、
W(λ)=｛x∣F(x)=λx｝W(λ)=｛x|F(x)=λx｝W(λ)=｛x∣F(x)=λx｝
を固有値λに対するFの固有空間という。を固有値λに対するFの固有空間という。を固有値λに対するFの固有空間という。
2. Ax=λx かつ x≠0 なる x∈Cn が存在するとき、2. Ax=λx かつ x≠0 なる x∈C^n が存在するとき、2. Ax=λx かつ x=0 なる x∈Cn が存在するとき、
λをAの固有値、xを固有値λに対するAの固有ベクトル、λをAの固有値、xを固有値λに対するAの固有ベクトル、λをAの固有値、xを固有値λに対するAの固有ベクトル、
W(λ)=｛x∣Ax=λx｝W(λ)=｛x|Ax=λx｝W(λ)=｛x∣Ax=λx｝
を固有値λに対するAの固有空間という。を固有値λに対するAの固有空間という。を固有値λに対するAの固有空間という。
平たく説明すると、A=(ajk)をn x n行列、xxx を列ベクトルとして、以下の一次方程式
Ax=λxAx=λxAx=λx
について、ベクトルxxx に行列Aを作用させても、もとのベクトルxxx のλλλ 倍になる特別なベクトルxxx を固有ベクトル、その根となるλλλ を固有値と言います。
固有値は一般に複素数です。また、固有値は重根となる場合があり、その場合を縮退していると言います。
量子力学について次の節で詳しく説明しますが、量子力学においてはしばしば、固有ベクトルを固有関数、または固有状態と言います。
量子力学では、固有値と固有関数を以下のように記述します。
Aを物理量とするとき、Aを物理量とするとき、Aを物理量とするとき、
Aφ(x)=αφ(x) ...式(1)Aφ(x)=αφ(x) ...式(1)Aφ(x)=αφ(x) ...式(1)
ならば、αをAの固有値、φ(x)を固有関数(固有状態)という。ならば、αをAの固有値、φ(x)を固有関数(固有状態)という。ならば、αをAの固有値、φ(x)を固有関数(固有状態)という。
量子力学とVQE
量子コンピュータは原理的に量子力学の理論を応用したコンピュータと言えます。
量子力学とは、物質を構成する原子や電子、イオン化されたこれらの粒子、さらにこれらを構成するさらに小さい各種の素粒子(まとめて量子と言います)の挙動(振る舞い)を説明する理論です。
これらの量子の振る舞いを記述する最も基本的な式がシュレーディンガーの方程式です。
以下に時間依存しないシュレーディンガー方程式を紹介します。
−ℏ22md2φ(x)dx2+V(x)φ(x)=Eφ(x)-\frac{\hbar^2}{2m}\frac{d^2\varphi(x)}{dx^2}+V(x)\varphi(x) = E\varphi(x)−2mℏ2dx2d2φ(x)+V(x)φ(x)=Eφ(x)
ℏ：プランク定数\hbar：プランク定数ℏ：プランク定数
V：ポテンシャルエネルギーV：ポテンシャルエネルギーV：ポテンシャルエネルギー
E：粒子のエネルギーE：粒子のエネルギーE：粒子のエネルギー
上記の式について、−ℏ22md2dx2+V(x)-\frac{\hbar^2}{2m}\frac{d^2}{dx^2}+V(x)−2mℏ2dx2d2+V(x) をハミルトニアンHHH とおくと以下のように簡潔に記述することができます。
Hφ(x)=Eφ(x) ...式(2)H\varphi(x)=E\varphi(x) ...式(2)Hφ(x)=Eφ(x) ...式(2)
上式を見ると、前節で説明した固有値・固有ベクトルの関係である式(1)と同じ形をしていることがわかります。シュレーディンガー方程式とはハミルトニアンHHH の固有関数φ(x)\varphi(x)φ(x) 、固有値EEE を求める式であると言えます。
シュレーディンガーの方程式【式(2)】を応用して粒子の状態である固有値(エネルギー状態)を求めることがVQEの理論の裏付けとなります。
期待値の公式
変分原理についての説明に入る前に、期待値の求め方について説明します。期待値はVQEでは以下のように書かれますが、なぜそのように書くのかを確認します(この部分の説明は飛ばしても構いません)。
期待値=〈ψ∣Aψ〉 ...式(3)期待値 = 〈\psi|A\psi〉 ...式(3)期待値=〈ψ∣Aψ〉 ...式(3)
Aは物理量Aは物理量Aは物理量
一般に物理量とは、位置や運動エネルギーなどの測定可能な量のことです。物理量の測定値は実数となります。
ここで以下の量子力学の公理を確認します。
量子系ではHilbert空間(ヒルベルト空間)上で考察を行います。
物理量はHilbert空間上のエルミート演算子で表現され、測定値はその固有値とします。量子状態ψ\psiψ の下で物理量Aを測定するとき、測定値が固有値aaa となる確率PrPrPr は以下の式で与えられます。
Pr(A=a∣∣ψ〉)=〈ψ∣Paψ〉 ...式(4)Pr(A=a | |\psi〉)=〈\psi|P_a\psi〉 ...式(4)Pr(A=a∣∣ψ〉)=〈ψ∣Paψ〉 ...式(4)
PaはAの固有値aに対する固有射影演算子P_aはAの固有値aに対する固有射影演算子PaはAの固有値aに対する固有射影演算子
式(4)は確率を求める式です。この式をBornの確率規則と言います。
上式は、エルミート演算子Aに対応する物理量の測定を行った時のAの固有値のうち一つが得られる確率を意味しています。
ここでは有限次元を考えますので、Hilbert空間は一般のベクトル空間と考えて良いです。
固有値aaa は実数値です。また、PaP_aPa は固有値aaa に対する固有射影演算子です。
固有射影演算子とは、固有値を与える射影演算子で、さらに射影演算子とは以下を満たす演算子です。
Pa=Pa†=Pa2 ...式(5)P_a=P_a†=P_a^2 ...式(5)Pa=Pa†=Pa2 ...式(5)
∑aPa=I ...式(6)\sum_{a}^{}P_a=\mathbb{I} ...式(6)∑aPa=I ...式(6)
式(5)を用いて、式(4)について以下が言えます。
これは、PrPrPr が0以上を取ることを意味します。
〈ψ∣Paψ〉=〈ψ∣Pa2ψ〉=〈Paψ∣Paψ〉≧0 ...式(7)〈\psi|P_a\psi〉=〈\psi|P_a^2\psi〉=〈P_a\psi|P_a\psi〉≧0 ...式(7)〈ψ∣Paψ〉=〈ψ∣Pa2ψ〉=〈Paψ∣Paψ〉≧0 ...式(7)
同様に、式(6)を用いて式(4)について以下が言えます。
これは、各確率の和を取ると1になることを意味します。
∑a〈ψ∣Paψ〉=\sum_{a}^{}〈\psi|P_a\psi〉=∑a〈ψ∣Paψ〉=
〈ψ∣(∑aPa)ψ〉=〈ψ∣Iψ〉=∣∣ψ∣∣2=1 ...式(8)〈\psi|(\sum_{a}^{}P_a)\psi〉=〈\psi|\mathbb{I}\psi〉=||\psi||^2 = 1 ...式(8)〈ψ∣(∑aPa)ψ〉=〈ψ∣Iψ〉=∣∣ψ∣∣2=1 ...式(8)
式(7)式(8)より、式(4)は確率分布となることが分かります。
上記の公理を用いて期待値を算出することができます。
状態∣ψ〉|\psi〉∣ψ〉 のもとで物理量Aを測定するときの測定値の期待値EEE は固有値aaa と確率PrPrPr の積aPaaP_aaPa の和を取り、
Eψ[A]=∑aaPr(A=a∣∣ψ〉)E_\psi[A] = \sum_{a}^{}aPr(A=a||\psi〉)Eψ[A]=∑aaPr(A=a∣∣ψ〉)
=∑aa〈ψ∣Paψ〉=\sum_{a}^{}a〈\psi|P_a\psi〉=∑aa〈ψ∣Paψ〉
=〈ψ∣(∑aaPa)ψ〉=〈\psi|(\sum_{a}^{}aP_a)\psi〉=〈ψ∣(∑aaPa)ψ〉
=〈ψ∣Aψ〉 ...式(9)=〈\psi|A\psi〉 ...式(9)=〈ψ∣Aψ〉 ...式(9)
上式において、以下のスペクトル分解定理を使用しました。
任意の正規演算子は以下のように表すことができる。任意の正規演算子は以下のように表すことができる。任意の正規演算子は以下のように表すことができる。
A=∑a∊σ(A)aPaA=\sum_{a∊\sigma(A)}^{}aP_aA=∑a∊σ(A)aPa
よって、式(3)の期待値を求めることができました。
これを期待値の公式と呼びます。
変分原理を理解する
変分原理と言っても古典力学や電磁気学、量子力学などいくつかの分野で定義された異なる理論がありますが、ここでは量子力学での変分原理であるリッツの変分原理を取り上げます。
期待値の公式【式(3)】より、期待値を状態∣ψ〉|\psi〉∣ψ〉 と物理量Aで求めことができることがわかりました。
すでに説明した通り、物理量Aは位置やエネルギーなどさまざまな量を表す演算子です。ここで、その物理量の一つでありエネルギー量を表す演算子のHHH (ハミルトニアン)をAの代わりに用います。
ハミルトニアンHHH と量子状態∣ψ〉|\psi〉∣ψ〉 にあるときの期待値は、ある固有値(固有エネルギー)E0E_0E0 を用いて、以下のように書き表すことができます。
〈ψ∣H∣ψ〉≧E0 ...式(10)〈\psi|H|\psi〉≧E_0 ...式(10)〈ψ∣H∣ψ〉≧E0 ...式(10)
これは、量子力学の観点から以下のように説明できます。
すなわち、さまざまな状態において、電子などの粒子は量子力学の法則に従いいくつかのエネルギー状態を取ります。
最もエネルギーの低い状態を基底状態と呼びます。また、エネルギーが高くなるにつれて、順に第一励起状態、第二励起状態、さらに上位に状態が変化します(下図)。
よって、以下が成り立ちます。
適当な境界条件を持つ任意の状態ψ\psiψ に対するハミルトニアンHの期待値は、基底状態のエネルギーE0E_0E0 よりも常に大きいか等しい。
これは式(10)が成立することを意味します。これが変分原理です。
ψ\psiψ へのパラメータθ\thetaθ の導入
上記の変分原理を応用し、エネルギーが最小となる固有値E0E_0E0 を求めるためにはどうすればよいでしょうか。この答えは単純で、様々な異なる状態のを与えたときに、期待値が最も小さくなるときのψ\psiψ を探し出せば良いわけです。そのときの期待値が固有値E0E_0E0 となります。
様々な異なるψ\psiψ を与えるために、新たに角度のパラメータθ\thetaθ を用意します。そして、ψ\psiψ がパラメータθ\thetaθ に依存する一種の関数とみなします。
期待値は以下のように書き直すことができます。
〈ψ(θ)∣H∣ψ(θ)〉〈\psi(\theta)|H|\psi(\theta)〉〈ψ(θ)∣H∣ψ(θ)〉
VQEアルゴリズムを用いたプログラミングの実践
変分原理を理解した上で、VQEを用いたプログラミングに入りたいと思います。
プログラミングには、量子コンピューターシミュレータとして評価の高い、blueqat株式会社(東京・本郷)が提供するblueqat(ブルーキャット)を使用します。
また、解く問題としては組合せ最適化問題の一つである頂点被覆問題を取り扱います。
(1)頂点被覆問題について
頂点被覆問題(Vertex Cover)とは、グラフ・ネットワーク問題の一つで、複数の頂点とそれを結ぶ辺があり、どの辺も少なくともその始点または終点のどちらかが、二つにクラス分けした頂点群のうちの一方のクラスに属する頂点(のいずれか)と接するもので、そのクラスを最小にする問題です。
数学的に説明すると以下となります。
頂点の部分集合X⊆Vに対して、どの枝も少なくとも一方の端点がXに含まれるとき、すなわち、任意の(i,j)∈Eに対してX∩i,j≠0が成り立つとき、Xを頂点被覆という。頂点の部分集合X⊆Vに対して、どの枝も少なくとも一方の端点がXに含まれるとき、すなわち、任意の(i,j)∈Eに対してX∩i,j≠0が成り立つとき、Xを頂点被覆という。頂点の部分集合X⊆Vに対して、どの枝も少なくとも一方の端点がXに含まれるとき、すなわち、任意の(i,j)∈Eに対してX∩i,j=0が成り立つとき、Xを頂点被覆という。
以下に図を示します。
ここで、v1,v2,v3,v4,v5は頂点を、e1,e2,e3,e4,e5は辺を表します。図1の通り、どの辺についてもその端点の一方が赤丸の頂点に属しています。赤丸の数が最小になるケースを求める問題が頂点被覆問題です。
例えば、下の図2のような頂点を選択した場合、頂点の数は3個なので最適解にはなりません。
この例題を解く際には、イジング模型に基づくハミルトニアンの定式化を行います。今回の例題では、ハミルトニアンHHH を以下のように定義します。
H=5−32x1−32x2−32x3−52x4−12x5+x1x2+x1x4+x2x3+x3x4+x4x5H=5-\frac{3}{2}x_1-\frac{3}{2}x_2-\frac{3}{2}x_3-\frac{5}{2}x_4-\frac{1}{2}x_5+x_1x_2+x_1x_4+x_2x_3+x_3x_4+x_4x_5H=5−23x1−23x2−23x3−25x4−21x5+x1x2+x1x4+x2x3+x3x4+x4x5
上記の定式化についての詳細は、以下の技術ブログをご参照下さい。
■Blueqatで簡単な頂点被覆問題を解く
https://qiita.com/ttabata/items/b4ab222ef9c5ee99233c
(2)変分原理によるVQEアルゴリズムの計算手順
変分原理を用いたVQEアルゴリズムは以下の手順で計算を行います。
①パラメータθ\thetaθ を用いた量子回路を準備(ansatz：アンザッツ)
②量子ビットを測定し、期待値を算出する
③古典最適化計算により、期待値を最小にするパラメータθ\thetaθ を探索する
VQEでは、量子コンピュータによる期待値の測定と、古典最適化計算(古典コンピュータの計算)によるパラメータθ\thetaθ の探索処理の両方を用いて解を算出します。
処理①~③を模式図で示すと以下となります。こちらはVQEの論文からの引用です。
引用：(https://arxiv.org/pdf/1304.3061.pdf)
処理①に記載のansatzとは、量子計算に必要な量子状態∣ψ〉|\psi〉∣ψ〉 を生成するための量子回路のことで、任意回転ゲートと入力パラメータθ\thetaθ を用いて量子状態を生成します。この量子状態を用いて期待値を計算します。
(3)blueqatによるVQEアルゴリズム計算のプログラミング
blueqatには上記の計算手順①~③を内部関数で実行してくれる機能が備わっており、プログラマーはパラメータや期待値計算、最小値の探索に手をかけなくても良くなっています。これはblueqatが優れている点の一つと言えます。
プログラマーがすべきことは、パラメータとハミルトニアンHHH を定義することのみです。
では、さっそくパラメータ及びハミルトニアンHHH の定義を行いましょう。前掲の例題である頂点被覆問題から、ハミルトニアンを以下とします。
H=5−32q0−32q1−32q2−52q3−12q4+q0q1+q0q3+q1q2+q2q3+q3q4H=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_4H=5−23q0−23q1−23q2−25q3−21q4+q0q1+q0q3+q1q2+q2q3+q3q4
例題では変数をxxx としていますが、ここでは量子ビット(quantum bit)の意味で、アルファベットqqq を使用しました。また、インデックス番号を0始まりとしています。
基本的にはこのハミルトン式を組み込めば良いです。プログラムのソースコードは以下となります。
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個分を宣言します。
「ハミルトニアンの定義」にハミルトニアンを定義します。
プログラミングは以上です。
(4)実行結果の確認
実行結果は以下となりました。結果の見方は例題と同じです。
q0～q4は以下の値を取ります。
ここでは、0,1はクラス分けを意味しています。0を取る頂点のクラスと、1を取る頂点のクラスに分かれます。
このクラス分けが今回の問題の解答となります。
つまり、q1とq3が同じクラス、q0,q2,q4が別の同じクラスに属することがわかります。
q1,q3は図1の赤丸の頂点、q0,q2,q4はその他の丸の頂点に分けられます。
この分け方が問題の解となります。
以上でプログラミングの説明を終わります。