例題1

今回は2×2の正方行列という簡単な例にしよう。

$$ A = \begin{pmatrix} 2 & -1\\ -1 & 5 \end{pmatrix}
$$

この固有値と固有ベクトルは簡単に求まる。固有値の解析解は以下の通り。
$$ A - \lambda I = 0 $$

$$\Leftrightarrow
\lambda^2 +3 \lambda -11 = 0 $$

$$\Leftrightarrow\
\lambda = -\frac{1}{2} (3 \pm \sqrt{53} )
$$

$ \sqrt{53}$は7より少し大きい位の値だとすると、最小の固有値の方は-5位の値と目算がつく。

厳密な値をpythonで書いて確かめよう。

import numpy as np
matrix = np.array([[2, -1],
                   [-1, -5]]
                   )
# np　で固有値と固有ベクトルを計算
eigenvalues, eigenvectors = np.linalg.eig(matrix)

# 結果の表示
print("Matrix:")
print(matrix)
print("\nEigenvalues: ", eigenvalues)
print("\nEigenvectors: \n", eigenvectors)

出力は以下の通り。

Matrix:
[[ 2 -1]
 [-1 -5]]

Eigenvalues:  [ 2.14005494 -5.14005494]

Eigenvectors: 
 [[ 0.99033427  0.13870121]
 [-0.13870121  0.99033427]]

ということで、最小の固有値は-5.140となり、上の解析解の結果と整合的である。

VQEによる解法

この固有値問題を量子コンピューティングの手法で解こうとするとどうなるか。ここではVQE（Variational Quantum Eigensolver）という手法を用いる。VQEは、量子コンピュータを使用して化学や物理学などの問題の固有値（最小エネルギー）を求めるためのアルゴリズムの一種。VQEはノイズのある中規模の量子デバイス上で動作することが期待されており、量子コンピューティングの利用例の一つ。

VQEは変分法を基盤としており、次のような基本的な手順で動作する。

• 1. アンザッツの選択: 問題に対して適切な変分フォーム（アンザッツ）を選択。アンザッツは、量子回路内のパラメータを調整することで構成され、問題に応じて最適な形状を見つけることが期待される
• 1. 期待値の評価: アンザッツでパラメータを調整した後、期待値を計算。これにより、アンザッツで表現された状態のエネルギーが得られる
• 1. 古典的最適化: 計算された期待値を最小化するようにアンザッツのパラメータを更新。これは通常、古典的な最適化アルゴリズムを使用して行う
• 1. 収束確認: アルゴリズムが収束するか、または所望の精度に達するまで、ステップ2-3を繰り返す

Qiskitで量子コンピューティング

上記の行列をっ量子コンピューティングの方法で求めてみよう。

まず前提となるパッケージを以下のようにインストールしておこう。

pip install qiskit-algorithms
pip install qiskit-aer
pip install pylatexenc

この上で、まずは必要なライブラリを呼び出す。

from qiskit import Aer
from qiskit.algorithms import VQE
from qiskit.algorithms.optimizers import COBYLA
from qiskit.circuit.library import RealAmplitudes
from qiskit.opflow import MatrixOp
import numpy as np

ここでVQEというのが固有値問題を解くためのソルバーになる。
COBYLAがオプティマイザーの選択だが本記事では深入りはしない。
MatrixOpが与えられた行列を複素化してハミルトニアンを生成するモジュール。

# 行列を定義
matrix = np.array([[2, -1],
                   [-1, -5]]
                   )

# MatrixOpオブジェクトを作成
hamiltonian = MatrixOp(matrix)
print(hamiltonian)

を行うと、以下の結果を得る。

Operator([[ 2.+0.j, -1.+0.j],
          [-1.+0.j, -5.+0.j]],
         input_dims=(2,), output_dims=(2,))

実数の行列が複素化されていることが確認出来る。但し、複素成分はゼロなので実質なにも変わっていない。このハミルトニアンの固有値問題をVQEで解くことで元の行列の固有値問題を解くことを目指す。


以下のコードでVQE用の量子回路を構築する。VQEモジュールの良いところは、量子回路の詳細について自分で設計しなくても済む点であるが、悪い点は何をやっているのか分からない点だ。上手く収束しない場合など意味を考える際に何をやっているか分からないと想像力も働かないので利用する場合は周辺領域を良く勉強しておかなければならない。

# VQE用の量子回路を構築
num_qubits = len(matrix)
ansatz = RealAmplitudes(num_qubits, reps=1)
optimizer = COBYLA(rhobeg=0.5, tol=0.0001)

from qiskit.utils import QuantumInstance
seed = 170
backend = Aer.get_backend('qasm_simulator')
qi = QuantumInstance(backend=backend, seed_simulator=seed, seed_transpiler=seed)

vqe_algorithm = VQE(ansatz=ansatz, optimizer=optimizer, quantum_instance=qi)

ここで最適化のロジックはCOBYLA（Constrained Optimization BY Linear Approximations）を用いている。COBYLAは制約最適化問題を解くための数値計算アルゴリズムで制約の近似を用い、各反復で線形制約のもとで目的関数を最小化しようとする。オプティマイザーの選択方法によって、解への収束が大分異なるのでいくつかの方法を試してみる方が無難だ。


また、バックエンドのシミュレーターは'qasm_simulator'を用いている。これは、古典的なシミュレーションを提供するシミュレーターで、量子ゲートのエラーやノイズを考慮しない理想的な量子回路の挙動をシミュレートする。このシミュレーターは、量子アルゴリズムの基本的な理解や初期の検証に使用される。よって実際に実機に投入するとこのような理想的な結果は出てこない点に注意が必要だ。


固有値問題を解いて結果を表示しているのが以下のコードになる。
今回は500回収束実験を行い、その結果となる固有値をEVaに格納した。

# 固有値問題を解く
EVa = []

for i in range(500):
    result = vqe_algorithm.compute_minimum_eigenvalue(hamiltonian)
    eigenvalue = result.eigenvalue.real
    eigenvector = result.eigenstate
    EVa.append(eigenvalue)
    
# 結果の表示
print("Matrix:")
print(matrix)
print("\nEigenvalue:")
print(eigenvalue)
print("\nEigenvector:")
print(eigenvector)

以下が結果となる。但し、表示されているのか500回目の試行結果である。

Matrix:
[[ 2 -1]
 [-1 -5]]

Eigenvalue:
-5.140046712248724

Eigenvector:
{'0': 0.13975424859373686, '1': 0.9901862198596787}

理想的な固有値の最小値が-5.14005494　　
量子回路による計算結果が-5.14004671
なのでかなり高い精度で結果が得られたことが分かる。

次に得られた結果を描画してみよう。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
    
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.hist(EVa, bins=50, histtype='barstacked', ec='black')
plt.xlabel("eigenvalue")
plt.ylabel("number")
plt.show()
print("Max:", np.array(EVa).max())
print("Min:", np.array(EVa).min())

上記コードによって得られた描画結果が以下となる。
-5.14の近傍に殆どの試行結果が集中しているが、時々外れた値が見られるが、おしなべて良い精度で収束していることが確認される。

発展的問題について

上記は非常に良好で分かり易い結果を表示している。
しかし、様々な条件次第では簡単には収束しない場合もあるので上記コードを安易に用いても期待通りの結果は得られない可能性がある。その場合はよく意味を考えて単年に取り組むことが必要になる。あるいは専門家の助けを求めた方が現実的な場合もあるだろう。

付録

ここまでの結果をまとめたのが以下のコードになる。

from qiskit import Aer
from qiskit.algorithms import VQE
from qiskit.algorithms.optimizers import COBYLA
from qiskit.circuit.library import RealAmplitudes
from qiskit.opflow import MatrixOp
import numpy as np

# 行列を定義
matrix = np.array([[2, -1],
                   [-1, -5]]
                   )

# MatrixOpオブジェクトを作成
hamiltonian = MatrixOp(matrix)

# VQE用の量子回路を構築
num_qubits = len(matrix)
ansatz = RealAmplitudes(num_qubits, reps=1)
optimizer = COBYLA(rhobeg=0.5, tol=0.0001)

from qiskit.utils import QuantumInstance
seed = 170
backend = Aer.get_backend('qasm_simulator')
qi = QuantumInstance(backend=backend, seed_simulator=seed, seed_transpiler=seed)

vqe_algorithm = VQE(ansatz=ansatz, optimizer=optimizer, quantum_instance=qi)

# 固有値問題を解く
EVa = []

for i in range(500):
    result = vqe_algorithm.compute_minimum_eigenvalue(hamiltonian)
    eigenvalue = result.eigenvalue.real
    eigenvector = result.eigenstate
    EVa.append(eigenvalue)
    
# 結果の表示
print("Matrix:")
print(matrix)
print("\nEigenvalue:")
print(eigenvalue)
print("\nEigenvector:")
print(eigenvector)