lab-2-ja.ipynb

Kernel: Python 3

Lab 2: 変分アルゴリズム入門

このLabでは、以下を学びます。

• 特定のVQE(変分量子固有値ソルバー）における変分量子アルゴリズム

• パラメーター付き回路と二次計画法のQiskitでの作成方法と実行方法

• QAOAを使った最適化問題の解法

import networkx as nx
import numpy as np
import plotly.graph_objects as go
import matplotlib as mpl
import pandas as pd
from IPython.display import clear_output
from plotly.subplots import make_subplots
from matplotlib import pyplot as plt
from qiskit import Aer
from qiskit import QuantumCircuit
from qiskit.visualization import plot_state_city
from qiskit.algorithms.optimizers import COBYLA, SLSQP, ADAM
from time import time
from copy import copy
from typing import List
from qc_grader.graph_util import display_maxcut_widget, QAOA_widget, graphs
mpl.rcParams['figure.dpi'] = 300

1) はじめに

この章は変分量子アルゴリズムのトピックの入門であり、また以前の講義の復習です。

変分量子アルゴリズム(VQA)

最も一般的な意味では、変分量子回路は、パラメーターのセット$\theta$に依存する回路です。一般的に、変分量子アルゴリズムは、このパラメーター化された量子回路の出力を量子コンピューターに問い合わせ、その出力に基づいて、与えられたコスト関数$C(\theta)$を評価します。その後、古典的なオプティマイザーを用いて、目的関数$C$を最大化または最小化するために、回路のパラメータを更新します。これらのステップは、量子-古典ハイブリッドループで繰り返され、最終的には古典的な最適化によって最適なパラメータ$\theta^*$が見つかると終了します。

変分量子アルゴリズム(VQA)は、近未来のデバイスで量子的優位性を実現するための有望な手法とみなされることが多いです。多くのケースでは、深い量子回路の実行を必要とせず、古典的なオプティマイザーに最適化手法を委託することで、システマティックなエラーを部分的に軽減することができます。それにもかかわらず、VQAは多くの課題を抱えています。特に、VQAが効率的に学習可能であり、古典的なアルゴリズムで得られたものよりも実際に優れたソリューションを生み出すことができるかどうかという問題です。これらの課題にもかかわらず、VQAは以下のような様々な問題設定に対して提案されています。

• 変分量子固有値ソルバー(VQE): VQEは、ハミルトニアン$H$で表現された量子システムの基底状態とそれに対応するエネルギーを近似します。（つまり、対応するエルミート行列の最小の固有値と固有ベクトルです。）（こちらをご覧ください。）

• QAOA: 組み合わせ最適化問題に使われる、近似最適化アルゴリズムです。QAOAは、問題ハミルトニアンとしてコスト関数をエンコードすることで最適化問題を解くVQEとみなすことができます。（Section 4をご覧ください。）

• 変分分類器: 変分分類器は、古典的な機械学習分類器を彷彿とさせる、見えないデータサンプルを分類するためにデータセットで学習する量子回路です。

• 変分量子線形ソルバー(VQLS): VQLSは、VQEの元となる基本的なアイディアを活用して線形方程式によるシステムを解きます。

このLabでは、VQEの特別な場合であるQAOAにフォーカスします。

変分法

基底状態$\vert \psi^* \rangle$と基底状態エネルギー$E_0$をもつ量子システムを記述するエルミート演算子$H$を考えます。 変分法は、$\vert \psi^* \rangle$と$E_0$を近似する技術手法です。 これは、パラメーター化された試行状態$\vert \psi(\theta) \rangle$を選ぶことで行います。ここで、$\theta$ は、パラメーターのセットまたはベクトルです。 状態$\vert \psi \rangle$のシステムのエネルギーは以下の$H$に関する期待値によって与えられることを思い出しましょう。

$$E(\vert \psi \rangle)= \langle \psi \vert H \vert \psi \rangle$$

システムの基底状態は、最低のエネルギーの固有状態であるので、定義上、以下のようになります。すべてのパラメーター$\theta$に対して、 $$E_0 = \langle \psi^* \vert H \vert \psi^* \rangle \leq \langle \psi(\theta) \vert H \vert \psi(\theta) \rangle,$$

したがって、試行状態$\psi(\theta)$の期待値を最小化することによって、つまり期待値$\langle \psi(\theta) \vert H \vert \psi(\theta) \rangle$ ができるだけ小さくなるようなパラメーター$\theta$を見つけることによって、基底状態エネルギー$E_0$の上限を得ることができ、基底状態そのものの近似値を得ることができます。 当然ながら、良い試行状態$\psi(\theta)$の選択が、変分法の成功の鍵です。

VQE（変分量子固有値ソルバー）

変分量子固有値ソルバーは、変分法を使って、基底状態とハミルトニアン$H$の最小の固有値を近似的に求めます。試行状態は、変分量子回路によって準備された量子状態に相当し、量子コンピューターによる回路の実行をすることで、相当する期待値が測定されます。そして、古典オプティマイザーは、回路のパラメーターの調整に使われ、測定される期待値を最小化します。

与えられたハミルトニアンの基底状態に直接興味のある、化学や量子力学そのものの問題への適用から離れて、最適化問題において、ハミルトニアンの基底状態が問題の最適解に相当するコスト関数をエンコーディングすることで変分量子固有値ソルバーの概念を使うことができます。このアイディアはQAOAの核心となっています。

パラメーター付き量子回路

このセクションでは、Qiskitでパラメーター付き回路を構築し、回路のパラメーターに値を割り当てる方法を学び、独自の変分形式を構築します。

パラメーター付き回路の構築

Qiskitでパラメーター付き回路を構築することは、標準的な量子回路を作成するのとあまり違いはありません。QiskitのParameterクラスを使ってパラメーターを初期化し、回路を構築するゲートが追加された時にそれに応じて使います。以下の例では、回転量子ゲートの回転角をパラメーターとします。

from qiskit.circuit import Parameter, ParameterVector
# パラメーターを文字列識別子で初期化します。
parameter_0 = Parameter('θ[0]')
parameter_1 = Parameter('θ[1]')
circuit = QuantumCircuit(1)
# Rx、Ryゲートの回転角として初期化されたパラメーターを渡すことができます
circuit.ry(theta = parameter_0, qubit = 0)
circuit.rx(theta = parameter_1, qubit = 0)
circuit.draw('mpl')

同じパラメーターは、同じ回路の中で何回も使うことができます。上記の回路で、同じパラメーターを異なるゲートに繰り返し使うことを考えてみてください。

parameter = Parameter('θ')
circuit = QuantumCircuit(1)
circuit.ry(theta = parameter, qubit = 0)
circuit.rx(theta = parameter, qubit = 0)
circuit.draw('mpl')

便宜上、QiskitにはParameterVectorクラスが存在し、一度に複数のパラメーターを作ることができます。 以下のRealAmplitudes回路について考えてみてください。この回路は、パラメーター付き$R_Y$ゲートとエンタングルする$CX$ゲートの交互の層で構成されています。RealAmplitudes 変分形式は量子機械学習でよく使われ、Qiskitのcircuit libraryにもあります。

# 層の数と量子ビットの数を設定
n=3
num_layers = 2

# ParameterVectorsが文字列識別子とベクトルの長さを特定する整数で初期化されます
parameters = ParameterVector('θ', n*(num_layers+1))


circuit = QuantumCircuit(n, n)
for layer in range(num_layers):
    # パラメーター付きRyゲートを上記で構築したベクトルによるパラメーターを使って追加します
    for i in range(n):
        circuit.ry(parameters[n*layer+i], i)
    circuit.barrier()

    # エンタングルするCNOTゲートを追加します
    for i in range(n):
        for j in range(i):
                circuit.cx(j,i)
    circuit.barrier()

# 追加のパラメーター付きRyゲートの層を追加します
for i in range(n):
    circuit.ry(parameters[n*num_layers+i], i)
circuit.barrier()

circuit.draw('mpl')

量子回路の一部であるパラメーターを調べることができます。

print(circuit.parameters)

パラメーターに値を割り当てる

パラメーター付き回路は、パラメーターの値が割り当てられないと量子バックエンドで実行できません。そのためにQuantumCircuitメソッドを使うことができます。

assign_parameters(parameters, inplace = False)
bind_parameters(values)

bind_parametersは回路のパラメーターに数値を割り当て、新しい回路を作ります。 assign_parametersで、数値を割り当て、または他のパラメーター式を用いて、パラメーターを代用することができます。 さらに、assign_parametersでは、新しい回路を作る代わりにパラメーターを代用することができます。 回路のパラメーターに割り当てられるべき値またはパラメーター式は、ディクショナリーとして提供できます。ディクショナリー・キーは回路のパラメーターに相当し、ディクショナリーの値が値に紐づけられるか、または値の反復として提供されます。 後者においては、値は、回路に追加されたパラメーターから順にパラメーターに割り当てられます。

# バインドするランダムな値によるパラメーター・ディクショナリーを作成
param_dict = {parameter: np.random.random() for parameter in parameters}
print(param_dict)
# assign_parametersメソッドを使ってパラメーターを割り当て
bound_circuit = circuit.assign_parameters(parameters = param_dict)
bound_circuit.draw('mpl')

元の回路のパラメーターを別のパラメーター式で置き換えた次の回路を考えてみましょう。

new_parameters = ParameterVector('Ψ',9)
new_circuit = circuit.assign_parameters(parameters = [k*new_parameters[k] for k in range(9)])
new_circuit.draw('mpl')

制限付きの回路が、量子デバイス上で実行できるようになりました。パラメーター付き量子回路を、パラメーターが割り当てられていない状態で実行しようとすると、エラーが発生します。

# Aerのstatevector simulatorで割り当てられたパラメーターとともに回路を実行
simulator = Aer.get_backend('statevector_simulator')
result = simulator.run(bound_circuit).result()
statevector = result.get_statevector(bound_circuit)
plot_state_city(statevector)

# 次の行は、'circuit'に割り当てられていないパラメーターが含まれているため、実行するとエラーになります。
#result = simulator.run(circuit).result()

二次計画法

この章では、MaxCut問題を振り返り、Qiskitで二次計画法を構築する方法を学びます。

MaxCut

MaxCut問題では、入力を（おそらく重みのある）グラフとし、異なるカットからの頂点を結ぶエッジの累積重みが最大となるような、グラフの頂点を2つの離散したセットに分割する方法を見つけます。

$G=(V,E)$を$n$個の頂点を持つ（重み付き）グラフとします。頂点$v_1, \dots, v_n$に番号を付け、各頂点に値0か1のどちらかを割り当てることで、グラフの任意のカットをベクトル$x \in \{0,1\}^n$で識別することができます。そして、カット$x$の重みは、以下のコスト関数により決まります。

$$C(x) = \sum_{i,j=1}^n W_{ij} x_i (1-x_j)$$

ここで、$W_{ij} = W_{ji}$は頂点$v_i$と$v_j$を結ぶエッジの重みで、$v_i$と$v_j$がエッジでつながっていない場合は$W_{ij} = 0$とします。$x_i \neq x_j$の場合に限り、各エッジの重みが合計にちょうど1回寄与することは簡単にわかります。

以下に定義されているディクショナリーのgraphsには、様々なグラフのインスタンスが含まれているので、調べてみてください。各グラフの末尾の名前は、使われている重みの種類を指定しています。

for key in graphs.keys():
    print(key)

networkxモジュールを使ってもう一つのグラフ・インスタンスを作成し、ディクショナリーに追加してみましょう。下のウィジェットを使って、異なるグラフ・インスタンスのカットを調べることができます。 ２グループの頂点がそれぞれ水色と紺で示され、異なるグループ間の頂点の間につながっている辺（つまりカットの重みに寄与する部分）が赤で示されています。

graph = nx.Graph()
#Add nodes and edges
graph.add_nodes_from(np.arange(0,6,1))
edges = [(0,1,2.0),(0,2,3.0),(0,3,2.0),(0,4,4.0),(0,5,1.0),(1,2,4.0),(1,3,1.0),(1,4,1.0),(1,5,3.0),(2,4,2.0),(2,5,3.0),(3,4,5.0),(3,5,1.0)]
graph.add_weighted_edges_from(edges)
graphs['custom'] = graph
#Display widget
display_maxcut_widget(graphs['custom'])

Exercise 1: MaxCut

上で新たに定義したグラフについて、最大のカットを見つけて、ビットのリスト'$x=$[$x_0, x_1, ..., x_n$]'としてgraderに渡します。可能なすべてのビット列に対するコスト関数の値を棒グラフで表示する以下の関数を使うことができますが、特定のビット列$x$に対するMaxcutのコスト関数の値を計算するコードを追加する必要があります。numpy行列のweight_matrixは$W$を参照しており、weight_matrix[i,j]を介してその要素$W_{ij}$にアクセスできます。

def maxcut_cost_fn(graph: nx.Graph, bitstring: List[int]) -> float:
    """
    与えられたグラフとビット列で表現されたカットに対して、maxcutコスト関数の値を計算します
    Args:
    Args:
        graph: カット値を計算するためのグラフです
        bitstring: グラフのカットを指定する整数値'0'または'1'のリスト
    Returns:
        カットの値
    """
    # グラフの重み行列を得ます
    weight_matrix = nx.adjacency_matrix(graph).toarray()
    size = weight_matrix.shape[0]
    value = 0.

    # ここにカットの重みを計算するコードを記入します

    return value

def plot_maxcut_histogram(graph: nx.Graph) -> None:
    """
    与えられたグラフのすべての可能なカットに対する値を棒グラフにプロットします。
    Args:
        graph: カット値を計算するためのグラフです
    """
    num_vars = graph.number_of_nodes()
    # ビット列のリストと対応するカットの値を作ります
    bitstrings = ['{:b}'.format(i).rjust(num_vars, '0')[::-1] for i in range(2**num_vars)]
    values = [maxcut_cost_fn(graph = graph, bitstring = [int(x) for x in bitstring]) for bitstring in bitstrings]
    # 大きなカットの値からリストをソートします
    values, bitstrings = zip(*sorted(zip(values, bitstrings)))
    # 棒グラフを表示します
    bar_plot = go.Bar(x = bitstrings, y = values, marker=dict(color=values, colorscale = 'plasma', colorbar=dict(title='Cut Value')))
    fig = go.Figure(data=bar_plot, layout = dict(xaxis=dict(type = 'category'), width = 1500, height = 600))
    fig.show()

plot_maxcut_histogram(graph = graphs['custom'])

from qc_grader import grade_lab2_ex1
bitstring = [0, 0, 0, 0, 0, 0] # ここに正しいMAXCUTのビット列を入れてください

# grade関数は整数'0'と'1'のリストを期待していることに注意してください
grade_lab2_ex1(bitstring)

二次計画法

二次制約付き二次計画法とは、二次目的関数と、最適化変数に対する線形制約および二次制約を持つ最適化問題です。つまり、次のような形で書くことができます。 $$\begin{align} \text{minimize} &&x^T Q x + c^T x &&\\ && && \\ \text{subject to} &&Ax \leq b &&\\ && x^TQ_ix + a_i^Tx \leq r_i \\ && l_j \leq x_j \leq u_j \\ \end{align}$$

ここで、$Q \in \mathbb{R}^{n \times n}$は対称的な実数の値の$n \times n$-行列、$c \in \mathbb{R}^n$は、目的関数の二次関数部分と一次関数部分を指定する実数値のベクトルです。最適化変数$x_i, i \in \{1, \dots, n\}$は、問題に応じて、バイナリー、整数、実数の値のいずれかです。

Qiskitでは、qiskit_optimizationモジュールにあるQuadraticProgramクラスを使って二次計画法を作成することができます。新しいモデルを作成するには、単に希望する名前でイニシャライザーを呼び出します。

from qiskit_optimization import QuadraticProgram
quadratic_program = QuadraticProgram('sample_problem')
print(quadratic_program.export_as_lp_string())

binary_var、integer_var、continuous_varというそれぞれのメソッドを呼び出すことで、バイナリー、整数、連続の最適化変数を追加することができます。二次計画法に追加された変数はすべて名前で指定されます。また、整数型や連続型の変数の境界を指定するには、lowerboundとupperboundという引数を使います。

quadratic_program.binary_var(name = 'x_0')
quadratic_program.integer_var(name = 'x_1')
quadratic_program.continuous_var(name = 'x_2', lowerbound = -2.5, upperbound = 1.8)

目的関数を設定するには、最小化問題か最大化問題かに応じて、minimize または maximizeというメソッドを使用します。linearとquadraticの引数には、リスト、行列、ディクショナリーのいずれかを渡すことで、二次項と線形項を指定できます。また、引数constantで定数のオフセットを指定することも可能です。

最後に、二次計画法を読みやすい形式で表示するには、export_as_lp_stringメソッドを使用します。

quadratic = [[0,1,2],[3,4,5],[0,1,2]]
linear = [10,20,30]
quadratic_program.minimize(quadratic = quadratic, linear = linear, constant = -5)
print(quadratic_program.export_as_lp_string())

二次計画法を一から作るのではなく、既存のdocplexモデルを変換することもできます。詳しくは、二次計画法のチュートリアル で説明しています。

二次計画法としてのMaxCut

任意のMaxCut問題を二次計画法で定式化することができます。対称的な重み行列$W$を持つMaxCut問題のコスト関数をもう一度考えてみましょう。

$$C(x) = \sum_{i,j=1}^n W_{ij} x_i (1-x_j).$$

これは明らかに二次コスト関数であり、上に書いたように二次計画法の標準的な形式で再構成することができます。 $$\begin{align} \sum_{i,j=1}^n W_{ij} x_i (1-x_j) &= \sum_{i,j=1}^n W_{ij} x_i - W_{ij}x_i x_j \\ &= \sum_{i=1}^n \left( \sum_{j=1}^n W_{ij} \right) x_i - \sum_{i,j = 1}^n W_{ij}x_i x_j \\ &= c^T x + x^T Q x, \\ \end{align}$$ $Q$と$c$については、以下のように定義されます。 $$Q_{ij} = -W_{ij} \qquad c_i = \sum_{j=1}^n W_{ij}.$$ このようにして、バイナリー変数、二次目的関数を持ち、変数制約のない最適化インスタンスが得られます。このような形式の二次計画法は、二次制約なし二値最適化インスタンス（quadratic unconstrained binary optimization instance）、略してQUBOとも呼ばれます。次のセクションでは、QAOAを使ってこのタイプの問題を最適化する方法を学びます。

Exercise 2: MaxCutからQUBOへ

以下の関数は、graphで定義されたMaxCut問題インスタンスから二次計画法を構築するものです。前のセクションで説明した方法を使って、コードを完成させて下さい。二次計画法にバイナリー変数を'x_0', 'x_1', ..., 'x_n'という名前をつけて追加する必要があります。ここでは、重み行列$W$をweight_matrixと呼び、qubo行列$Q$をqubo_matrixと呼びます。

def quadratic_program_from_graph(graph: nx.Graph) -> QuadraticProgram:
    """MaxCut問題インスタンスに対して、与えられたグラフから二次計画法を構築します。
    Args:
        graph: 問題となるグラフ
    Returns:
        QuadraticProgram
    """
    # グラフの重み行列を得る
    weight_matrix = nx.adjacency_matrix(graph)
    shape = weight_matrix.shape
    size = shape[0]
    # 重み行列Wからqubo行列Qを構築
    qubo_matrix = np.zeros((size, size))
    qubo_vector = np.zeros(size)
    for i in range(size):
        for j in range(size):
            qubo_matrix[i, j] -= weight_matrix[i, j]
    for i in range(size):
        for j in range(size):
            qubo_vector[i] += weight_matrix[i,j]

    # ここにコードを記入します
    
    return quadratic_program

quadratic_program = quadratic_program_from_graph(graphs['custom'])
print(quadratic_program.export_as_lp_string())

from qc_grader import grade_lab2_ex2

# grade関数は二次計画法を期待していることに注意してください
grade_lab2_ex2(quadratic_program)

QAOA

このセクションでは、Qiskitに実装されているQuantum Approximate Optimization Algorithm (QAOA)を用いて、QAOAの変分形式の独自のバージョンを実装し、上のセクションで定義したMaxCut問題を解きます。

QAOA の変分形式

QAOAの変分形式は、次のような層状構造になっていることを思い出してください。

ここで、$H_C$は最適化問題のコスト関数を符号化したコストハミルトニアン、$H_M$はミキサーハミルトニアンを指します。
オリジナル版のQAOAでは、初期準備状態は、均等な重ね合わせ状態 $$\vert + \rangle = \sum_{x \in \{0,1\}^{n}} \frac{1}{\sqrt{2^n}} \vert x \rangle$$ であり、ミキサーハミルトニアンは、すべての量子ビットの単一パウリ$X$演算子の和として以下のように定義されます。 $$H_M = \sum_{i=1}^n{X_i}.$$

対応する行列$Q$を持つQUBOインスタンスに対して、我々はコスト関数$C(x)$をコストハミルトニアン$H_C$として符号化したいです。つまり、すべての計算基底状態$x \in \{0,1\}^n$に対して以下のような演算子$H_C$を見つけようとしています。

$$H_C \vert x \rangle = \left( x^T Q x + c^T x \right) \vert x \rangle = \left( \sum_{i,j=1}^n x_i Q_{ij} x_j + \sum_{i=1}^n c_ix_i \right) \vert x \rangle$$

以下の事実を使います。 $$Z_i\vert x \rangle = (-1)^{x_i} \vert x \rangle = \left( 1-2x_i \right) \vert x \rangle \implies x_i \vert x \rangle = \frac{\mathbb{I}-Z_i}{2}\vert x \rangle,$$

ここで、$\mathbb{I}$は単位演算子、$Z_i$ は量子ビット$i$のPauli-Z行列を表します。コスト関数の式$C(x)$において、$x_i \to \frac{\mathbb{I}-Z_i}{2}$という置換を行うことで、$H_C$をPauli-Z演算子の総和として書くことができます。これにより、次のような$H_C$の式が得られます。

$$H_C = \sum_{i,j=1}^n \frac{1}{4}Q_{ij} Z_iZ_j - \sum_{i=1}^n \frac{1}{2} \left(c_i + \sum_{j=1}^n Q_{ij} \right) Z_i + \left( \sum_{i,j=1}^n \frac{Q_{ij}}{4} + \sum_{i=1}^n \frac{c_i}{2} \right)$$

両方のハミルトニアンを指数化すると、コスト層に$R_Z$と$R_{ZZ}$のゲート、ミキサー層に$R_X$のゲートからなる変分法が得られます。 $$e^{- i \beta H_M} = \prod_{i=1}^n R_{X_i}(2\beta) \quad \quad e^{- i \gamma H_C} = \prod_{i,j=1 \\ i \neq j}^n R_{Z_iZ_j} \left( \frac{1}{2} Q_{ij} \gamma \right) \prod_{i=1}^nR_{Z_i} \left( \left( c_i + \sum_{j=1}^n Q_{ij} \right) \gamma \right)$$ 例として、上で定義したMaxCutインスタンスのQUBOに対するQAOA変分形式（$p=1$の場合）を考えます。

上の回路の$R_{ZZ}$ゲートは、2つの$CNOT$ゲートと1つのシングル回転$R_Z$ゲートの組み合わせとして分解されています。また、MaxCut QUBOインスタンスのQAOA回路では、シングル回転$R_Z$ゲートの角度は常に0に等しいことに注意してください。

Exercise 3: QAOAの変分形式

二次計画法とパラメータ$p$を入力として受け取り、出力として対応する$p$層のQAOA回路を生成する関数を書いてください。 関数の基本的な構成は以下の通りですが、コスト層やミキサー層として適用される部分を作成する必要があります。 また、上の説明の最後の式にあるように、回転ゲートの角度を計算する必要があります。 以下のコードでは、qubo行列$Q$をqubo_matrix、ベクトル$c$をqubo_linearityと表記しています。回路の中に測定が入らないように注意してください。

def qaoa_circuit(qubo: QuadraticProgram, p: int = 1):
    """
    QUBOインスタンスと層数pが与えられ、対応するパラメータ付きQAOA回路をp層で構築します。
    Args:
        qubo: 二次計画法のインスタンス
        p: QAOA回路の層の数
    Returns:
        パラメーター化されたQAOA回路
    """
    size = len(qubo.variables)
    qubo_matrix = qubo.objective.quadratic.to_array(symmetric=True)
    qubo_linearity = qubo.objective.linear.to_array()

    # 量子レジスターと古典レジスターを準備
    qaoa_circuit = QuantumCircuit(size,size)
    # 全ての量子ビットにアダマールゲートの初期層を適用
    qaoa_circuit.h(range(size))

    # 回路中で使われるパラメーターを作成
    gammas = ParameterVector('gamma', p)
    betas = ParameterVector('beta', p)

    # 各層を作成する外側のループ    
    for i in range(p):

        # コスト層からR_Z回転ゲートを適用
        # ここにコードを記入します

        # コスト層からエンタングルした量子ビットの回転としてR_ZZ回転ゲートを適用
        # ここにコードを記入します
                        
        # 全ての量子ビットに2*beta_iの角度のX回転を適用
        # ここにコードを記入します

    return qaoa_circuit

上記の演習を終えたら、さまざまなグラフインスタンスのQAOA回路を調べてみましょう。

quadratic_program = quadratic_program_from_graph(graphs['custom'])
custom_circuit = qaoa_circuit(qubo = quadratic_program)

test = custom_circuit.assign_parameters(parameters=[1.0]*len(custom_circuit.parameters))

from qc_grader import grade_lab2_ex3

# grade関数は量子回路を期待していることに注意してください
grade_lab2_ex3(test)

QAOAクラスとMinimumEigenOptimizerクラス

これまでの演習を完了された方、おめでとうございます！ あなたは、一般的なQUBOインスタンス、特にMaxCut問題を量子コンピューター上で実行して解くことができる独自QAOAを実装しました。Qiskitは、たった数行のコードで最適化問題を解く独自のQAOAの実装も提供しています。

QiskitにおけるQAOAのクラスは、qiskit.algorithmsの中にあり、変分固有値ソルバーのクラスVQEを直接継承しています。 QAOAのイニシャライザーは、特に以下の入力パラメーターを受け取ります。

• optimizer: この引数は、回路のパラメーターをアップデートするための古典オプティマイザーをさします。Qiskitは、数多くの異なるoptimizersを提供しており、QiskitのNative Optimizer クラスをサブクラス化することで独自のオプティマイザーを定義することもできます。Qiskitで提供されている基本的なオプティマイザーは以下のものがあります。
  

  • COBYLA : 線形近似最適化による制限付きオプティマイザー。

  • SLSQP : シーケンシャル最小二乗プログラミング・オプティマイザー。

  • ADAM : 勾配ベースの最適化アルゴリズムで、低次モーメントの適応推定に依存しています。

• reps: QAOAの変分形式における層の数$p$、つまり、アルゴリズムの深さ（デプス）です。$p$の値が高いと、QAOAは理論的には、より良い結果を得ますが、量子回路は深くなり、最適化すべきパラメーターが増えます。

• initial_point: 最適化を始めるための、最適化された初期パラメーターの値（$\beta$と$\gamma$の値）。

• quantum_instance:アルゴリズムを実行させる量子インスタンス。これは、シミュレーターかまたは、実際のハードウェアデバイスです。

• callback: 最適化の過程をモニターするための呼び出し関数。

以下で、QAOAの最適化過程において、これらの引数の影響をテストすることができます。

MinimumEigenOptimizerオブジェクトは、最適化プロセスにおけるラッパーを提供し、対応する最適化の結果を得るために、二次計画法から、QAOAなどの与えられたMinimumEigenSolverと呼ばれる量子ビット演算子への変換を処理します。イニシャライザーは、最適化に使用するMinimumEigenSolverを入力として受け取り、二次計画法を入力とするoptimizeメソッドを呼び出すことで最適化が実行されます。全てを一つにまとめると、上で定義したMaxCut問題の解をたった数行のコードで求めることができます。

from qiskit.algorithms import QAOA
from qiskit_optimization.algorithms import MinimumEigenOptimizer
backend = Aer.get_backend('statevector_simulator')
qaoa = QAOA(optimizer = ADAM(), quantum_instance = backend, reps=1, initial_point = [0.1,0.1])
eigen_optimizer = MinimumEigenOptimizer(min_eigen_solver = qaoa)
quadratic_program = quadratic_program_from_graph(graphs['custom'])
result = eigen_optimizer.solve(quadratic_program)
print(result)

ご覧のように、わずかな回数の最適化ステップで最適値を得ることができます。最適化された状態ベクトルを表示することで、パラメーターの最適化が完了した後に、どのビット列がよりサンプリングされる可能性が高いかを確認できます。

def plot_samples(samples):
    """
    量子アルゴリズムのサンプルに対する棒グラフをプロットする
    Args:
        samples
    """
    # サンプルを確率でソートします
    samples = sorted(samples, key = lambda x: x.probability)
    # 確率、関数の値、ビット列のリストを得ます
    probabilities = [sample.probability for sample in samples]
    values = [sample.fval for sample in samples]
    bitstrings = [''.join([str(int(i)) for i in sample.x]) for sample in samples]
    # 棒グラフをプロットします
    sample_plot = go.Bar(x = bitstrings, y = probabilities, marker=dict(color=values, colorscale = 'plasma',colorbar=dict(title='Function Value')))
    fig = go.Figure(
        data=sample_plot, 
        layout = dict(
            xaxis=dict(
                type = 'category'
            )
        )
    )
    fig.show()

plot_samples(result.samples)

QAOAの最適化プロセスの探究

以下のウィジェットは、上で定義されたMaxCut問題についてすでに計算されたエネルギーの分布を表し、深さ$p=1$におけるQAOAの最適化の過程を可視化しています。パラメーターを更新するたびに、状態ベクトルと平均関数値がどのように変化するかを示しています。 別の古典オプティマイザーや初期ポイントを使って、最適化プロセスにどれだけ影響があるかみてみましょう。initial_pointセットにおける最初の値は、パラメーター$\gamma$の初期値を設定し、2個目の値はパラメーター$\beta$の値を設定することに注意してください。

graph_name = 'custom'
quadratic_program = quadratic_program_from_graph(graph=graphs[graph_name])

trajectory={'beta_0':[], 'gamma_0':[], 'energy':[]}
offset = 1/4*quadratic_program.objective.quadratic.to_array(symmetric  = True).sum() + 1/2*quadratic_program.objective.linear.to_array().sum()
def callback(eval_count, params, mean, std_dev):
    trajectory['beta_0'].append(params[1])
    trajectory['gamma_0'].append(params[0])
    trajectory['energy'].append(-mean + offset)

optimizers = {
    'cobyla': COBYLA(),
    'slsqp': SLSQP(),
    'adam': ADAM()
}
qaoa = QAOA(optimizer = optimizers['cobyla'], quantum_instance = backend, reps=1, initial_point = [6.2,1.8],callback = callback)
eigen_optimizer = MinimumEigenOptimizer(min_eigen_solver = qaoa)
result = eigen_optimizer.solve(quadratic_program)
fig = QAOA_widget(landscape_file=f'./resources/energy_landscapes/{graph_name}.csv', trajectory = trajectory, samples = result.samples)
fig.show()

上にプロットされたコスト関数のランドスケープに周期性があることに注目してください。 これは、QAOAの変分形式で使われる回転ゲートに周期性があることの直接的な結果です。 全ての回転ゲートが$2\pi$の周期性を持ち、我々が調べているMaxCutの例が、整数の重みのみをエッジに持つので、パラメーター$\gamma$の$2 \pi$の周期性が、パラメーター$\beta$の$\pi$の周期性が認められます。

$p$の値がより高い時

より大きな値の$p$においては、QAOAは、理論的にコストハミルトニアンに対してより良いエネルギー値を持つ量子状態に到達することができます。なぜなら、層の数を増やすと、到達できる量子状態の空間が厳密に増えるからです。つまり、もし

$$M_p = \max_{\beta, \gamma \in \mathbb{R}^p} \langle \psi_p(\beta, \gamma) \vert H_C \vert \psi_p(\beta, \gamma) \rangle$$

が深さ$p$のQAOA変分形式で最大エネルギー値に到達することを示す場合、以下が成り立ちます。 $$M_{p+1} \geq M_p.$$

実際、断熱定理を使えば、無限の極限では、以下のコスト関数の最大値$C_{max}$に到達することを証明できます。 $$\lim\limits_{p \to \inf} M_p = C_{max}.$$

断熱量子コンピューターとの関連性は、QAOAが良い結果をもたらすと期待する理由の正当性を示しますが、それが成立するのは無限の限界においてのみであり、それはつまり、無限の長い量子アルゴリズムに対応するものです。 また、層の数の増加とともにパラメーターの数が増え、パラメーターの数が増えれば増えるほど大域的な最適解を見つけるのが難しくなるという事実にも注意する必要があります。

以下のプロットは、$p$の値の増加に伴う最適なQAOAの状態の進化を示し、回路の深さがQAOAのパフォーマンスに与える影響について示しています。

graph_name = 'custom'
quadratic_program = quadratic_program_from_graph(graphs[graph_name])
# 進化のすべての数を記録するために呼び出しを作ります
max_evals = 0
def callback(eval_count, params, mean, std_dev):
    global max_evals
    max_evals = eval_count

# 値を記録するための空のリストを作ります
energies = []
runtimes = []
num_evals=[]
    
# 異なるpの値においてQAOAを実行します
for p in range(1,10):
    print(f'Evaluating for p = {p}...')
    qaoa = QAOA(optimizer = optimizers['cobyla'], quantum_instance = backend, reps=p, callback=callback)
    eigen_optimizer = MinimumEigenOptimizer(min_eigen_solver = qaoa)
    start = time()
    result = eigen_optimizer.solve(quadratic_program)
    runtimes.append(time()-start)
    num_evals.append(max_evals)
    # サンプルから最終状態のエネルギーを計算します
    avg_value = 0.
    for sample in result.samples:
        avg_value += sample.probability*sample.fval
    energies.append(avg_value)

# プロットを作成して表示します
energy_plot = go.Scatter(x = list(range(1,10)), y =energies, marker=dict(color=energies, colorscale = 'plasma'))
runtime_plot = go.Scatter(x = list(range(1,10)), y =runtimes, marker=dict(color=runtimes, colorscale = 'plasma'))
num_evals_plot = go.Scatter(x = list(range(1,10)), y =num_evals, marker=dict(color=num_evals, colorscale = 'plasma'))
fig = make_subplots(rows = 1, cols = 3, subplot_titles = ['Energy value', 'Runtime', 'Number of evaluations'])
fig.update_layout(width=1800,height=600, showlegend=False)
fig.add_trace(energy_plot, row=1, col=1)
fig.add_trace(runtime_plot, row=1, col=2)
fig.add_trace(num_evals_plot, row=1, col=3)
clear_output()
fig.show()

CVaR

最近提案されたQAOAの適応作の一つは、リスクのある条件付き価値（CVaR）を計算するアイディアです。ここでは、コスト関数$f$を計算するために、回路の測定結果から、全てのカットの値$c_i$の平均を計算する代わりに、アルゴリズムは、測定された最大のカット値の一部$\alpha$のみを考慮します。

$$f(\theta) = \frac{1}{n}\sum_{i=1}^n c_i \to f(\theta) = \frac{1}{\lceil \alpha n \rceil}\sum_{i=1}^{\lceil \alpha n \rceil} c_i,$$

ここで、$c_i$は大きさが多い側に並んでいます。 我々が興味があるのは、最適化問題の最適解を一つ得ることだけなので、最適なカットの一部のみを見ることで最適化プロセスをスピードアップさせることができます。 CVaRを使うことで、QAOAインスタンスのエネルギーランドスケープがどのように変化するか見てみましょう。 以下のコードは、上のexerciseであなたが作成したコードを使って、与えられたMaxCutインスタンスについてQAOAランドスケープを作成しプロットし、グリッドサーチにおいて最適なパラメーターを返します。 エネルギーランドスケープの計算には少し時間がかかります。

def plot_qaoa_energy_landscape(graph: nx.Graph, cvar: float = None):
    num_shots = 1000
    seed = 42
    simulator = Aer.get_backend('qasm_simulator')
    simulator.set_options(seed_simulator = 42)

    # 回路の作成
    circuit = qaoa_circuit(qubo = quadratic_program_from_graph(graph), p=1)
    circuit.measure(range(graph.number_of_nodes()),range(graph.number_of_nodes()))

    # 全てのビット列に対して既に計算されたカット値でディクショナリーを作成 
    cut_values = {}
    size = graph.number_of_nodes()
    for i in range(2**size):
        bitstr = '{:b}'.format(i).rjust(size, '0')[::-1]
        x = [int(bit) for bit in bitstr]
        cut_values[bitstr] = maxcut_cost_fn(graph, x)

    # 全てのパラメーターにおいてグリッドサーチを実行
    data_points = []
    max_energy = None
    for beta in np.linspace(0,np.pi, 50):
        for gamma in np.linspace(0, 4*np.pi, 50):
            bound_circuit = circuit.assign_parameters([beta, gamma])
            result = simulator.run(bound_circuit, shots = num_shots).result()
            statevector = result.get_counts(bound_circuit)
            energy = 0
            measured_cuts = []
            for bitstring, count in statevector.items():
                measured_cuts =  measured_cuts + [cut_values[bitstring]]*count

            if cvar is None:
                # 全てのカット値の平均を計算します
                energy = sum(measured_cuts)/num_shots
            else:
                raise NotImplementedError()
                # ここにコードを入れます

            # 最適なパラメーターを更新します
            if max_energy is None or energy > max_energy:
                max_energy = energy
                optimum = {'beta': beta, 'gamma': gamma, 'energy': energy}

            # データを更新します
            data_points.append({'beta': beta, 'gamma': gamma, 'energy': energy})


    # data_pointsから面を作成し表示します
    df = pd.DataFrame(data_points)
    df = df.pivot(index='beta', columns='gamma', values='energy')
    matrix = df.to_numpy()
    beta_values = df.index.tolist()
    gamma_values = df.columns.tolist()

    surface_plot = go.Surface(
        x=gamma_values, 
        y=beta_values,
        z=matrix,
        coloraxis = 'coloraxis'
    )
    fig = go.Figure(data = surface_plot)
    fig.show()

    # 最適値を返します
    return optimum

graph = graphs['custom']
optimal_parameters = plot_qaoa_energy_landscape(graph = graph)
print('Optimal parameters:')
print(optimal_parameters)

Exercise 4: CVaR

上記のコードを適応して、サンプルされたすべてのカット値の平均を取る代わりに、リスク時の条件付き値を計算するアルゴリズムで、パラメータ$\alpha$の設定を変化させることでQAOAのエネルギー分布がどのように変わるかを観察します。1000ショットを用いてCVaRパラメータを$\alpha = 0.2$に設定し、qasm simulatorのrandom seedを42に設定した場合（上記のコードで既に設定しています）、グリッド探索によって得られるエネルギーと最適なパラメーターはどのようになるでしょうか。

optimal_parameters = plot_qaoa_energy_landscape(graph = graph, cvar = 0.2)
print(optimal_parameters)

from qc_grader import grade_lab2_ex4

# grade関数は、エントリー'beta','gamma','energy'を持つPythonディクショナリーを
# 期待していることに注意してください。
grade_lab2_ex4(optimal_parameters)

その他の参考資料

• QiskiteテキストブックにおけるVQEのチュートリアル

• QiskiteテキストブックにおけるQAOAによる組み合わせ最適化のチュートリアル

• E. Farhi, J. Goldstone, and S. Gutmann, A Quantum Approximate Optimization Algorithm (2014): QAOAの最初の論文

• Qiskit Optimizationチュートリアル 以下を含む

  • 二次計画法

  • 二次計画法コンバーター

  • 最小固有オプティマイザー

• サイモンズ計算機理論研究所におけるE. Farhiによる QAOAの最近の結果 についてのスピーチ