注釈

このページは docs/tutorials/06_examples_max_cut_and_tsp.ipynb から生成されました。

最大カット問題とセールスマン巡回問題#

はじめに#

金融や工学などの定量的な分野における多くの問題は最適化問題です。最適化問題は、複雑な意思決定と戦略定義の中心にあります。

最適化(または組合せ最適化) とは、有限あるいは加算無限の潜在的解決策集合から、最適な解決策を探索することを意味します。最適性は、最小化または最大化するいくつかの基準関数に対して定義されます。これは通常、コスト関数あるいは目的関数と呼ばれています。

典型的な最適化問題

最小化: コスト、距離、トラバーサルの長さ、重量、処理時間、資材、エネルギー消費量、オブジェクト数

最大化: 利益、価値、生産高、リターン、利回り、有用性、効率、容量、オブジェクトの数

ここでは、多くの分野で実用的関心のある最大カット問題を考えます。それらをどのようにすれば量子コンピューターに手動でマッピングできるか、あるいはQiskitの最適化モジュールがこれをどのようにサポートしているかを示します。

重み付き最大カット#

最大カットはNP完全問題であり、クラスタリング、ネットワーク・サイエンス、統計物理学に応用されています。実用的なアプリケーションが特定の最大カット・インスタンスにどのようにマッピングされるかを把握するために、多くの人々が相互に作用し影響を及ぼしあうことのできるシステムを考えてみてください。個人はグラフの頂点によって表され、それらの相互作用はグラフの頂点またはエッジ間のペアワイズ接続として見られます。この表現を念頭に置くと、典型的なマーケティングの問題を簡単にモデル化できます。たとえば、個人がお互いの購買決定に影響を与えると想定し、それらがお互いにどれほど強い影響を与えるかについての知識が与えられているとします。影響は、グラフの各エッジに割り当てられた重みによってモデル化できます。次に、製品を一部の個人に無料で提供するマーケティング戦略の結果を予測し、収益を最大化するために、無料の製品を入手する必要がある個人の最適なサブセットを問うことができます。

この問題の正式な定義は以下のとおりです。:

\(n\)-node の無向グラフ G = (V, E) を考えます。ここでは、 |V| = n\((i, j)\in E\) におけるエッジの重みは \(w_{ij}>0\)\(w_{ij}=w_{ji}\) です。カットは、元のセットVを2つのサブセットにするパーティションとして定義されます。この場合、最適化されるコスト関数は、2つの異なるサブセットの点を接続しており、カットと 交差 しているエッジの重みの合計です。 \(x_i=0\) または \(x_i=1\) を各ノード \(i\) に割り当てることで、全体の利潤関数を最大化しようとします(ここと次の合計では、インデックス 0,1,…n-1 に対して実行されます)。

\[\tilde{C}(\textbf{x}) = \sum_{i,j} w_{ij} x_i (1-x_j).\]

単純なマーケティングモデルでは、 \(w_{ij}\) は人 \(i\) が無料品を手に入れた後に、 \(j\) が製品を購入する確率を表します。重み \(w_{ij}\) は原則として \(1\) よりも大きくなる(または負の値になる)ことがあり、これは個々の \(j\) が複数の製品を購入する場合に対応しています。総購入確率を最大化することは、将来の総利益を最大化することに対応します。利益率が最初の無料サンプルのコストよりも大きくなる場合、この戦略は好都合なものとなります。このモデル拡張ではノード自体に重みが付けられます。これは、私たちのマーケティングモデルでは、製品の無料サンプルを提供された人が将来再び購入する可能性と見なすことができます。モデルのこの追加情報により、最大化する目的関数は次のようになります。

\[C(\textbf{x}) = \sum_{i,j} w_{ij} x_i (1-x_j)+\sum_i w_i x_i.\]

量子コンピューターにおいてこの問題の解決策を見つけるには、まずそれをイジング・ハミルトニアンにマッピングする必要があります。これは、 \(x_i\rightarrow (1-Z_i)/2\) の代入により実行できます。ここで、 \(Z_i\) は、固有値 \(\pm 1\) を持つパウリZ演算子です。 これを行うと、

\[C(\textbf{Z}) = \sum_{i,j} \frac{w_{ij}}{4} (1-Z_i)(1+Z_j) + \sum_i \frac{w_i}{2} (1-Z_i) = -\frac{1}{2}\left( \sum_{i<j} w_{ij} Z_i Z_j +\sum_i w_i Z_i\right)+\mathrm{const},\]

ここで、 \(\mathrm{const} = \sum_{i<j}w_{ij}/2+\sum_i w_i/2\) です。言い換えると、重み付き最大カット問題は、イジング・ハミルトニアンの最小化と等価です。

\[H = \sum_i w_i Z_i + \sum_{i<j} w_{ij} Z_iZ_j.\]

Qiskitの最適化モジュールは、最初の利益関数 \(\tilde{C}\) のイジングハミルトニアンを生成できます。 この点で、関数 \(\tilde{C}\) は、 to_ising() メソッドを提供する QuadraticProgram としてモデル化できます。

最適化問題のための近似ユニバーサル量子コンピューティング#

量子コンピューターの使用については、組合せ最適化問題の解決策を見つけるため、近年かなりの関心が向けられています。組合せ問題の古典的な性質を考えると、最良の古典的アルゴリズムとの比較において、量子コンピューターを使用する指数関数的な加速は保証されないということが重要と言えます。しかしながら、対象とする問題の性質と重要性によっては、問題のインスタンスを高速化し得る量子コンピューターでのヒューリスティックなアプローチを調査することは価値をみいだせます。ここでは、Farhi、 Goldstone、および Gutmann (2014) による Quantum Approximate Optimization Algorithm (QAOA) に基づくアプローチを示します。ヒューリスティックな性質を前提として、アルゴリズムを 近似量子コンピューティング のコンテキストで構成します。

アルゴリズムは次のように機能します。

  1. ターゲット・イジング問題で \(w_i\)\(w_{ij}\) を選択します。 原則として、Zのさらに高い累乗が許可されます。

  2. 量子回路の深さ \(m\) を選択します。 深さは適応的に変更できることに注意してください。

  3. コントロール \(\theta\) のセットを選択し、 \(\theta\) の成分によってパラメータ化された制御位相ゲートと単一量子ビットY回転で構成される量子回路を使用して構築された試行関数 \(|\psi(\boldsymbol\theta)\rangle\) を作成します。

  4. 評価します。

    \[C(\boldsymbol\theta) = \langle\psi(\boldsymbol\theta)~|H|~\psi(\boldsymbol\theta)\rangle = \sum_i w_i \langle\psi(\boldsymbol\theta)~|Z_i|~\psi(\boldsymbol\theta)\rangle+ \sum_{i<j} w_{ij} \langle\psi(\boldsymbol\theta)~|Z_iZ_j|~\psi(\boldsymbol\theta)\rangle\]

    Z基底で回路の結果をサンプリングし、個々のイジング項の期待値を合計して、 上記を評価します。 一般に、選択した古典的なオプティマイザーに応じて、 \(\boldsymbol\theta\) 周辺のさまざまな制御点を推定する必要があります。

  5. 古典的なオプティマイザーを使用して、新しいコントロールのセットを選択します。

  6. \(C(\boldsymbol\theta)\) が最小値に達するまで続け、解 \(\boldsymbol\theta^*\) に十分に近づけます。

  7. 最後の \(\boldsymbol\theta\) を使用して、分布 \(|\langle z_i~|\psi(\boldsymbol\theta)\rangle|^2\;\forall i\) からサンプルの最終セットを生成し、答えを取得します。

優れたヒューリスティック・アルゴリズムを見つけることの難しさを、適切な試行波動関数の選択に帰着することが私たちの信念です。たとえば、エンタングルメントがターゲットの問題と最もよく一致する試行関数を検討するか、単にエンタングルメントの量を変数にすることができます。このチュートリアルでは、次の形式の簡単な試行関数を検討します

\[|\psi(\theta)\rangle = [U_\mathrm{single}(\boldsymbol\theta) U_\mathrm{entangler}]^m |+\rangle\]

ここで、 \(U_\mathrm{entangler}\) は制御位相ゲート(完全に絡み合うゲート)のコレクションであり、 \(U_\mathrm{single}(\theta) = \prod_{i=1}^n Y(\theta_{i})\) です。ここで \(n\) は量子ビットの数、 \(m\) は量子回路の深さです。この選択の動機は、これらの古典的な問題の場合、この選択により、実際の係数のみを持つ量子状態の空間を検索でき、エンタングルメントを利用して、より速く解に収束する可能性があるということです。

断熱アプローチと比較して、このサンプリング方法を使用する利点の1つは、ターゲットのイジング・ハミルトニアンをハードウェアに直接実装する必要がないため、このアルゴリズムをデバイスの接続に限定する必要がないことです。さらに、 \(Z_iZ_jZ_k\) などのコスト関数の高次項も効率的にサンプリングできますが、断熱的アプローチまたはアニーリング・アプローチでは、それらを処理することは一般的に非現実的です。

参考文献:

  •  A. Lucas, Frontiers in Physics 2, 5 (2014)

  •  E. Farhi, J. Goldstone, S. Gutmann, e-print arXiv 1411.4028 (2014)

  •  D. Wecker, M. B. Hastings, M. Troyer, Phys. Rev. A 94, 022309 (2016)

  •  E. Farhi, J. Goldstone, S. Gutmann, H. Neven, e-print arXiv 1703.06199 (2017)

アプリケーション・クラス#

このページでは、最大カット問題と巡回セールスマンの問題のためにアプリケーション・クラスを使用します。 他の最適化問題用のアプリケーション・クラスもそれぞれあります。詳細については、「最適化問題のアプリケーション・クラス 」を参照してください。

[1]:
# useful additional packages
import matplotlib.pyplot as plt
import numpy as np
import networkx as nx

from qiskit.tools.visualization import plot_histogram
from qiskit.circuit.library import TwoLocal
from qiskit_optimization.applications import Maxcut, Tsp
from qiskit_algorithms import SamplingVQE, NumPyMinimumEigensolver
from qiskit_algorithms.optimizers import SPSA
from qiskit_algorithms.utils import algorithm_globals
from qiskit.primitives import Sampler
from qiskit_optimization.algorithms import MinimumEigenOptimizer

最大カット問題#

[2]:
# Generating a graph of 4 nodes

n = 4  # Number of nodes in graph
G = nx.Graph()
G.add_nodes_from(np.arange(0, n, 1))
elist = [(0, 1, 1.0), (0, 2, 1.0), (0, 3, 1.0), (1, 2, 1.0), (2, 3, 1.0)]
# tuple is (i,j,weight) where (i,j) is the edge
G.add_weighted_edges_from(elist)

colors = ["r" for node in G.nodes()]
pos = nx.spring_layout(G)


def draw_graph(G, colors, pos):
    default_axes = plt.axes(frameon=True)
    nx.draw_networkx(G, node_color=colors, node_size=600, alpha=0.8, ax=default_axes, pos=pos)
    edge_labels = nx.get_edge_attributes(G, "weight")
    nx.draw_networkx_edge_labels(G, pos=pos, edge_labels=edge_labels)


draw_graph(G, colors, pos)
../_images/tutorials_06_examples_max_cut_and_tsp_5_0.png
[3]:
# Computing the weight matrix from the random graph
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)
[[0. 1. 1. 1.]
 [1. 0. 1. 0.]
 [1. 1. 0. 1.]
 [1. 0. 1. 0.]]

力まかせのアプローチ#

可能なすべての \(2^n\) の組み合わせを試してください。 \(n = 4\) の場合、この例のように16の組み合わせのみを扱いますが、n = 1000の場合、1.071509e+30の組み合わせを持っているため、力まかせのアプローチを使用して扱うのは非現実的です。

[4]:
best_cost_brute = 0
for b in range(2**n):
    x = [int(t) for t in reversed(list(bin(b)[2:].zfill(n)))]
    cost = 0
    for i in range(n):
        for j in range(n):
            cost = cost + w[i, j] * x[i] * (1 - x[j])
    if best_cost_brute < cost:
        best_cost_brute = cost
        xbest_brute = x
    print("case = " + str(x) + " cost = " + str(cost))

colors = ["r" if xbest_brute[i] == 0 else "c" for i in range(n)]
draw_graph(G, colors, pos)
print("\nBest solution = " + str(xbest_brute) + " cost = " + str(best_cost_brute))
case = [0, 0, 0, 0] cost = 0.0
case = [1, 0, 0, 0] cost = 3.0
case = [0, 1, 0, 0] cost = 2.0
case = [1, 1, 0, 0] cost = 3.0
case = [0, 0, 1, 0] cost = 3.0
case = [1, 0, 1, 0] cost = 4.0
case = [0, 1, 1, 0] cost = 3.0
case = [1, 1, 1, 0] cost = 2.0
case = [0, 0, 0, 1] cost = 2.0
case = [1, 0, 0, 1] cost = 3.0
case = [0, 1, 0, 1] cost = 4.0
case = [1, 1, 0, 1] cost = 3.0
case = [0, 0, 1, 1] cost = 3.0
case = [1, 0, 1, 1] cost = 2.0
case = [0, 1, 1, 1] cost = 3.0
case = [1, 1, 1, 1] cost = 0.0

Best solution = [1, 0, 1, 0] cost = 4.0
../_images/tutorials_06_examples_max_cut_and_tsp_8_1.png

イジング問題へのマッピング#

Qiskit optimizationは、問題の仕様から該当するイジング・ハミルトニアンを作成するのと同様に QuadraticProgram を生成する機能を提供します。

[5]:
max_cut = Maxcut(w)
qp = max_cut.to_quadratic_program()
print(qp.prettyprint())
Problem name: Max-cut

Maximize
  -2*x_0*x_1 - 2*x_0*x_2 - 2*x_0*x_3 - 2*x_1*x_2 - 2*x_2*x_3 + 3*x_0 + 2*x_1
  + 3*x_2 + 2*x_3

Subject to
  No constraints

  Binary variables (4)
    x_0 x_1 x_2 x_3

[6]:
qubitOp, offset = qp.to_ising()
print("Offset:", offset)
print("Ising Hamiltonian:")
print(str(qubitOp))
Offset: -2.5
Ising Hamiltonian:
SparsePauliOp(['IIZZ', 'IZIZ', 'IZZI', 'ZIIZ', 'ZZII'],
              coeffs=[0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j])
[7]:
# solving Quadratic Program using exact classical eigensolver
exact = MinimumEigenOptimizer(NumPyMinimumEigensolver())
result = exact.solve(qp)
print(result.prettyprint())
objective function value: 4.0
variable values: x_0=1.0, x_1=0.0, x_2=1.0, x_3=0.0
status: SUCCESS

問題は最小化問題に変換されたため、 \(-4\) の解が最適と一致します。

すべてのハミルトニアンが適切なコストを与えることの確認#

[8]:
# Making the Hamiltonian in its full form and getting the lowest eigenvalue and eigenvector
ee = NumPyMinimumEigensolver()
result = ee.compute_minimum_eigenvalue(qubitOp)

x = max_cut.sample_most_likely(result.eigenstate)
print("energy:", result.eigenvalue.real)
print("max-cut objective:", result.eigenvalue.real + offset)
print("solution:", x)
print("solution objective:", qp.objective.evaluate(x))

colors = ["r" if x[i] == 0 else "c" for i in range(n)]
draw_graph(G, colors, pos)
energy: -1.5
max-cut objective: -4.0
solution: [1. 0. 1. 0.]
solution objective: 4.0
../_images/tutorials_06_examples_max_cut_and_tsp_16_1.png

量子コンピューターでの実行#

Y量子ビット回転、 \(U_\mathrm{single}(\theta) = \prod_{i=1}^n Y(\theta_{i})\) および エンタングラー・ステップ \(U_\mathrm{entangler}\) で構築された試行関数を使用する量子コンピューターで、フィードバックループを使用して最適化ルーチンを実行します。

[9]:
algorithm_globals.random_seed = 123
seed = 10598
[10]:
# construct SamplingVQE
optimizer = SPSA(maxiter=300)
ry = TwoLocal(qubitOp.num_qubits, "ry", "cz", reps=5, entanglement="linear")
vqe = SamplingVQE(sampler=Sampler(), ansatz=ry, optimizer=optimizer)

# run SamplingVQE
result = vqe.compute_minimum_eigenvalue(qubitOp)

# print results
x = max_cut.sample_most_likely(result.eigenstate)
print("energy:", result.eigenvalue.real)
print("time:", result.optimizer_time)
print("max-cut objective:", result.eigenvalue.real + offset)
print("solution:", x)
print("solution objective:", qp.objective.evaluate(x))

# plot results
colors = ["r" if x[i] == 0 else "c" for i in range(n)]
draw_graph(G, colors, pos)
energy: -1.4996861455587311
time: 8.708028078079224
max-cut objective: -3.999686145558731
solution: [0 1 0 1]
solution objective: 4.0
../_images/tutorials_06_examples_max_cut_and_tsp_19_1.png
[11]:
# create minimum eigen optimizer based on SamplingVQE
vqe_optimizer = MinimumEigenOptimizer(vqe)

# solve quadratic program
result = vqe_optimizer.solve(qp)
print(result.prettyprint())

colors = ["r" if result.x[i] == 0 else "c" for i in range(n)]
draw_graph(G, colors, pos)
objective function value: 4.0
variable values: x_0=1.0, x_1=0.0, x_2=1.0, x_3=0.0
status: SUCCESS
../_images/tutorials_06_examples_max_cut_and_tsp_20_1.png

巡回セールスマン問題#

2世紀以上にわたり、コンピューター科学者や数学者の注意を引いた悪名高いNP完全問題であることに加え、巡回セールスマン問題(Traveling Salesman Problem、TSP)という名前が示唆するように、財務とマーケティングに重要な影響を与えます。通俗的に言えば、巡回セールスマンは商品を販売するために都市から都市へ行く人のことです。この場合の目的は、セールスマンがすべての都市を訪れ、その故郷、つまり旅行を始めた都市に戻ることができる最短の道を見つけることです。これを行うことにより、セールスマンは最短時間で潜在的な売上を最大化することができます。

この問題の重要性は、その “難易度” と、実際に発生する他の関連する組合せ最適化問題との普遍的な等価性に起因します。

19世紀初頭に、W.R.ハミルトンによっていくつかの初期分析を伴う数学的定式化が提案されました。数学的には、問題は最大カットの場合のように、グラフの観点から最も抽象化されています。グラフのノードのTSPは、各ノードを通過できる最短の ハミルトニアン・サイクル を要求します。ハミルトン・サイクルは、グラフのすべての頂点を一度使用する閉じた経路です。一般的な解決策は不明であり、それを効率的に(たとえば、多項式時間で)見つけるアルゴリズムは存在しないと予想されています。

\(n=|V|\) ノードと、距離 \(w_{ij}\) (頂点 \(i\) から頂点 \(j\) までの距離)におけるグラフ \(G=(V,E)\) 最短のハミルトニアンサイクルを見つけます。ハミルトニアンサイクルは、 \(N^2\) 変数 \(x_{i,p}\) によって記述されます。ここで、 \(i\) はノードを表し、 \(p\) はプロスペクティブ内のその順序を表します。 時間変数 \(p\) 、ノード \(i\) で解が発生した場合、決定変数は値1をとります。すべてのノードがサイクルに一度だけ出現し、ノードが発生するたびに出現する必要があります。これは、2つの制約に相当します(これ以降、指定されていない場合、加和は常に 0、1、… N-1 に渡ります)

\[\sum_{i} x_{i,p} = 1 ~~\forall p\]
\[\sum_{p} x_{i,p} = 1 ~~\forall i.\]

予想される順序のノードに対し、 \(x_{i,p}\)\(x_{j,p+1}\) が両方とも1の場合、 \((i,j) \notin E\) (グラフで接続されていない)であれば、エネルギー・ペナルティが発生するはずです。 このペナルティの形式は以下の通りです。

\[\sum_{i,j\notin E}\sum_{p} x_{i,p}x_{j,p+1}>0,\]

ここで、ハミルトニアン・サイクル \((p=N)\equiv (p=0)\) の境界条件が仮定されています。 ただし、ここでは完全に接続されたグラフが想定され、この項は含まれていません。 最小化する必要がある距離は以下の通りです。

\[C(\textbf{x})=\sum_{i,j}w_{ij}\sum_{p} x_{i,p}x_{j,p+1}.\]

これらすべてを1つの目標関数にまとめて最小化すると、次のようになります。

\[C(\textbf{x})=\sum_{i,j}w_{ij}\sum_{p} x_{i,p}x_{j,p+1}+ A\sum_p\left(1- \sum_i x_{i,p}\right)^2+A\sum_i\left(1- \sum_p x_{i,p}\right)^2,\]

ここで、 \(A\) は自由パラメーターです。 これらの制約が尊重されるように、 \(A\) が十分に大きいことを確認する必要があります。 これを行う1つの方法は、 \(A > \mathrm{max}(w_{ij})\) となるように \(A\) を選択することです。

繰り返しになりますが、この形式の問題を量子コンピューターにマッピングするのは簡単であり、解決策はイジング・ハミルトニアンを最小化することによって見つけられます。

[12]:
# Generating a graph of 3 nodes
n = 3
num_qubits = n**2
tsp = Tsp.create_random_instance(n, seed=123)
adj_matrix = nx.to_numpy_array(tsp.graph)
print("distance\n", adj_matrix)

colors = ["r" for node in tsp.graph.nodes]
pos = [tsp.graph.nodes[node]["pos"] for node in tsp.graph.nodes]
draw_graph(tsp.graph, colors, pos)
distance
 [[ 0. 48. 91.]
 [48.  0. 63.]
 [91. 63.  0.]]
../_images/tutorials_06_examples_max_cut_and_tsp_22_1.png

力まかせのアプローチ#

[13]:
from itertools import permutations


def brute_force_tsp(w, N):
    a = list(permutations(range(1, N)))
    last_best_distance = 1e10
    for i in a:
        distance = 0
        pre_j = 0
        for j in i:
            distance = distance + w[j, pre_j]
            pre_j = j
        distance = distance + w[pre_j, 0]
        order = (0,) + i
        if distance < last_best_distance:
            best_order = order
            last_best_distance = distance
            print("order = " + str(order) + " Distance = " + str(distance))
    return last_best_distance, best_order


best_distance, best_order = brute_force_tsp(adj_matrix, n)
print(
    "Best order from brute force = "
    + str(best_order)
    + " with total distance = "
    + str(best_distance)
)


def draw_tsp_solution(G, order, colors, pos):
    G2 = nx.DiGraph()
    G2.add_nodes_from(G)
    n = len(order)
    for i in range(n):
        j = (i + 1) % n
        G2.add_edge(order[i], order[j], weight=G[order[i]][order[j]]["weight"])
    default_axes = plt.axes(frameon=True)
    nx.draw_networkx(
        G2, node_color=colors, edge_color="b", node_size=600, alpha=0.8, ax=default_axes, pos=pos
    )
    edge_labels = nx.get_edge_attributes(G2, "weight")
    nx.draw_networkx_edge_labels(G2, pos, font_color="b", edge_labels=edge_labels)


draw_tsp_solution(tsp.graph, best_order, colors, pos)
order = (0, 1, 2) Distance = 202.0
Best order from brute force = (0, 1, 2) with total distance = 202.0
../_images/tutorials_06_examples_max_cut_and_tsp_24_1.png

イジング問題へのマッピング#

[14]:
qp = tsp.to_quadratic_program()
print(qp.prettyprint())
Problem name: TSP

Minimize
  48*x_0_0*x_1_1 + 48*x_0_0*x_1_2 + 91*x_0_0*x_2_1 + 91*x_0_0*x_2_2
  + 48*x_0_1*x_1_0 + 48*x_0_1*x_1_2 + 91*x_0_1*x_2_0 + 91*x_0_1*x_2_2
  + 48*x_0_2*x_1_0 + 48*x_0_2*x_1_1 + 91*x_0_2*x_2_0 + 91*x_0_2*x_2_1
  + 63*x_1_0*x_2_1 + 63*x_1_0*x_2_2 + 63*x_1_1*x_2_0 + 63*x_1_1*x_2_2
  + 63*x_1_2*x_2_0 + 63*x_1_2*x_2_1

Subject to
  Linear constraints (6)
    x_0_0 + x_0_1 + x_0_2 == 1  'c0'
    x_1_0 + x_1_1 + x_1_2 == 1  'c1'
    x_2_0 + x_2_1 + x_2_2 == 1  'c2'
    x_0_0 + x_1_0 + x_2_0 == 1  'c3'
    x_0_1 + x_1_1 + x_2_1 == 1  'c4'
    x_0_2 + x_1_2 + x_2_2 == 1  'c5'

  Binary variables (9)
    x_0_0 x_0_1 x_0_2 x_1_0 x_1_1 x_1_2 x_2_0 x_2_1 x_2_2

[15]:
from qiskit_optimization.converters import QuadraticProgramToQubo

qp2qubo = QuadraticProgramToQubo()
qubo = qp2qubo.convert(qp)
qubitOp, offset = qubo.to_ising()
print("Offset:", offset)
print("Ising Hamiltonian:")
print(str(qubitOp))
Offset: 7581.0
Ising Hamiltonian:
SparsePauliOp(['IIIIIIIIZ', 'IIIIIIIZI', 'IIIIIIZII', 'IIIIIZIII', 'IIIIZIIII', 'IIIZIIIII', 'IIZIIIIII', 'IZIIIIIII', 'ZIIIIIIII', 'IIIIIIIZZ', 'IIIIIIZIZ', 'IIIIIIZZI', 'IIIIIZIIZ', 'IIIIIZIZI', 'IIIIIZZII', 'IIIIZIIIZ', 'IIIIZIIZI', 'IIIIZIZII', 'IIIIZZIII', 'IIIZIIIIZ', 'IIIZIIIZI', 'IIIZIIZII', 'IIIZIZIII', 'IIIZZIIII', 'IIZIIIIIZ', 'IIZIIIIZI', 'IIZIIIZII', 'IIZIIZIII', 'IIZIZIIII', 'IIZZIIIII', 'IZIIIIIIZ', 'IZIIIIIZI', 'IZIIIIZII', 'IZIIIZIII', 'IZIIZIIII', 'IZIZIIIII', 'IZZIIIIII', 'ZIIIIIIIZ', 'ZIIIIIIZI', 'ZIIIIIZII', 'ZIIIIZIII', 'ZIIIZIIII', 'ZIIZIIIII', 'ZIZIIIIII', 'ZZIIIIIII'],
              coeffs=[-1282.5 +0.j, -1282.5 +0.j, -1282.5 +0.j, -1268.5 +0.j, -1268.5 +0.j,
 -1268.5 +0.j, -1290.  +0.j, -1290.  +0.j, -1290.  +0.j,   606.5 +0.j,
   606.5 +0.j,   606.5 +0.j,   606.5 +0.j,    12.  +0.j,    12.  +0.j,
    12.  +0.j,   606.5 +0.j,    12.  +0.j,   606.5 +0.j,    12.  +0.j,
    12.  +0.j,   606.5 +0.j,   606.5 +0.j,   606.5 +0.j,   606.5 +0.j,
    22.75+0.j,    22.75+0.j,   606.5 +0.j,    15.75+0.j,    15.75+0.j,
    22.75+0.j,   606.5 +0.j,    22.75+0.j,    15.75+0.j,   606.5 +0.j,
    15.75+0.j,   606.5 +0.j,    22.75+0.j,    22.75+0.j,   606.5 +0.j,
    15.75+0.j,    15.75+0.j,   606.5 +0.j,   606.5 +0.j,   606.5 +0.j])
[16]:
result = exact.solve(qubo)
print(result.prettyprint())
objective function value: 202.0
variable values: x_0_0=1.0, x_0_1=0.0, x_0_2=0.0, x_1_0=0.0, x_1_1=1.0, x_1_2=0.0, x_2_0=0.0, x_2_1=0.0, x_2_2=1.0
status: SUCCESS

すべてのハミルトニアンが適切なコストを与えることの確認#

[17]:
# Making the Hamiltonian in its full form and getting the lowest eigenvalue and eigenvector
ee = NumPyMinimumEigensolver()
result = ee.compute_minimum_eigenvalue(qubitOp)

print("energy:", result.eigenvalue.real)
print("tsp objective:", result.eigenvalue.real + offset)
x = tsp.sample_most_likely(result.eigenstate)
print("feasible:", qubo.is_feasible(x))
z = tsp.interpret(x)
print("solution:", z)
print("solution objective:", tsp.tsp_value(z, adj_matrix))
draw_tsp_solution(tsp.graph, z, colors, pos)
energy: -7379.0
tsp objective: 202.0
feasible: True
solution: [0, 1, 2]
solution objective: 202.0
../_images/tutorials_06_examples_max_cut_and_tsp_30_1.png

量子コンピューターでの実行#

Y量子ビット回転、 \(U_\mathrm{single}(\theta) = \prod_{i=1}^n Y(\theta_{i})\) および エンタングラー・ステップ \(U_\mathrm{entangler}\) で構築された試行関数を使用する量子コンピューターで、フィードバックループを使用して最適化ルーチンを実行します。

[18]:
algorithm_globals.random_seed = 123
seed = 10598
[19]:
optimizer = SPSA(maxiter=300)
ry = TwoLocal(qubitOp.num_qubits, "ry", "cz", reps=5, entanglement="linear")
vqe = SamplingVQE(sampler=Sampler(), ansatz=ry, optimizer=optimizer)

result = vqe.compute_minimum_eigenvalue(qubitOp)

print("energy:", result.eigenvalue.real)
print("time:", result.optimizer_time)
x = tsp.sample_most_likely(result.eigenstate)
print("feasible:", qubo.is_feasible(x))
z = tsp.interpret(x)
print("solution:", z)
print("solution objective:", tsp.tsp_value(z, adj_matrix))
draw_tsp_solution(tsp.graph, z, colors, pos)
energy: -7326.024699521861
time: 28.526036977767944
feasible: True
solution: [1, 2, 0]
solution objective: 202.0
../_images/tutorials_06_examples_max_cut_and_tsp_33_1.png
[20]:
algorithm_globals.random_seed = 123
seed = 10598
[21]:
# create minimum eigen optimizer based on SamplingVQE
vqe_optimizer = MinimumEigenOptimizer(vqe)

# solve quadratic program
result = vqe_optimizer.solve(qp)
print(result.prettyprint())

z = tsp.interpret(x)
print("solution:", z)
print("solution objective:", tsp.tsp_value(z, adj_matrix))
draw_tsp_solution(tsp.graph, z, colors, pos)
objective function value: 202.0
variable values: x_0_0=0.0, x_0_1=0.0, x_0_2=1.0, x_1_0=1.0, x_1_1=0.0, x_1_2=0.0, x_2_0=0.0, x_2_1=1.0, x_2_2=0.0
status: SUCCESS
solution: [1, 2, 0]
solution objective: 202.0
../_images/tutorials_06_examples_max_cut_and_tsp_35_1.png
[22]:
import qiskit.tools.jupyter

%qiskit_version_table
%qiskit_copyright

Version Information

Qiskit SoftwareVersion
qiskit-terra0.25.0.dev0+1d844ec
qiskit-aer0.12.0
qiskit-ibmq-provider0.20.2
qiskit-nature0.7.0
qiskit-optimization0.6.0
System information
Python version3.10.11
Python compilerClang 14.0.0 (clang-1400.0.29.202)
Python buildmain, Apr 7 2023 07:31:31
OSDarwin
CPUs4
Memory (Gb)16.0
Thu May 18 16:57:58 2023 JST

This code is a part of Qiskit

© Copyright IBM 2017, 2023.

This code is licensed under the Apache License, Version 2.0. You may
obtain a copy of this license in the LICENSE.txt file in the root directory
of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.

Any modifications or derivative works of this code must retain this
copyright notice, and modified files need to carry a notice indicating
that they have been altered from the originals.

[ ]: