Help us understand the problem. What is going on with this article?

VQEセミナー2周目まとめ

はじめに

量子コンピュータの勉強会をしていますが、VQEと言うアルゴリズムを勉強して1周終わりました。2周目のためにまとめておきます。

VQE#2-1(理論編)

主に加藤くんがまとめてくれたものを勉強会向けにまとめなおしてみます。より詳しく学びたい方は、こちらをみてみてください。

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

2-1-1 VQEとは?

VQEはVariational Quantum Eigensolverと言うアルゴリズムです。現在の量子コンピュータは理想的な量子コンピュータに対してエラーがかなりあるので長い回路の実行が厳しいです。従来期待されていた位相推定アルゴリズムというアルゴリズムに対して、量子古典ハイブリッドアルゴリズムとして現在の量子コンピュータでも利用できるものとして2013年に開発されました。

量子回路は短いものを実行し、それを複数実行します。その結果を古典コンピュータで集計し、最適化をかけます。その繰り返しを行うことで、短い回路で特定の機能を実現できます。

2-1-2 VQEの目的

VQEはエルミート行列の固有値(の期待値)を求めるためのアルゴリズムです。これがもとまると様々な問題に応用することができ、現在でも量子化学計算や組合せ最適化問題を解くことができます。

2-1-3 VQEの計算手順

手順は簡単です。

1、パラメータ化された量子回路を準備する
2、量子ビットの測定値からハミルトニアンの期待値を計算する
3、パラメータを探索して上記の期待値を最小化する

ハミルトニアンは固有値を求める行列に対応します。

2-1-4 量子変分原理

量子変分原理より、任意の波動関数におけるハミルトニアンの期待値は基底状態を上回るというのがあります。

\langle \Psi (\theta) \mid H \mid \Psi (\theta) \rangle \geq E_0

ハミルトニアンの値を最小にするパラメータ$\theta$を探すのが今回の目的です。

2-1-5 ansatz

実際に利用されるのはansatz(アンザッツ)と呼ばれる量子状態を作るための回路で、角度でパラメータ化されています。任意の量子状態を回転ゲートを使って表現でき、良い量子状態が見つかってハミルトニアンの期待値を求めれば基底状態を探すことができます。

探索するための量子回路は任意の量子回路をかければいいので、下記のようにblueqatを使う場合には、

Circuit().ry(a)[0].rz(b)[0]

のようにかけたりします。

from strangeworks.blueqat import StrangeworksBackend #quantumcomputing.comの時だけ必要
import numpy as np
from blueqat import Circuit
from blueqat.pauli import X, Y, Z, I
from blueqat.vqe import AnsatzBase, Vqe

class OneQubitAnsatz(AnsatzBase):
    def __init__(self, hamiltonian):
        super().__init__(hamiltonian.to_expr(), 2)
        self.step = 1

    def get_circuit(self, params):
        a, b = params
        return Circuit().ry(a)[0].rz(b)[0]


# この定式化が大事
h = Z(0)
runner = Vqe(OneQubitAnsatz(h))
result = runner.run()

print('Result by VQE')
print(runner.ansatz.get_energy(result.circuit, runner.sampler))

2-1-6 1量子ビットのansatzを使ってハミルトニアンZの期待値を求める

上記ansatzを使って任意の量子状態を作れたとして、量子状態は、

\mid \psi \rangle = \alpha \mid 0 \rangle + \beta \mid 1 \rangle

で表されます。例えば求めたいハミルトニアンが、Zの場合、

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

求めたい式は、

\langle \psi \mid H \mid \psi \rangle = 

\begin{bmatrix}
\alpha&\beta
\end{bmatrix}

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

\begin{bmatrix}
\alpha\\
\beta
\end{bmatrix}

= |\alpha|^2 - |\beta|^2

となります。ここで、$\alpha^2$は0がでる確率、$\beta^2$は1が出る確率に対応しますので、最小になるのは、$\beta$=1となるときで、-1になります。

上記ansatzでは、Ryに180度回転が適用された時に最小となり、計算結果は最適化された状態$\mid 1 \rangle$の確率振幅が1に近づきます。

2-1-7 Xも

ハミルトニアンがXの時には、

\langle \psi \mid H \mid \psi \rangle = 

\begin{bmatrix}
\alpha&\beta
\end{bmatrix}

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

\begin{bmatrix}
\alpha\\
\beta
\end{bmatrix}

= 2 \alpha \beta

となりますが、通常の測定ではこの値を求めることができません。そこで測定の基底を変更します。

\mid \psi \rangle = \frac{\alpha + \beta}{\sqrt{2}}\mid + \rangle + \frac{\alpha - \beta}{\sqrt{2}}\mid - \rangle

+と-状態というのは、通常の状態から軸を90度回転させてX軸上で取られた状態です。

+状態を測定すると、確率振幅から、出現確率が、

\frac{\alpha^2 + 2\alpha\beta + \beta^2}{2}

となります。一方-状態は、

\frac{\alpha^2 - 2\alpha\beta + \beta^2}{2}

となりますので、前者から後者を引くと、

2\alpha\beta

となります。これが求めたかったハミルトニアンの期待値になります。通常、ハミルトニアンがXであると分かっている場合には、測定する前にアダマールゲートをかけてX軸の基底にしてから測定をします。

2-1-8 ハミルトニアンの変形

ハミルトニアンの要素がZの場合には、そのまま計算ができますが、XやYの場合には上記の基底を変換するという考え方がわかりづらい場合、

X = HZH

より、

\langle \psi \mid X \mid \psi \rangle \\
= \langle \psi \mid HZH \mid \psi \rangle\\
= \langle \psi' \mid Z \mid \psi' \rangle

のように波動関数を変形することでも対応できます。

Y = RX(-\pi/2) Z RX(\pi/2)

なので、これも同様に、

\langle \psi \mid Y \mid \psi \rangle \\
= \langle \psi \mid RX(-\pi/2) Z RX(\pi/2) \mid \psi \rangle\\
= \langle \psi'' \mid Z \mid \psi'' \rangle

とすることができます。

2-1-9 線形結合?

ハミルトニアンは下記のような場合、分解できます。

\langle \psi \mid aH_1 + bH_2 \mid \psi \rangle \\ = \langle \psi \mid aH_1 \mid \psi \rangle + \langle \psi \mid bH_2 \mid \psi \rangle \\ = a\langle \psi \mid H_1 \mid \psi \rangle + b\langle \psi \mid H_2 \mid \psi \rangle

これを使えば、複雑なハミルトニアンの場合、

H = 1.2 X_0 Z_2 + 2.5 Z_3 X_4 Y_2 - 3.4 Z_2 X_3

期待値は、

\langle \psi \mid 1.2 X_0 Z_2 + 2.5 Z_3 X_4 Y_2 - 3.4 Z_2 X_3 \mid \psi \rangle\\
= 1.2\langle \psi \mid X_0 Z_2 \mid \psi \rangle + 2.5 \langle \psi \mid Z_3 X_4 Y_2\mid \psi \rangle - 3.4 \langle \psi \mid Z_2 X_3 \mid \psi \rangle

のように個別にできます。

2-1-10 あとは計算

すべてblueqatで自動でやってくれます。pythonからpip install blueqatです。

from blueqat import vqe
from blueqat.pauli import Z

H = Z(0)*Z(1)

step = 2

result = vqe.Vqe(vqe.QaoaAnsatz(H, step)).run()
print(result.most_common(12))

VQE#2-2(実装初級編)

今回は前回学んだ内容を元に、実際にプログラムに落としてみたいと思います。まずはソフトウェアの導入からです。もちろん全て無料です。

2-2-1 インストール

今回はblueqatsdkを使います。

https://github.com/Blueqat/Blueqat

手元にpythonとpipのインストールが完了している、もしくはgoogle colabを使う場合には、

https://colab.research.google.com/notebook

インストールは簡単で、

pip3 install blueqat

で終わります。また、quantumcomputing.comというサービスがあります。こちらはいちばん簡単ですが、登録が必要です。

https://app.quantumcomputing.com/

にアクセスすればいいですが、アカウントの作成は、

https://app.quantumcomputing.com/signup

から。現在はβ運用中で今年製品版が出ます。無料です。こちらは最初からblueqatが入っているのでインストール不要です。

2-2-2 コード

早速コードを見てみましょう。下記のコードを入れてみます。大事なところは2カ所です。今回は1量子ビットの行列の固有値を求めてみます。

from strangeworks.blueqat import StrangeworksBackend #quantumcomputing.comのみ必要
import numpy as np
from blueqat import Circuit
from blueqat.pauli import X, Y, Z, I
from blueqat.vqe import AnsatzBase, Vqe

class OneQubitAnsatz(AnsatzBase):
    def __init__(self, hamiltonian):
        super().__init__(hamiltonian.to_expr(), 2)
        self.step = 1

    def get_circuit(self, params):
        a, b = params
        return Circuit().ry(a)[0].rz(b)[0]


# Calculate by VQE
cost = 1.23 * I - 4.56 * X(0) + 2.45 * Y(0) + 2.34 * Z(0)
runner = Vqe(OneQubitAnsatz(cost))
result = runner.run()

print('Result by VQE')
print(runner.ansatz.get_energy(result.circuit, runner.sampler))

# Hamiltonian to matrix
mat = cost.to_matrix()

# Calculate by numpy
print('Result by numpy')
print(np.linalg.eigh(mat)[0][0])

上記は、costで決められた行列の期待値を求め、その答え合わせをnumpyというよく使うツールで確かめをしています。

上記のポイントは2カ所。Circuit()....のところと、一番大事なのは、costと書いてあるところです。

2-2-3 VQEの定式化手順

VQEのアルゴリズムは通常簡単に使うことができます。

1、定式化
2、計算

これだけです。定式化は量子ゲートと呼ばれるものを使って計算をします。blueqatの場合、QAOAと呼ばれるアルゴリズムもとても簡単に利用できます。QAOAの方がVQEよりも簡単に使えるように実装されています。今回はVQEを使います。

from blueqat.pauli import X, Y, Z, I

h = 1.23 * I - 4.56 * X(0) + 2.45 * Y(0) + 2.34 * Z(0)

定式化は、パウリ行列と呼ばれる基本の量子ゲートと単位行列を使って実現できます。上記は、

「1量子ビットVQEを世界一簡単に書く 2020年2月Ver.」
https://qiita.com/gyu-don/items/65f830ae3777534f9a04

の例題からですが、IとXとYとZというゲートを使って2*2の行列を作っています。XやYのあとの0という数字は何番目の量子ビットを使うかという通し番号になっています。

量子ビットの通し番号は通常0からスタートします。

ちなみに上記の行列は、

[[ 3.57,   -4.56-2.45j]
 [-4.56+2.45j, -1.11  ]]

となっています。

2-2-4 Ansatzの準備

ansatzの準備も簡単です。

class OneQubitAnsatz(AnsatzBase):
    def __init__(self, hamiltonian):
        super().__init__(hamiltonian.to_expr(), 2)
        self.step = 1

    def get_circuit(self, params):
        a, b = params
        return Circuit().ry(a)[0].rz(b)[0]

この部分で、Circuit().ry(a)[0].rz(b)[0]がansatzです。今回は主にこれを計算します。

2-2-5 結果

結果は出力され、確かめもできます。

Result by VQE
-4.450591481563331
Result by numpy
-4.450818602983201

期待値なので、ほぼ近い値が出ています。

2-2-6 quantumcomputing.com

こちらの共有リンクも貼っておきます。

https://app.quantumcomputing.com/public/projects/udFDCw

こちらを実行してもらうと、

式の途中にパラメータが使われていて、上記で設定したansatzのパラメータが最適化されて答えが出ています。

Circuit().ry(a)[0].rz(b)[0]

中身を理解するのはちょっと時間がかかるので、とりあえず

1、行列を作る
2、実行して計算をする

を中心に進めました。

VQE#2-3(組合せ最適化イジング+QUBO編)

量子コンピュータは正直どっから手をつけていいかわからんと思います(もしくは思う人もいると思います)。実際1から量子計算を覚えたりしてプログラムをして仕事に生かしたいという人もいますし、手っ取り早くプログラムを覚えたいという人も。

今回は、

1、まずプログラムを覚える
2、壁に当たったら理論に戻る

という順でやることを想定します。疑問があったらslackでいつでも聞けるので。

https://join.slack.com/t/blueqat/shared_invite/enQtNzc3NjgxNjA5OTA2LTNmNjg0YzAzNGUxYmYwNDIxNWFlZGVkODI0MTE2NDMwM2I3OTRjNGExZjMxZWFmZDI3ZTRlNDU3NTcwZDUyOTk

2-3-1 組合せ最適化問題

組合せ最適化問題は多くの選択肢からベストな答えを選ぶ問題です。社会問題でよく例題として利用されます。

VQEで解ける問題の一つに、組合せ最適化問題があります。社会問題をうまく式の形にし、その式をVQEで解くことで答えを取ることができます。

主に選択肢を選んだ状態を1、選ばない状態を0としますが、様々な条件をつけながら問題を解きます。

今回は社会問題を量子コンピュータに応用する手筈としてVQEに注目してみたいと思います。

2-3-2 定式化

今回は定式化と言って式を作ってみます。ルールは、

・問題をイジングモデルと呼ばれるモデルに落とし込む
・Zだけを使う

ということです。

一応例題のリンクはこちらです。

https://app.quantumcomputing.com/public/projects/udFDCw

今回は、

cost = -Z(0) - Z(0)*Z(1)

という問題をやってみます。まずはこれを解きます。

from strangeworks.blueqat import StrangeworksBackend #quantumcomputing.comの時だけ必要
import numpy as np
from blueqat import Circuit
from blueqat.pauli import X, Y, Z, I
from blueqat.vqe import AnsatzBase, Vqe

class OneQubitAnsatz(AnsatzBase):
    def __init__(self, hamiltonian):
        super().__init__(hamiltonian.to_expr(), 4)
        self.step = 1

    def get_circuit(self, params):
        a, b, c, d = params
        return Circuit().ry(a)[0].rz(b)[0].ry(c)[1].rz(d)[1]


# この定式化が大事
cost = -Z(0) - Z(0)*Z(1)
runner = Vqe(OneQubitAnsatz(cost))
result = runner.run()

print('Result by VQE')
print(runner.ansatz.get_energy(result.circuit, runner.sampler))

そうすると、

Result by VQE
-1.9999943170511056

約-2が答えとして出てきました。これで終わりです。ルールを説明します。

2-3-3 Zは1か-1をとる

h = -Z(0) - Z(0)*Z(1)

Zの後ろの数字は、量子ビットの通し番号を表します。0番目と1番目の量子ビットの2つを使っています。また、問題設定で大事なのは、Zの前の係数です。

Z(0)の前は-1
Z(0)*Z(1)の前は-1

がついています。業界では、それぞれ「局所磁場」と「相互作用」とか会社によって呼び名が違いますが、ここではあえて、「バイアス」と「ウェイト」という名前にしてみます。

Zは-1か+1のどちらかをとります。
また、hはより小さい値を取ると正解になります。
上記VQEではhの一番小さい数字を取るように計算されます。

手計算で考えてみると、
Z(0)のバイアスが-1で、Z(0)*Z(1)のウェイトが-1の時、
Z(0)=+1で、
Z(1)=+1の時に

h = -1*1 - 1*1 = -2

で最小になるはずです。これを自動でやってくれます。

Z(0),Z(1)の値で場合わけすると、

Z(0) Z(1) h
-1 -1 0
-1 1 2
1 -1 0
1 1 -2

このZ(0)とZ(1)を探すのが、今回の目的なので、達成されました。

2-3-4 この問題を大きくするだけで社会問題がたくさん解ける

そして実際にはより大きな問題を解きます。

上記問題は役にたたなさそうですが、Z(0)をAさん、Z(1)をBさんに見立てて、Aさんはグループ1に属し、BさんはAさんと同じグループに属するという条件をつけた分類問題と同じです。

ただ、このままでは毎回問題を解くのが大変なので様々な工夫がされています。-1と+1の代わりに0と1を使うなど覚えるべきテクニックがたくさんあります。

2-3-5 定式化は01で

定式化とは、上記の記事でもちょっと触れていますが、組合せ最適化問題を+1と-1の値の組合せで表される式に落とし込みます。ただ、課題があります。物理学で使われるイジングモデルは-1,+1を利用しますが、通常の産業用では0と1を計算として利用します(計算基底)。幸い-1と+1は自動的に変換ができますので、01を使っても組合せ最適化問題としては問題がありません。それがこの8年ほど発達してきました。

01で定式化をするのをQUBOと呼びます。QUBOは量子ビットの01を利用できます。係数はまだ名前が確定していませんが、「バイアス」と「ウェイト(結合荷重)」でいいでしょう。定式化で作るのは「コスト関数」と呼びたいと思います。

定式化は社会問題をQUBO形式で、コスト関数を作ることで実現でき、コスト関数は01の値をとる量子ビットに設定するバイアスとウェイトを設定することで実現できます。

この手法は量子コンピュータに限ったことではないので普通の計算に慣れている人でも受け入れられやすいでしょう。

2-3-6 早速プログラミング

今回は前回から少し進んで、

from strangeworks.blueqat import StrangeworksBackend #quantumcomputing.comを利用する時に必要
import numpy as np
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, 4)
        self.step = 1

    def get_circuit(self, params):
        a, b, c, d = params
        return Circuit().ry(a)[0].rz(b)[0].ry(c)[1].rz(d)[1]

cost = -3*q(0)-3*q(1)-2*q(0)*q(1)
cost = cost.to_expr().simplify()
runner = Vqe(QubitAnsatz(cost))
result = runner.run()

print('Result by VQE')
print(runner.ansatz.get_energy(result.circuit, runner.sampler))

# Hamiltonian to matrix
mat = cost.to_matrix()

# Calculate by numpy
print('Result by numpy')
print(np.linalg.eigh(mat)[0][0])

print('Matrix is')
print(mat)

quantumcomputing.comのリンクは、
https://app.quantumcomputing.com/public/projects/gxUVVo

今回のと期待コスト関数は、

cost = -3*q(0)-3*q(1)-2*q(0)*q(1)

明らかにq(0)=1,q(1)=1の時に最小値-3-3-2=-8を取りますね。上記のコスト関数は以前のように-1と+1ではなく、0と1で考えることができるので簡単です。

設定するのは、q(0)とq(1)に-3のバイアスを設定し、q(0)*q(1)に-2のウェイトを設定します。

これをVQEで解くことで、最小値-8が得られます。

ポイントは、

class QubitAnsatz(AnsatzBase):
    def __init__(self, hamiltonian):
        super().__init__(hamiltonian, 4)
        self.step = 1

    def get_circuit(self, params):
        a, b, c, d = params
        return Circuit().ry(a)[0].rz(b)[0].ry(c)[1].rz(d)[1]

となっていて、パラメータをa,b,c,dの4つを用意し、回路も2量子ビットに対応しています。この状態ですと、最適化するパラメータはN量子ビットに対して、2N必要になりそうです。すいません、本当はエンタングルメントを考えるべきなのでansatzはもうちょい複雑な式になるはずです。いい例題を今度探しますが、今回は勘弁してください。

Result by VQE
-7.981654210576551
Result by numpy
-8.0
Matrix is
[[ 0.+0.j  0.+0.j  0.+0.j  0.+0.j]
 [ 0.+0.j -3.+0.j  0.+0.j  0.+0.j]
 [ 0.+0.j  0.+0.j -3.+0.j  0.+0.j]
 [ 0.+0.j  0.+0.j  0.+0.j -8.+0.j]]

あとは、定式化の方法を工夫して覚えることになります。

2-3-7 QAOA

上記のVQEは行列の固有値を求めるだけです。組合せ最適化問題では定式化にpauli-Zしか使わないので、QAOAというアルゴリズムでも解くことができます。

class QaoaAnsatz(AnsatzBase):
    def __init__(self, hamiltonian, step=1, init_circuit=None):
        super().__init__(hamiltonian, step * 2)
        self.hamiltonian = hamiltonian.to_expr().simplify()
        if not self.check_hamiltonian():
            raise ValueError("Hamiltonian terms are not commutable")

        self.step = step
        self.n_qubits = self.hamiltonian.max_n() + 1
        if init_circuit:
            self.init_circuit = init_circuit
            if init_circuit.n_qubits > self.n_qubits:
                self.n_qubits = init_circuit.n_qubits
        else:
            self.init_circuit = Circuit(self.n_qubits).h[:]
        self.init_circuit.make_cache()
        self.time_evolutions = [term.get_time_evolution() for term in self.hamiltonian]

    def check_hamiltonian(self):
        """Check hamiltonian is commutable. This condition is required for QaoaAnsatz"""
        return self.hamiltonian.is_all_terms_commutable()

    def get_circuit(self, params):
        c = self.init_circuit.copy()
        betas = params[:self.step]
        gammas = params[self.step:]
        for beta, gamma in zip(betas, gammas):
            beta *= np.pi
            gamma *= 2 * np.pi
            for evo in self.time_evolutions:
                evo(c, gamma)
            c.rx(beta)[:]
        return c

中身のコードを見てみると、時間発展+変分パラメータ設定を組み合わせており、ちょっと複雑なことをしているのが分かります。QAOAを詳しく知りたい方は是非その辺りのキーワードで検索をしてみてください。

最後にほぼSDK側で処理されたQAOAを使った計算をしてみます。
設定するのはcostです。

from strangeworks.blueqat import StrangeworksBackend #quantumcomputing.comの時だけ必要
from blueqat import vqe
from blueqat.pauli import qubo_bit as q

cost = -3*q(0)-3*q(1)-2*q(0)*q(1)
step = 2

result = vqe.Vqe(vqe.QaoaAnsatz(cost, step)).run()
print(result.most_common(12))

回路はちょっと複雑になりました。

答えは、

(((1, 1), 0.9922852962030535), ((0, 0), 0.0047331664261775755), ((1, 0), 0.0014907686853841205), ((0, 1), 0.0014907686853841168))

https://app.quantumcomputing.com/public/projects/XMtzWE

最初の数字が答えの候補で、次の数字はそれが現れる確率です。11が最大確率になっているのでOKでした。

2-3-8 組合せ最適化問題は難しくない

ということで、01の量子ビットでの定式化をして、あとはツールに任せればいいというのが分かりました。最初はツールに任せて定式化を頑張りましょう!

2-3-8 実例:交通最適化

交通最適化は度々出てくるテーマです。

youtubeでは、
交通最適化
https://www.youtube.com/watch?v=JU2omOzuOQg

がわかりやすいです。記事は、

「D-WaveでVW社の交通最適化アプリケーションの実装を解く」
https://qiita.com/YuichiroMinato/items/19dd89e381f14ad24134

「交通最適化問題を量子アニーリングで解く」
https://qiita.com/sychocola1/items/646d2c9504e34ed1589f

今回は、交通最適をときます。やることは、

1、1台の車に3つのルート選択を許容する(古典)
2、現在の選択ルートから混雑状況を計算(古典)
3、混雑度を最小にするようにルート選択を最小化(QAOA)

となります。問題は使い回しをして、上記のリンクから。

「Quantum Computing at Volkswagen:
Traffic Flow Optimization using the D-Wave Quantum Annealer」
引用:https://www.dwavesys.com/sites/default/files/VW.pdf

一つ目は混雑度をコスト関数として、

(q_0+q_1+q_3+q_4)^2+(q_0+q_1+q_3+q_4)^2+(q_0+q_3)^2+(q_0+q_3)^2+(q_1+q_4)^2+(q_1+q_2+q_4+q_5)^2+(q_2+q_5)^2+(q_2+q_5)^2+(q_2+q_5)^2\\ 
=2(q_0+q_1+q_3+q_4)^2+2(q_0+q_3)^2+(q_1+q_4)^2+(q_1+q_2+q_4+q_5)^2+3(q_2+q_5)^2\\ 
=4q_0^2+4q_0q_1+8q_0q_3+4q_0q_4+4q_1^2+2q_1q_2+4q_1q_3+8q_1q_4+2q_1q_5+4q_2^2 +2q_2q_4+8q_2q_5+4q_3^2+4q_3q_4+4q_4^2+2q_4q_5+4q_5^2\\
=4q_0+4q_0q_1+8q_0q_3+4q_0q_4+4q_1+2q_1q_2+4q_1q_3+8q_1q_4+2q_1q_5+4q_2 +2q_2q_4+8q_2q_5+4q_3+4q_3q_4+4q_4+2q_4q_5+4q_5

3つのルートのうちに1つをえらぶ制約条件をこちらに、

K(q_0+q_1+q_2-1)^2+K(q_3+q_4+q_5-1)^2 = -Kq_0+2Kq_0q_1+2Kq_0q_2-Kq_1+2Kq_1q_2-Kq_2-Kq_3+2Kq_3q_4+2Kq_3q_5-Kq_4+2Kq_4q_5-Kq_5+2K
from blueqat import vqe
from blueqat.pauli import qubo_bit as q

K=2
cost2 = -K*q(0)+2*K*q(0)*q(1)+2*K*q(0)*q(2)-K*q(1)+2*K*q(1)*q(2)-K*q(2)-K*q(3)+2*K*q(3)*q(4)+2*K*q(3)*q(5)-K*q(4)+2*K*q(4)*q(5)-K*q(5)

step = 4

result = vqe.Vqe(vqe.QaoaAnsatz(cost2, step)).run()
print(result.most_common(12))

この二つを繋げてQUBO式とします。

そして、実行します。

from strangeworks.blueqat import StrangeworksBackend #quantumcomputing.comを利用する時に必要

from blueqat import vqe
from blueqat.pauli import qubo_bit as q

K=10
cost1 = 4*q(0)+4*q(0)*q(1)+8*q(0)*q(3)+4*q(0)*q(4)+4*q(1)+2*q(1)*q(2)+4*q(1)*q(3)+8*q(1)*q(4)+2*q(1)*q(5)+4*q(2)+2*q(2)*q(4)+8*q(2)*q(5)+4*q(3)+4*q(3)*q(4)+4*q(4)+2*q(4)*q(5)+4*q(5)
cost2 = -K*q(0)+2*K*q(0)*q(1)+2*K*q(0)*q(2)-K*q(1)+2*K*q(1)*q(2)-K*q(2)-K*q(3)+2*K*q(3)*q(4)+2*K*q(3)*q(5)-K*q(4)+2*K*q(4)*q(5)-K*q(5)

step = 8

result = vqe.Vqe(vqe.QaoaAnsatz(cost1+cost2, step)).run()
print(result.most_common(12))

なかなか答えが出ませんが、、、、

(((0, 1, 0, 0, 0, 1), 0.12986970172628198), ((0, 0, 1, 0, 1, 0), 0.1298697017262819), ((0, 0, 0, 0, 0, 0), 0.08232027492914641), ((0, 1, 0, 0, 1, 0), 0.05463053577530196), ((0, 0, 0, 1, 0, 0), 0.04609411831952029), ((1, 0, 0, 0, 0, 0), 0.04609411831952023), ((0, 0, 1, 1, 0, 1), 0.037740176136034545), ((1, 0, 1, 0, 0, 1), 0.0377401761360345), ((1, 0, 0, 0, 1, 1), 0.03131986281902896), ((0, 1, 1, 1, 0, 0), 0.031319862819028946), ((0, 0, 1, 0, 0, 1), 0.03071276141604334), ((1, 0, 0, 1, 0, 0), 0.02140326399538612))

うまくいけば制約条件を満たす解を与えます。

2-4 VQEセミナー#2-4(量子化学計算H2編)

前回の一連のセミナーには量子化学計算は入っていませんでしたが、入れてみましょう。やるべきことはこれまでの文脈からすると、ハミルトニアン・コスト関数の作成をパウリ行列を使って行います。

2-4-1 追加のインストール

追加でツールのインストールが必要です。
今回は、blueqatに追加で、openfermionblueqatとopenfermionを追加します。

pip install blueqat openfermionblueqat openfermion

準備OKです。

2-4-2 コードをざっくり

今回はopenfermionに準備されたコードを使います。

from blueqat import *
from openfermion import *
from openfermionblueqat import*
import numpy as np

def get_molecule(bond_len):
  geometry = [('H',(0.,0.,0.)),('H',(0.,0.,bond_len))]

  description = format(bond_len)
  molecule = MolecularData(geometry, "sto-3g",1,description=description)

  molecule.load()
  return molecule

x = [];e=[];fullci=[]
for bond_len in np.arange(0.2,2.5,0.1):
  m = get_molecule("{:.2}".format(bond_len))
  h = bravyi_kitaev(get_fermion_operator(m.get_molecular_hamiltonian()))
  runner = vqe.Vqe(UCCAnsatz(h,6,Circuit().x[0]))
  result = runner.run()
  x.append(bond_len)
  e.append(runner.ansatz.get_energy(result.circuit,runner.sampler))
  fullci.append(m.fci_energy)

%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(x,fullci)
plt.plot(x,e,"o")

こちらが正しく実行されると下記のような画像が出力されます。

ダウンロード.png

今回は水素分子の結合間の距離に応じて、各状態の安定エネルギーをVQEを用いて求めています。

2-4-3 中身確認

MoluculearDataというのでOpenFermionに最初から入っているデータを使えるようです。H2でSto-3g系を読み込んでみます。

def get_molecule(bond_len):
  geometry = [('H',(0.,0.,0.)),('H',(0.,0.,bond_len))]

  description = format(bond_len)
  molecule = MolecularData(geometry, "sto-3g",1,description=description)

  molecule.load()
  return molecule

チェックしてみると、原子間距離が0.4のものを読み込んで、

m = get_molecule(0.4)

そのあと、

m.get_molecular_hamiltonian()

これによってハミルトニアンと呼ばれる基本の式が取得できます。

() 1.322943021475
((0, 1), (0, 0)) -1.4820918858979102
((1, 1), (1, 0)) -1.4820918858979102
((2, 1), (2, 0)) -0.1187350527865787
((3, 1), (3, 0)) -0.1187350527865787
((0, 1), (0, 1), (0, 0), (0, 0)) 0.36843967630348756
((0, 1), (0, 1), (2, 0), (2, 0)) 0.08225771204699692
((0, 1), (1, 1), (1, 0), (0, 0)) 0.36843967630348756
((0, 1), (1, 1), (3, 0), (2, 0)) 0.08225771204699692
((0, 1), (2, 1), (0, 0), (2, 0)) 0.082257712046997
((0, 1), (2, 1), (2, 0), (0, 0)) 0.3626667179796745
((0, 1), (3, 1), (1, 0), (2, 0)) 0.082257712046997
((0, 1), (3, 1), (3, 0), (0, 0)) 0.3626667179796745
((1, 1), (0, 1), (0, 0), (1, 0)) 0.36843967630348756
((1, 1), (0, 1), (2, 0), (3, 0)) 0.08225771204699692
((1, 1), (1, 1), (1, 0), (1, 0)) 0.36843967630348756
((1, 1), (1, 1), (3, 0), (3, 0)) 0.08225771204699692
((1, 1), (2, 1), (0, 0), (3, 0)) 0.082257712046997
((1, 1), (2, 1), (2, 0), (1, 0)) 0.3626667179796745
((1, 1), (3, 1), (1, 0), (3, 0)) 0.082257712046997
((1, 1), (3, 1), (3, 0), (1, 0)) 0.3626667179796745
((2, 1), (0, 1), (0, 0), (2, 0)) 0.36266671797967454
((2, 1), (0, 1), (2, 0), (0, 0)) 0.08225771204699726
((2, 1), (1, 1), (1, 0), (2, 0)) 0.36266671797967454
((2, 1), (1, 1), (3, 0), (0, 0)) 0.08225771204699726
((2, 1), (2, 1), (0, 0), (0, 0)) 0.08225771204699728
((2, 1), (2, 1), (2, 0), (2, 0)) 0.38272169831413727
((2, 1), (3, 1), (1, 0), (0, 0)) 0.08225771204699728
((2, 1), (3, 1), (3, 0), (2, 0)) 0.38272169831413727
((3, 1), (0, 1), (0, 0), (3, 0)) 0.36266671797967454
((3, 1), (0, 1), (2, 0), (1, 0)) 0.08225771204699726
((3, 1), (1, 1), (1, 0), (3, 0)) 0.36266671797967454
((3, 1), (1, 1), (3, 0), (1, 0)) 0.08225771204699726
((3, 1), (2, 1), (0, 0), (1, 0)) 0.08225771204699728
((3, 1), (2, 1), (2, 0), (3, 0)) 0.38272169831413727
((3, 1), (3, 1), (1, 0), (1, 0)) 0.08225771204699728
((3, 1), (3, 1), (3, 0), (3, 0)) 0.38272169831413727

こちらはパウリ行列の表現ではないので、今回VQEで計算できるように変換をします。

2-4-4 量子コンピュータ向け変換

ブラビキタエフ変換やジョルダンウィグナー変換がありますが、今回は前者のBK変換をしてみると、

h = bravyi_kitaev(get_fermion_operator(m.get_molecular_hamiltonian()))

確認してみると、

print(h)

(0.7407724940116754+0j)*I + (0.23528824284103544+0j)*Z[0] + (0.23528824284103542+0j)*Z[0]*Z[1] + (-0.45353118471995524+0j)*Z[2] + (-0.45353118471995524+0j)*Z[1]*Z[2]*Z[3] + (0.18421983815174378+0j)*Z[1] + (0.041128856023498556+0j)*Y[0]*Z[1]*Y[2]*Z[3] + (0.041128856023498556+0j)*X[0]*Z[1]*X[2] + (0.041128856023498556+0j)*X[0]*Z[1]*X[2]*Z[3] + (0.041128856023498556+0j)*Y[0]*Z[1]*Y[2] + (0.14020450296633868+0j)*Z[0]*Z[2] + (0.18133335898983727+0j)*Z[0]*Z[1]*Z[2]*Z[3] + (0.18133335898983727+0j)*Z[0]*Z[1]*Z[2] + (0.14020450296633868+0j)*Z[0]*Z[2]*Z[3] + (0.19136084915706864+0j)*Z[1]*Z[3]

このようにきちんと変換されました。ちょっと複雑な形をしていますが、これこそがこれまで学んできたパウリ行列の和で書かれたハミルトニアンです。そして、このhをVQEにかけます。

runner = vqe.Vqe(UCCAnsatz(h,2,Circuit().x[0]))
result = runner.run(verbose = True)

そうすると、

params: [0.00803084 0.88310007] val: -0.3685674394653079
params: [0.00803084 0.88310007] val: -0.3685674394653079
params: [1.00803084 0.88310007] val: 0.10544430487391798
params: [-1.61000316  0.88310007] val: 0.48593858308326665
params: [0.00803084 0.88310007] val: -0.3685674394653079
params: [-0.61000314  0.88310007] val: -0.24619628358384726
params: [0.38999684 0.88310007] val: -0.27563514593504546
params: [-0.07664849  0.88310007] val: -0.371744879738929
params: [-0.07664849  0.88310007] val: -0.371744879738929
params: [-0.07664849  1.88310007] val: 1.0541144286439048
params: [-0.07664849 -0.73493393] val: -0.44860542849648677
params: [-0.07664849  0.0289703 ] val: -0.8993554487133357
params: [-0.07664849  0.0289703 ] val: -0.8993554487133357
params: [-0.16132782 -0.82515946] val: -0.32853423118712466
params: [-0.07664849  0.0289703 ] val: -0.8993554487133357
params: [0.92335151 0.0289703 ] val: -0.3234253290431143
params: [-1.69468249  0.0289703 ] val: 0.7219628224155241
params: [-0.07664849  0.0289703 ] val: -0.8993554487133357
params: [-0.69468246  0.0289703 ] val: -0.5631470828007092
params: [0.30531751 0.0289703 ] val: -0.8360012977529356
params: [-0.0024937  0.0289703] val: -0.903736928997239
params: [-0.0024937  0.0289703] val: -0.903736928997239
params: [-0.0024937  1.0289703] val: -0.19293944513988606
params: [-0.0024937 -1.5890637] val: 0.5730186565494813
params: [-0.0024937  0.0289703] val: -0.903736928997239
params: [-0.0024937  -0.58906367] val: -0.6556036819150325
params: [-0.0024937  0.4109363] val: -0.7807595396151468
params: [-0.0024937  -0.00256402] val: -0.9043519838381495
params: [-0.0024937  -0.16023561] val: -0.885812805585636

出てきました。パラメータが最適化されながら最小基底状態が探索されています。

今回はansatzは量子化学の理論に基づいて電子の配置をするUCCansatzと呼ばれるものが利用されました。ansatzもきちんと理論化されているのがいいですね。

runner.ansatz.get_energy(result.circuit,runner.sampler)

-0.9043519838381495

求まりました。このように0.4でやったものを全体でやったのが最初のコードです。大きな問題はそれなりに大変そうですが、いろいろ工夫がされているようです。

VQEセミナー#2-5(中級設定編)

VQEを使い続けていると、様々な設定をしたくなってきます。今回は、以前紹介したansatzに加えて、古典最適化を見てみたいと思います。

2-5-1 古典最適化

VQEの内部構造は大きく分けて、

1、パラメータ付きの短い量子回路
2、量子回路を集計し、期待値を計算、次のパラメータの提案

次のトライアルのパラメータを渡すのは最適化アルゴリズムですが、それは古典で実行されます。

from strangeworks.blueqat import StrangeworksBackend #quantumcomputing.comを利用する時に必要
import numpy as np
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, 4)
        self.step = 1

    def get_circuit(self, params):
        a, b, c, d = params
        return Circuit().ry(a)[0].rz(b)[0].ry(c)[1].rz(d)[1]

cost = -3*q(0)-3*q(1)-2*q(0)*q(1)
cost = cost.to_expr().simplify()
minimizer=vqe.get_scipy_minimizer(method="COBYLA",options={"tol":5.0e-4})
runner = Vqe(QubitAnsatz(cost),minimizer=minimizer)
result = runner.run()

print('Result by VQE')
print(runner.ansatz.get_energy(result.circuit, runner.sampler))

# Hamiltonian to matrix
mat = cost.to_matrix()

# Calculate by numpy
print('Result by numpy')
print(np.linalg.eigh(mat)[0][0])

print('Matrix is')
print(mat)

2-5-2 最適化アルゴリズムの選択

主にscipyから選べるようです。

https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html

中には使えないものもアリそうですが、
(開発者の加藤くんが下の質問で補足してくれました!)

'Nelder-Mead', 'Powell', 'CG', 'BFGS', 'L-BFGS-B', 'TNC', 'COBYLA', 'SLSQP'

この辺りを見れば良さそうです。自分の好きなソルバーを組み込むことができます。

def hyperopt_minimizer(objective, n_params):
    from hyperopt import fmin, Trials, tpe, hp
    trials = Trials()
    best = fmin(objective, [hp.uniform(f'p{i}', 0., 2 * np.pi) for i in range(n_params)],
            algo=tpe.suggest, max_evals=10000, trials=trials, verbose=1)
    return list(best.values())

2-5-3 Qiskitでも

IBMのqiskitも参考に見てみます。インストールは、

pip install qiskit-aqua

Qiskitはちょっと長いです。

import numpy as np
import networkx as nx
import matplotlib.pyplot as plt

from qiskit.aqua.translators.ising import max_cut
from qiskit.aqua.input import EnergyInput

from qiskit import BasicAer
from qiskit.aqua.algorithms import VQE
from qiskit.aqua.components.optimizers import SPSA
from qiskit.aqua.components.variational_forms import RY
from qiskit.aqua import QuantumInstance


n=4

G=nx.Graph()
G.add_nodes_from(np.arange(0,n,1))

#エッジijに対して、(i,j,結合荷重)で設定します。
elist=[(0,1,1.0),(0,2,1.0),(0,3,1.0),(1,2,1.0),(2,3,1.0)]

G.add_weighted_edges_from(elist)
colors = ['r' for node in G.nodes()]
pos = nx.spring_layout(G)
default_axes = plt.axes(frameon=True)
nx.draw_networkx(G, node_color=colors, node_size=600, alpha=.8, ax=default_axes, pos=pos)

w = np.zeros([n,n])
for i in range(n):
    for j in range(n):
        temp = G.get_edge_data(i,j,default=0)
        if temp != 0:
            w[i,j] = temp['weight'] 
print(w)


qubitOp, offset = max_cut.get_max_cut_qubitops(w)
algo_input = EnergyInput(qubitOp)


seed = 10598

spsa = SPSA(max_trials=300)
ry = RY(qubitOp.num_qubits, depth=5, entanglement='linear')
vqe = VQE(qubitOp, ry, spsa, 'matrix')

backend = BasicAer.get_backend('statevector_simulator')
quantum_instance = QuantumInstance(backend, seed=seed, seed_transpiler=seed)

result = vqe.run(quantum_instance)


x = max_cut.sample_most_likely(result['eigvecs'][0])
print('energy:', result['energy'])
print('time:', result['eval_time'])
print('max-cut objective:', result['energy'] + offset)
print('solution:', max_cut.get_graph_solution(x))
print('solution objective:', max_cut.max_cut_value(x, w))

colors = ['r' if max_cut.get_graph_solution(x)[i] == 0 else 'b' for i in range(n)]
nx.draw_networkx(G, node_color=colors, node_size=600, alpha = .8, pos=pos)

ここで、上記はmaxcutを解いていますが、

spsa = SPSA(max_trials=300)
ry = RY(qubitOp.num_qubits, depth=5, entanglement='linear')
vqe = VQE(qubitOp, ry, spsa, 'matrix')

のところで、最適化アルゴとansatzを準備してVQEを解いているのが分かります。どのツールでも基本的には使っている回路や最適化は大きくは変わらないので、安心して量子コンピュータを学びましょう。

最後に

これはVQEセミナーという勉強会の2周目です。だいぶ資料も纏まってきましたが、もう一周して最後また資料をわかりやすくしたいと思います。

mdrft
量子コンピュータのアプリケーション、ミドルウェア、ハードウェアをフルスタックで開発
https://mdrft.com/
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした