{ "cells": [ { "cell_type": "markdown", "id": "d2c31ae8", "metadata": {}, "source": [ "---\n", "title: \"最適化マッパー Qiskit アドオンを使用したウォームスタート型 QAOA\"\n", "description: \"qiskit-addon-opt-mapper パッケージを使用して、連続緩和問題の解を初期値として設定することで、最大カット問題における QAOA の収束性を向上させる。\"\n", "---\n", "\n", "{/* cspell:ignore prereqs Mareček rhobeg skyblue steelblue edgecolor fontsize Farhi */}\n", "\n", "\n", "\n", "# 最適化マッパー Qiskit アドオンを使用したウォームスタート型 QAOA\n", "\n", "*推定使用時間：Heron r3 で 9 分（注：これはあくまで目安です。 （実行時間は状況によって異なる場合があります。）*\n", "\n" ] }, { "cell_type": "markdown", "id": "8bf80006", "metadata": {}, "source": [ "\n", "\n", "## 学習成果\n", "\n", "このチュートリアルを修了すると、以下の内容を理解できるようになります：\n", "\n", "* 以下を用いて、最大切断問題を量子二次無制約二値最適化（QUBO）の定式化に写像する方法 `qiskit-addon-opt-mapper`\n", "* シミュレータ上で標準QAOAを実装・実行する方法\n", "* 二次計画法（QP）による緩和問題を計算し、ウォームスタート回路を構築することで、WS-QAOAを適用する方法\n", "* 標準的なQAOAとWS-QAOAのエネルギー収束性と解の質を比較する方法\n", "\n", "\n", "\n", "## 前提条件\n", "\n", "以下のトピックについて、あらかじめ理解しておくことをお勧めします：\n", "\n", "* [QAOAチュートリアル](/docs/tutorials/quantum-approximate-optimization-algorithm)\n", "* [QAOA上級チュートリアル](/docs/tutorials/advanced-techniques-for-qaoa)\n", "\n", "\n", "\n", "## 背景\n", "\n", "量子近似最適化アルゴリズム（QAOA）は、最大カットや一般的なQUBO定式化といった組み合わせ最適化問題を解くために設計された、量子・古典ハイブリッドアルゴリズムである。 Qiskit における QAOA の基礎的な概要については、 [QAOA チュートリアル](/docs/tutorials/quantum-approximate-optimization-algorithm)を参照してください。より高度な回路構築手法については、 [上級者向け QAOA チュートリアル](/docs/tutorials/advanced-techniques-for-qaoa)を参照してください。\n", "\n", "標準的なQAOAでは：\n", "\n", "* 初期状態は、一様重ね合わせ $|+\\rangle^{\\otimes n}$ である。\n", "* 変分パラメータはランダムに初期化されます。\n", "* 従来の最適化アルゴリズムは、コスト関数を最小化するパラメータを探索する。\n", "\n", "しかし、現実的な問題規模やノイズの多い量子ハードウェアの場合、ランダムな初期化を行うと、収束が遅くなったり、局所極小に陥ったり、最適化コストが増大したりする可能性がある。\n", "\n", "**ウォームスタートQAOA** （WS-QAOA）は、古典的な最適化の知見を量子回路に直接組み込むことで、この問題を改善する。 このチュートリアルでは、Egger、Mareček、およびWoernerが『 [*Warm-starting quantum optimization*](https://arxiv.org/abs/2009.10095) 』で紹介した手法に従っています。 重要なポイントは次の通りです：\n", "\n", "1. 元の二値問題の**連続緩和問題** （ $\\{0,1\\}^n$ の代わりに $[0,1]^n$ 上の二次計画問題）を解く。\n", "2. $Y$ -rotation angles $\\theta_i = 2\\arcsin(\\sqrt{c^*_i})$ を用いて、 **緩和解** $c^*_i \\in [0,1]$ をカスタム初期状態にエンコードし、量子ビット $i$ が、 $|1\\rangle$ を測定した際の確率が $c^*_i$ となる状態で開始するようにする。\n", "3. **標準の $X$ -mixerを**、ウォームスタート初期状態を基底状態とするカスタムミキサーに置き換え、アルゴリズムが古典解の近くから開始され、その近傍を探索できるようにする。\n", "\n", "正則化パラメータ $\\varepsilon \\in [0, 0.5]$ は、到達可能性の問題を回避するために、 $c^*_i$ を0および1から切り離します。 $|0\\rangle$ または $|1\\rangle$ で初期化された量子ビットは、コストハミルトニアンによって移動させることはできません。 $\\varepsilon = 0.5$ において、WS-QAOA は標準的な QAOA に完全に帰着する。\n", "\n", "問題のモデル化には パッケージが [`qiskit-addon-opt-mapper`](https://qiskit.github.io/qiskit-addon-opt-mapper/) 使用されます。このパッケージの [`Maxcut`](https://qiskit.github.io/qiskit-addon-opt-mapper/stubs/qiskit_addon_opt_mapper.applications.Maxcut.html) アプリケーションクラスは、グラフから直接 QUBO を構築し、そのコンバータとトランスレータは、その結果として得られる問題を量子ハミルトニアンにマッピングします。\n", "\n" ] }, { "cell_type": "markdown", "id": "55b94021", "metadata": {}, "source": [ "\n", "\n", "## 要件\n", "\n", "このチュートリアルを始める前に、以下のものがインストールされていることを確認してください：\n", "\n", "* Qiskit SDK v2.0 またはそれ以降、 [可視化](/docs/api/qiskit/visualization)機能をサポートしたもの\n", "* Qiskit Runtime v0.43 またはそれ以降 (`pip install qiskit-ibm-runtime`)\n", "* 最適化マッパー Qiskit アドオン (`pip install qiskit-addon-opt-mapper`)\n", "* SciPy (`pip install scipy`)\n", "* NetworkX (`pip install networkx`)\n", "\n" ] }, { "cell_type": "markdown", "id": "7db2e559", "metadata": {}, "source": [ "\n", "\n", "## セットアップ\n", "\n", "必要なライブラリをすべてインポートし、このチュートリアル全体で使用されるヘルパー関数を定義します。\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "id": "bc380c46", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import networkx as nx\n", "from scipy.optimize import minimize\n", "\n", "from qiskit.circuit import QuantumCircuit, ParameterVector\n", "from qiskit.circuit.library import qaoa_ansatz\n", "from qiskit.quantum_info import Statevector\n", "from qiskit.primitives import StatevectorEstimator, StatevectorSampler\n", "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n", "from qiskit_ibm_runtime import (\n", " QiskitRuntimeService,\n", " Session,\n", " EstimatorOptions,\n", " EstimatorV2 as Estimator,\n", " SamplerV2 as Sampler,\n", ")\n", "\n", "from qiskit_addon_opt_mapper.applications import Maxcut\n", "from qiskit_addon_opt_mapper.converters import OptimizationProblemToQubo\n", "from qiskit_addon_opt_mapper.translators import to_ising" ] }, { "cell_type": "markdown", "id": "431a5bd2-e6ed-471b-ad9e-c4edd27784a8", "metadata": {}, "source": [ "\n", "\n", "# 小規模シミュレータの例\n", "\n", "ここでは、重み付きグラフ上の小さな**最大切断**問題を手本として用います。 Max-cut問題：辺の重みが $w_{ij}$ であるグラフ $G=(V,E)$ が与えられたとき、カットを横切る辺の総重みを最大化するように、頂点を2つの集合 $S$ と $\\bar{S}$ に分割する方法を求めよ。\n", "\n", "QUBO最小化問題として、最大カットは次のように表すことができる：\n", "$\\min_{x \\in \\{0,1\\}^n} -\\sum_{(i,j) \\in E} w_{ij}(x_i + x_j - 2x_i x_j)$\n", "\n", "シミュレータ上で処理しやすくするため、4ノードのグラフを用いて処理を行います。\n", "\n" ] }, { "cell_type": "markdown", "id": "988ee237", "metadata": {}, "source": [ "\n", "\n", "### ステップ1：古典的な入力を量子問題に写像する\n", "\n", "`qiskit-addon-opt-mapper`我々は、グラフから直接QUBO形式を構築する『』のアプリケーションクラ `Maxcut` スを用いて、最大カット問題を定義する。 次に、これをQUBOに変換し、QAOAに適したイジング・ハミルトニアン（`SparsePauliOp`）へと変換する。 また、QUBOの連続緩和（二値制約 $x_i \\in \\{0,1\\}$ を $x_i \\in [0,1]$ に置き換える）を解き、ウォームスタートの初期点 $c^*$ を求める。\n", "\n" ] }, { "cell_type": "code", "execution_count": 2, "id": "step1-graph-code", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Define a 4-node weighted graph for the max-cut problem\n", "n_nodes = 4\n", "edges = [(0, 1, 1.0), (0, 2, 1.0), (1, 2, 1.0), (1, 3, 1.0), (2, 3, 1.0)]\n", "\n", "G = nx.Graph()\n", "G.add_nodes_from(range(n_nodes))\n", "G.add_weighted_edges_from(edges)\n", "\n", "pos = nx.spring_layout(G, seed=42)\n", "edge_labels = {(u, v): d[\"weight\"] for u, v, d in G.edges(data=True)}\n", "\n", "fig, ax = plt.subplots(figsize=(4, 3))\n", "nx.draw(G, pos, with_labels=True, node_color=\"lightblue\", ax=ax)\n", "nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, ax=ax)\n", "ax.set_title(\"Max-Cut graph\")\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "step1-graph-md", "metadata": {}, "source": [ "このグラフには5本の辺があります。 最適なマックスカットは、ノードを $S = \\{0, 3\\}$ と $\\bar{S} = \\{1, 2\\}$ （またはその補集合）に分割し、5本の辺のうち4本を切断するため、カット値は4となる。\n", "\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "step1-qubo-code", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Problem name: Max-cut\n", "\n", "Maximize\n", " -2*x_0*x_1 - 2*x_0*x_2 - 2*x_1*x_2 - 2*x_1*x_3 - 2*x_2*x_3 + 2*x_0 + 3*x_1\n", " + 3*x_2 + 2*x_3\n", "\n", "Subject to\n", " No constraints\n", "\n", " Binary variables (4)\n", " x_0 x_1 x_2 x_3\n", "\n" ] } ], "source": [ "# Build the max-cut problem directly from the NetworkX graph using the\n", "# Maxcut application class. Internally it constructs the QUBO\n", "# minimize -sum_{(i,j) in E} w_ij * (x_i + x_j - 2*x_i*x_j)\n", "# (each edge contributes -w to the linear terms and +2w to the quadratic\n", "# term), so we get the same OptimizationProblem without the boilerplate.\n", "maxcut = Maxcut(G)\n", "prob = maxcut.to_optimization_problem()\n", "print(prob.prettyprint())" ] }, { "cell_type": "markdown", "id": "step1-qubo-md", "metadata": {}, "source": [ "この `Maxcut` クラスはQUBOの構築をラップしているため、最大カットの目的関数を手動で展開する必要がありません。 表示された目的関数には、各変数の線形係数（その変数が個別にカットにどれだけ寄与するか）と、各交差項の二次係数（隣接する2つのノードを同じ側に配置した場合のペナルティ）が示されています。 が `to_optimization_problem()` 返す基底 `OptimizationProblem` オブジェクトは、バイナリ、整数、連続、およびスピン型の変数をサポートしており、次のステップで使用されるコンバータやトランスレータが期待するオブジェクトと同じものです。\n", "\n" ] }, { "cell_type": "code", "execution_count": 4, "id": "step1-ising-code", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Cost Hamiltonian H_C (4 qubits):\n", "SparsePauliOp(['IIZZ', 'IZIZ', 'IZZI', 'ZIZI', 'ZZII'],\n", " coeffs=[0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j])\n", "\n", "Offset (constant shift): -2.5\n", " QUBO value = Ising energy + offset\n" ] } ], "source": [ "# Convert the OptimizationProblem to a QUBO, then translate to an Ising Hamiltonian\n", "#\n", "# The substitution x_i = (1 - z_i)/2 maps binary variables to spin operators,\n", "# yielding a Hamiltonian H_C = sum_i h_i Z_i + sum_{i to find the ground state, which encodes the optimal cut.\n", "converter = OptimizationProblemToQubo()\n", "qubo = converter.convert(prob)\n", "\n", "cost_operator, offset = to_ising(qubo)\n", "n_qubits = cost_operator.num_qubits\n", "\n", "print(f\"Cost Hamiltonian H_C ({n_qubits} qubits):\")\n", "print(cost_operator)\n", "print(f\"\\nOffset (constant shift): {offset}\")\n", "print(\" QUBO value = Ising energy + offset\")" ] }, { "cell_type": "markdown", "id": "step1-ising-md", "metadata": {}, "source": [ "この `to_ising` 変換器は、 $H_C$ を表す `SparsePauliOp` と、 $\\text{QUBO value} = \\langle H_C \\rangle + \\text{offset}$ を満たすスカラー `offset` を返す。すべての重みが 1 であるこの最大カット問題において、すべての量子ビットについて $h_i = 0$ が成り立つ（ $x_i \\to z_i$ の置換を行うと、グラフは線形的に対称となる）。また、各辺は、結合強度 $+0.5$ を持つ $Z_i Z_j$ の結合を寄与する。 $H_C$ の最小固有値は、最大カットに対応する。\n", "\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "step1-qp-code", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "QP relaxation solution c* = [1. 0. 0. 1.]\n", "QP objective value = -4.0000\n" ] } ], "source": [ "# Solve the continuous (QP) relaxation to obtain the warm-start point c*\n", "#\n", "# The QP relaxation replaces the binary constraint x_i in {0,1} with x_i in [0,1]\n", "# and minimizes the same quadratic objective. Its solution c*_i gives the\n", "# probability that variable i should be 1 according to the classical relaxation.\n", "#\n", "# The max-cut QUBO has a non-convex quadratic matrix (negative eigenvalues),\n", "# so the relaxed problem has multiple local minima. A naive single start from\n", "# [0.5,...,0.5] converges to the symmetric saddle point c* = [0.5,...,0.5],\n", "# which carries no useful structural information about the problem.\n", "# Multi-start optimization is used to reliably find the global minimum.\n", "Q = qubo.objective.quadratic.to_array(symmetric=True)\n", "mu = qubo.objective.linear.to_array()\n", "\n", "\n", "def qp_objective(x_cont):\n", " \"\"\"Continuous relaxation of the QUBO objective.\"\"\"\n", " return x_cont @ Q @ x_cont + mu @ x_cont + qubo.objective.constant\n", "\n", "\n", "bounds = [(0.0, 1.0)] * n_qubits\n", "\n", "rng = np.random.default_rng(42)\n", "best_val = np.inf\n", "c_star = None\n", "for _ in range(200):\n", " x0 = rng.uniform(0.0, 1.0, n_qubits)\n", " result = minimize(qp_objective, x0, method=\"L-BFGS-B\", bounds=bounds)\n", " if result.fun < best_val:\n", " best_val = result.fun\n", " c_star = result.x\n", "\n", "print(f\"QP relaxation solution c* = {np.round(c_star, 4)}\")\n", "print(f\"QP objective value = {best_val:.4f}\")" ] }, { "cell_type": "markdown", "id": "step1-qp-md", "metadata": {}, "source": [ "マルチスタートソルバーは、 $c^* = [1, 0, 0, 1]$ （またはその補集合である $[0, 1, 1, 0]$ ）を求め、これが実際の最適な二進解となります。 この問題では、QP緩和がタイトであり、連続最小値が整数最適解と一致するため、緩和によって最良のカットが即座に特定される。 ステップ2で $\\varepsilon = 0.25$ を用いて正則化した後、この解はウォームスタートの初期状態にエンコードされます。\n", "\n" ] }, { "cell_type": "markdown", "id": "ac6f36e3", "metadata": {}, "source": [ "\n", "\n", "### ステップ 2：量子ハードウェアでの実行に向けて問題を最適化する\n", "\n", "2つのQAOA回路を構築し、QP解からウォームスタート角を算出する。\n", "\n", "**標準的なQAOAでは**、初期状態として一様な重ね合わせ $|+\\rangle^{\\otimes n}$ を用い、層ごとに $\\prod_i R_X(-2\\beta)$ として実装された標準的な $X$ 混合器 $H_M = -\\sum_i X_i$ を採用しています。\n", "\n", "[\\[1\\]](#Reference1) に記載**されているウォームスタート型 QAOA（WS-QAOA）** では、1クビットあたり2つの構造変更が行われる $i$ ：\n", "\n", "* **初期状態：**$R_Y(\\theta_i)|0\\rangle$、 $\\theta_i = 2\\arcsin(\\sqrt{c^*_i})$ であるため、 $|1\\rangle$ が観測される確率は $c^*_i$ となる。\n", "* **カスタムミキサー：**$R_Y(\\theta_i)\\, R_Z(-2\\beta)\\, R_Y(-\\theta_i)$。その基底状態は $R_Y(\\theta_i)|0\\rangle$ である。 これは、WS-QAOAが自身のミキサーの基底状態から開始することを意味しており、これは標準的なQAOAが $|+\\rangle$ および $X$ ミキサーを用いて満たすのと同じ性質である。\n", "\n", "`p=1`層に関する注記：（単一のQAOA層において）標準的なQAOAは、三角形を含むグラフ（このグラフには0-1-2の三角形が含まれている）において、最適エネルギーの約49％に解析的に制限される。 ウォームスタートは、解に関する事前知識を初期状態に直接エンコードすることで、この制限を回避する。\n", "\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "step2-angles-code", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Continuous relaxation c* = [1. 0. 0. 1.]\n", "After regularization = [0.75 0.25 0.25 0.75]\n", "Warm-start angles theta = [2.0944 1.0472 1.0472 2.0944] radians\n", "\n", "Angle interpretation:\n", " theta = 0 <-> c* = 0 (qubit points toward |0>)\n", " theta = pi/2 <-> c* = 0.5 (qubit in equal superposition, like |+>)\n", " theta = pi <-> c* = 1 (qubit points toward |1>)\n" ] } ], "source": [ "# Number of QAOA layers (each layer = one cost unitary + one mixer unitary)\n", "p = 1\n", "\n", "# Regularization: clip c* to [epsilon, 1-epsilon] so no qubit is initialized\n", "# in |0> or |1>, which would freeze it under the cost Hamiltonian.\n", "epsilon = 0.25\n", "\n", "c_clipped = np.clip(c_star, epsilon, 1 - epsilon)\n", "thetas = 2 * np.arcsin(np.sqrt(c_clipped))\n", "\n", "print(f\"Continuous relaxation c* = {np.round(c_star, 4)}\")\n", "print(f\"After regularization = {np.round(c_clipped, 4)}\")\n", "print(f\"Warm-start angles theta = {np.round(thetas, 4)} radians\")\n", "print()\n", "print(\"Angle interpretation:\")\n", "print(\" theta = 0 <-> c* = 0 (qubit points toward |0>)\")\n", "print(\n", " \" theta = pi/2 <-> c* = 0.5 (qubit in equal superposition, like |+>)\"\n", ")\n", "print(\" theta = pi <-> c* = 1 (qubit points toward |1>)\")" ] }, { "cell_type": "markdown", "id": "step2-angles-md", "metadata": {}, "source": [ "クリッピング後、 $c^* = 1$ は $1 - \\varepsilon = 0.75$ となり、 $c^* = 0$ は $\\varepsilon = 0.25$ となる。その結果生じる角度 $\\theta \\approx [2.09, 1.05, 1.05, 2.09]$ ラジアンにより、キュービット0と3は $|1\\rangle$ の方向へ、キュービット1と2は $|0\\rangle$ の方向へと強く回転し、これにより最適カットの構造が初期量子状態に直接エンコードされる。\n", "\n" ] }, { "cell_type": "code", "execution_count": 7, "id": "step2-circuit-builders-code", "metadata": {}, "outputs": [], "source": [ "def apply_cost_unitary(qc, cost_op, gamma):\n", " \"\"\"Apply exp(-i * gamma * H_C) to the circuit.\n", "\n", " Each Pauli term in H_C contributes a rotation gate:\n", " - Single-Z term h_i * Z_i -> RZ(2 * gamma * h_i) on qubit i\n", " - Two-Z term J_ij * Z_i Z_j -> CNOT, RZ(2 * gamma * J_ij), CNOT\n", " \"\"\"\n", " for pauli_term, coeff in zip(cost_op.paulis, cost_op.coeffs):\n", " indices = [\n", " j for j, q in enumerate(pauli_term.to_label()[::-1]) if q == \"Z\"\n", " ]\n", " if len(indices) == 1:\n", " qc.rz(2 * gamma * coeff.real, indices[0])\n", " elif len(indices) == 2:\n", " qc.cx(indices[0], indices[1])\n", " qc.rz(2 * gamma * coeff.real, indices[1])\n", " qc.cx(indices[0], indices[1])\n", "\n", "\n", "def build_ws_qaoa(cost_op, n_layers, n_qubits, thetas):\n", " \"\"\"WS-QAOA: warm-start initial state + custom per-qubit mixer.\n", "\n", " Per Egger et al. (2021) Eq. (1)-(2):\n", " Initial state per qubit i: R_Y(theta_i) |0>\n", " Mixer gate per qubit i: R_Y(theta_i) R_Z(-2*beta) R_Y(-theta_i)\n", " \"\"\"\n", " gammas = ParameterVector(\"γ\", n_layers)\n", " betas = ParameterVector(\"β\", n_layers)\n", " qc = QuantumCircuit(n_qubits)\n", " for i, theta in enumerate(thetas):\n", " qc.ry(theta, i) # warm-start initial state\n", " for k in range(n_layers):\n", " apply_cost_unitary(qc, cost_op, gammas[k])\n", " for i, theta in enumerate(thetas):\n", " qc.ry(theta, i)\n", " qc.rz(-2 * betas[k], i)\n", " qc.ry(-theta, i)\n", " return qc, gammas, betas\n", "\n", "\n", "# Standard QAOA via the Qiskit built-in helper:\n", "# qaoa_ansatz prepares |+>^n, then alternates exp(-i*gamma*H_C) with the\n", "# default X-mixer for `reps` layers. The returned circuit exposes the\n", "# variational parameters via std_qc.parameters.\n", "std_qc = qaoa_ansatz(cost_operator, reps=p)\n", "\n", "# WS-QAOA: keep the custom builder. The per-qubit mixer\n", "# R_Y(theta_i) R_Z(-2*beta) R_Y(-theta_i) is implemented as an explicit gate\n", "# sequence rather than as a SparsePauliOp, so we construct the circuit\n", "# directly to stay close to the Egger et al. (2021) formulation.\n", "ws_qc, ws_gammas, ws_betas = build_ws_qaoa(cost_operator, p, n_qubits, thetas)" ] }, { "cell_type": "markdown", "id": "step2-circuit-builders-md", "metadata": {}, "source": [ "[`qaoa_ansatz`](/docs/api/qiskit/qiskit.circuit.library.qaoa_ansatz)標準的なアプローチについては、 $|+\\rangle^{\\otimes n}$ を構築し、コストユニタリーを適用し、各レイヤーに対して `reps` デフォルトの $X$ -ミキサーを適用する に委ねます。 WS-QAOA については、クビットごとのミキサー $R_Y(\\theta)\\,R_Z(-2\\beta)\\,R_Y(-\\theta)$ が、パウリ演算の和ではなくゲートシーケンスとして表現されるため、明示 `build_ws_qaoa` 的なヘルパーを残しています。 この `apply_cost_unitary` ヘルパーはハミルトニアン `SparsePauliOp` から直接読み込むため、手動で回路を構築することなく、あらゆるQUBO問題を処理することができます。\n", "\n" ] }, { "cell_type": "code", "execution_count": 8, "id": "step2-draw-code", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Standard QAOA circuit (p=1):\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "print(\"Standard QAOA circuit (p=1):\")\n", "std_qc.draw(\"mpl\", fold=-1)" ] }, { "cell_type": "code", "execution_count": 9, "id": "6fad9eda", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "WS-QAOA circuit (p=1):\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "print(\"\\nWS-QAOA circuit (p=1):\")\n", "ws_qc.draw(\"mpl\", fold=-1)" ] }, { "cell_type": "markdown", "id": "step2-draw-md", "metadata": {}, "source": [ "どちらの回路も、同じ構造を採用しています。すなわち、初期状態準備層に続き、 $p$、コストユニタリー層とミキサーユニタリー層が交互に配置されています。 WS-QAOA回路では、冒頭の $R_Y$ ゲートが $c^*$ を符号化し、ミキサーは各 $R_X$ を、共役な $R_Y$ – $R_Z$ – $R_Y$ の3つ組に置き換えます。 これら2つの回路の深さの差は、 $p$ に比例して直線的に大きくなりますが、深さが浅い場合は許容範囲内に収まります。\n", "\n" ] }, { "cell_type": "markdown", "id": "b4d480b3", "metadata": {}, "source": [ "\n", "\n", "### ステップ 3: `Qiskit primitives` を使用して実行する\n", "\n", "我々は、正確でノイズのないシミュレーションを行うために を使用 `StatevectorEstimator` しています。 SciPy の COBYLA オプティマイザ `minimize` を用いた関数が、変分ループを駆動し、各反復ごとに推定関数を呼び出して、与えられたパラメータセット $(\\gamma, \\beta)$ に対する $\\langle H_C \\rangle$ を評価します。\n", "\n", "これら2つのアルゴリズムは、最適化前にそれぞれが持っている知識を反映した、異なる初期パラメータを使用しています：\n", "\n", "* **標準的なQAOA：**$[0, \\pi]$ におけるランダム初期化 — 構造情報が得られないため、これは妥当である。\n", "* **WS-QAOA:** $\\gamma = 0$, $\\beta = \\pi/4$ — $\\gamma=0$ によると、コスト単位は恒等写像であるため、最初の回路評価ではウォームスタートの初期状態から直接サンプリングが行われる。 これにより、COBYLAは従来の解と整合した強力な出発点を得ることになる。\n", "\n" ] }, { "cell_type": "code", "execution_count": 10, "id": "step3-optimize-code", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Standard QAOA optimal energy : -0.5859\n", " optimal params: [0.6803 2.0533]\n", " optimizer calls: 47\n", "\n", "WS-QAOA optimal energy : -1.5000\n", " optimal params: gamma=[-0.0001], beta=[1.5708]\n", " optimizer calls: 42\n" ] } ], "source": [ "estimator = StatevectorEstimator()\n", "\n", "\n", "def make_cost_fn(circuit, param_order, cost_op, estimator, history):\n", " \"\"\"Return a scalar cost function compatible with scipy.optimize.minimize.\"\"\"\n", "\n", " def cost_fn(params):\n", " bound = circuit.assign_parameters(dict(zip(param_order, params)))\n", " job = estimator.run([(bound, cost_op)])\n", " energy = job.result()[0].data.evs.real\n", " history.append(energy)\n", " return energy\n", "\n", " return cost_fn\n", "\n", "\n", "# Standard QAOA: random initialization\n", "np.random.seed(42)\n", "std_param_order = list(std_qc.parameters)\n", "std_params0 = np.random.uniform(0, np.pi, len(std_param_order))\n", "std_history = []\n", "\n", "std_result = minimize(\n", " make_cost_fn(\n", " std_qc, std_param_order, cost_operator, estimator, std_history\n", " ),\n", " std_params0,\n", " method=\"COBYLA\",\n", " options={\"maxiter\": 300, \"rhobeg\": 0.5},\n", ")\n", "print(f\"Standard QAOA optimal energy : {std_result.fun:.4f}\")\n", "print(f\" optimal params: {std_result.x.round(4)}\")\n", "print(f\" optimizer calls: {len(std_history)}\")\n", "\n", "\n", "# WS-QAOA: informed initialization\n", "ws_params0 = np.concatenate([np.zeros(p), np.full(p, np.pi / 4)])\n", "ws_history = []\n", "ws_param_order = list(ws_gammas) + list(ws_betas)\n", "\n", "ws_result = minimize(\n", " make_cost_fn(ws_qc, ws_param_order, cost_operator, estimator, ws_history),\n", " ws_params0,\n", " method=\"COBYLA\",\n", " options={\"maxiter\": 300, \"rhobeg\": 0.5},\n", ")\n", "print(f\"\\nWS-QAOA optimal energy : {ws_result.fun:.4f}\")\n", "print(\n", " f\" optimal params: gamma={ws_result.x[:p].round(4)}, beta={ws_result.x[p:].round(4)}\"\n", ")\n", "print(f\" optimizer calls: {len(ws_history)}\")" ] }, { "cell_type": "markdown", "id": "step3-optimize-md", "metadata": {}, "source": [ "WS-QAOAの「情報に基づく出発点」という特徴により、COBYLAはウォームスタート解に近い有意義なエネルギー値から計算を開始するのに対し、標準的なQAOAはエネルギーランドスケープ上の実質的にランダムな点から開始することになります。 この初期品質の差こそが、ステップ4で見られる収束ギャップの主な要因である。\n", "\n" ] }, { "cell_type": "code", "execution_count": 11, "id": "step3-reference-code", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Exact optimal energy : -1.5000\n", "Standard QAOA approx. ratio : 0.3906\n", "WS-QAOA approx. ratio : 1.0000\n" ] } ], "source": [ "# Compute the exact optimal energy by brute-force over all 2^n bitstrings\n", "all_energies = [\n", " Statevector.from_label(format(k, f\"0{n_qubits}b\"))\n", " .expectation_value(cost_operator)\n", " .real\n", " for k in range(2**n_qubits)\n", "]\n", "optimal_energy = min(all_energies)\n", "\n", "print(f\"Exact optimal energy : {optimal_energy:.4f}\")\n", "print(f\"Standard QAOA approx. ratio : {std_result.fun / optimal_energy:.4f}\")\n", "print(f\"WS-QAOA approx. ratio : {ws_result.fun / optimal_energy:.4f}\")" ] }, { "cell_type": "markdown", "id": "step3-reference-md", "metadata": {}, "source": [ "近似比は、 $\\langle H_C \\rangle_{\\text{QAOA}} / E_{\\text{opt}}$ と定義される。 $E_{\\text{opt}} < 0$ となる最小化問題において、この比が 1 に近いほど、アルゴリズムがより低いエネルギー（より良い解）を見つけたことを意味する。 $2^n$ のすべての基底状態に対する総当たり探索は、 $n$ が小さい場合にのみ実行可能であり、真の値の参照として機能する。\n", "\n" ] }, { "cell_type": "markdown", "id": "50b94af2", "metadata": {}, "source": [ "\n", "\n", "### ステップ4：後処理を行い、結果を所望の従来の形式で出力する\n", "\n", "収束状況を可視化し、ビットストリング解について最適化された回路をサンプリングし、それらのビットストリングをデコードして最大切断分割に戻し、最終結果をまとめます。\n", "\n" ] }, { "cell_type": "code", "execution_count": 12, "id": "step4-convergence-code", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(figsize=(7, 4))\n", "ax.plot(std_history, label=\"Standard QAOA\", alpha=0.85)\n", "ax.plot(ws_history, label=\"WS-QAOA\", alpha=0.85)\n", "ax.axhline(\n", " optimal_energy,\n", " color=\"k\",\n", " linestyle=\"--\",\n", " label=f\"Exact optimal ({optimal_energy:.2f})\",\n", ")\n", "ax.set_xlabel(\"Optimizer call\")\n", "ax.set_ylabel(r\"$\\langle H_C \\rangle$\")\n", "ax.set_title(\"Convergence: Standard QAOA vs. WS-QAOA\")\n", "ax.legend()\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "step4-convergence-md", "metadata": {}, "source": [ "収束プロットには、COBYLA関数の各評価時点におけるエネルギー $\\langle H_C \\rangle$ が示されています。 $p=1$ における標準的なQAOAは、このグラフ上で最適エネルギーの約49％（三角形を含むグラフにおける $p=1$ 型QAOAの理論上の最大値）に制限されており、 $-0.74$ 付近で収束する。一方、最適解の近くで初期化されたWS-QAOAは、はるかに少ない反復回数で、 $-1.50$ （厳密な最適値）付近に素早く収束する。 これは、ウォームスタートの主な利点を示しています。すなわち、同じ回路の深さにおいて、ウォームスタートの方がはるかに優れた解に到達するのです。\n", "\n" ] }, { "cell_type": "code", "execution_count": 13, "id": "step4-sample-code", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Standard QAOA most-probable bitstring : 0110\n", " Partition: S=[0, 3], S̄=[1, 2] | cut value = 4.0\n", "\n", "WS-QAOA most-probable bitstring : 0110\n", " Partition: S=[0, 3], S̄=[1, 2] | cut value = 4.0\n" ] } ], "source": [ "# Sample the optimized circuits to recover the most probable bitstring solutions\n", "sampler = StatevectorSampler()\n", "shots = 1024\n", "\n", "\n", "def get_best_bitstring(circuit, param_order, optimal_params, sampler, shots):\n", " bound = circuit.assign_parameters(dict(zip(param_order, optimal_params)))\n", " bound.measure_all()\n", " job = sampler.run([bound], shots=shots)\n", " counts = job.result()[0].data.meas.get_counts()\n", " return max(counts, key=counts.get), counts\n", "\n", "\n", "def evaluate_cut(bitstring, G):\n", " \"\"\"Compute the Max-Cut value for a bitstring node assignment.\"\"\"\n", " x = [int(b) for b in bitstring]\n", " cut_val = sum(\n", " w for u, v, w in G.edges.data(\"weight\", default=1) if x[u] != x[v]\n", " )\n", " set0 = [i for i, b in enumerate(bitstring) if b == \"0\"]\n", " set1 = [i for i, b in enumerate(bitstring) if b == \"1\"]\n", " return cut_val, set0, set1\n", "\n", "\n", "# Qiskit bitstring ordering: rightmost character = qubit 0\n", "def decode_bitstring(bs):\n", " return bs[::-1]\n", "\n", "\n", "std_best, std_counts = get_best_bitstring(\n", " std_qc, std_param_order, std_result.x, sampler, shots\n", ")\n", "ws_best, ws_counts = get_best_bitstring(\n", " ws_qc, ws_param_order, ws_result.x, sampler, shots\n", ")\n", "\n", "std_cut, std_s0, std_s1 = evaluate_cut(decode_bitstring(std_best), G)\n", "ws_cut, ws_s0, ws_s1 = evaluate_cut(decode_bitstring(ws_best), G)\n", "\n", "print(f\"Standard QAOA most-probable bitstring : {std_best}\")\n", "print(f\" Partition: S={std_s0}, S̄={std_s1} | cut value = {std_cut}\")\n", "print()\n", "print(f\"WS-QAOA most-probable bitstring : {ws_best}\")\n", "print(f\" Partition: S={ws_s0}, S̄={ws_s1} | cut value = {ws_cut}\")" ] }, { "cell_type": "markdown", "id": "step4-sample-md", "metadata": {}, "source": [ "からの `Sampler` ビット文字列は、最右位置にクビット 0 が配置された状態で返されるため、この文字列を逆順にすると、インデックス $i$ が変数 $x_i$ に割り当てられる。カット値とは、パーティションを横切る辺の総重みであり、これが最大カット問題で最大化を目指す値である。 カット値が4の場合、利用可能な5つの辺のうち4つが使用され、これはこのグラフにおける理論上の最大値である。\n", "\n" ] }, { "cell_type": "code", "execution_count": 14, "id": "step4-visualize-code", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "=== Summary ===\n", "Method Ising energy Cut value Approx. ratio\n", "-----------------------------------------------------------------\n", "Standard QAOA -0.5859 4.0 0.3906\n", "WS-QAOA -1.5000 4.0 1.0000\n", "Exact optimal -1.5000 4 1.0000\n" ] } ], "source": [ "# Visualize the WS-QAOA solution on the graph\n", "fig, axes = plt.subplots(1, 2, figsize=(8, 3))\n", "\n", "for ax, s0, s1, cut, title in [\n", " (axes[0], std_s0, std_s1, std_cut, f\"Standard QAOA (cut = {std_cut})\"),\n", " (axes[1], ws_s0, ws_s1, ws_cut, f\"WS-QAOA (cut = {ws_cut})\"),\n", "]:\n", " colors = [\"skyblue\" if i in s0 else \"salmon\" for i in G.nodes()]\n", " nx.draw(G, pos, with_labels=True, node_color=colors, ax=ax)\n", " nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, ax=ax)\n", " ax.set_title(title)\n", "\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "# Summary\n", "# to_ising offset: QUBO value = Ising energy + offset, so Max-Cut value = -(Ising energy + offset)\n", "optimal_cut = -(optimal_energy + offset)\n", "print(\"=== Summary ===\")\n", "print(\n", " f\"{'Method':<20} {'Ising energy':>14} {'Cut value':>12} {'Approx. ratio':>15}\"\n", ")\n", "print(\"-\" * 65)\n", "print(\n", " f\"{'Standard QAOA':<20} {std_result.fun:>14.4f} {std_cut:>12} {std_result.fun/optimal_energy:>15.4f}\"\n", ")\n", "print(\n", " f\"{'WS-QAOA':<20} {ws_result.fun:>14.4f} {ws_cut:>12} {ws_result.fun/optimal_energy:>15.4f}\"\n", ")\n", "print(\n", " f\"{'Exact optimal':<20} {optimal_energy:>14.4f} {optimal_cut:>12.0f} {'1.0000':>15}\"\n", ")" ] }, { "cell_type": "markdown", "id": "step4-visualize-md", "metadata": {}, "source": [ "このグラフの可視化では、各ノードはパーティションの割り当てに応じて色分けされています（青＝ $S$、オレンジ＝ $\\bar{S}$ ）。パーティションをまたぐエッジ（異なる色のノードを結ぶエッジ）が、カットとしてカウントされます。\n", "\n", "どちらの方法も、カット値が4のビット文字列を見つけますが、その理由はまったく異なります。 **収束プロットとサンプリングされたビット列は、それぞれ異なるものを測定している**点に留意することが重要です：\n", "\n", "* **収束プロット**は、完全な量子状態の平均エネルギー $\\langle H_C \\rangle$ を追跡するもので、これは重ね合わせに含まれるすべてのビット列に対する重み付き平均である。 標準的なQAOAは、約 $-0.62$ に収束するが、これは最適値である $-1.50$ を大幅に上回っており、その量子状態は多くの次善のビット列に分散しており、正しい答えが含まれるのはごくまれであることを意味する。\n", "* **サンプリングされたビット列**は、その状態から1回抽出したものです。 この点で、標準的なQAOAは幸運だった。拡散状態からであっても、最適な分割がたまたま最も頻繁にサンプリングされる結果だったのだ。 難易度の高い問題や、ノイズの多いハードウェア、あるいは競合する候補解が多い場合、この「運」は尽きてしまう。\n", "\n", "対照的に、WS-QAOAは平均エネルギーを $-1.50$ まで収束させる。つまり、その量子状態は最適なビット列に集中しているということである。 ほぼすべての試行で正しい答えが得られるため、この解法は偶然ではなく、確実に正解を導き出すことができる。\n", "\n", "実際的な結果として、この小型でノイズのないシミュレータ上ではその違いは些細なものに見えるかもしれませんが、問題規模が大きくなったり、実際のハードウェア上で実行したりすると、平均エネルギーが最適値に近い状態は、拡散的な分布からたまに正しい答えをサンプリングするだけの状態よりも、はるかに堅牢です。\n", "\n" ] }, { "cell_type": "code", "execution_count": 15, "id": "b01696c2", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "P(cut = 4) | Standard QAOA = 0.4639 WS-QAOA = 1.0000\n" ] } ], "source": [ "# Compare the full probability distribution over cut values for both\n", "# algorithms. The most-probable bitstring above only reveals the mode;\n", "# this histogram exposes how much of the quantum state's probability mass\n", "# lands on the optimal cut versus on suboptimal partitions.\n", "def cut_value_distribution(counts, G, shots):\n", " dist = {}\n", " for bs, c in counts.items():\n", " cut, _, _ = evaluate_cut(decode_bitstring(bs), G)\n", " dist[cut] = dist.get(cut, 0.0) + c / shots\n", " return dist\n", "\n", "\n", "std_cut_dist = cut_value_distribution(std_counts, G, shots)\n", "ws_cut_dist = cut_value_distribution(ws_counts, G, shots)\n", "\n", "cut_values = sorted(set(std_cut_dist) | set(ws_cut_dist))\n", "std_probs = [std_cut_dist.get(c, 0.0) for c in cut_values]\n", "ws_probs = [ws_cut_dist.get(c, 0.0) for c in cut_values]\n", "\n", "fig, ax = plt.subplots(figsize=(7, 4))\n", "x = np.arange(len(cut_values))\n", "width = 0.4\n", "ax.bar(\n", " x - width / 2, std_probs, width, label=\"Standard QAOA\", color=\"steelblue\"\n", ")\n", "ax.bar(x + width / 2, ws_probs, width, label=\"WS-QAOA\", color=\"salmon\")\n", "ax.axvline(\n", " cut_values.index(optimal_cut),\n", " color=\"k\",\n", " linestyle=\"--\",\n", " alpha=0.4,\n", " label=f\"Optimal cut = {optimal_cut:g}\",\n", ")\n", "ax.set_xticks(x)\n", "ax.set_xticklabels([f\"{c:g}\" for c in cut_values])\n", "ax.set_xlabel(\"Cut value\")\n", "ax.set_ylabel(\"Probability\")\n", "ax.set_title(f\"Probability of measuring each cut value ({shots} shots)\")\n", "ax.legend()\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "print(\n", " f\"P(cut = {optimal_cut:g}) | Standard QAOA = \"\n", " f\"{std_cut_dist.get(optimal_cut, 0):.4f} \"\n", " f\"WS-QAOA = {ws_cut_dist.get(optimal_cut, 0):.4f}\"\n", ")" ] }, { "cell_type": "markdown", "id": "0763ae8e", "metadata": {}, "source": [ "このヒストグラムは、収束プロットが示唆していたことを数値的に表したものである。 標準的なQAOAでは、確率が複数の次善のカット値に分散しているため、1回の試行で4という最適なカットをサンプリングできる確率は、総質量のほんの一部に過ぎません。 WS-QAOAは、その確率のほぼすべてを最適カットに集中させているため、ほぼすべての試行で正しい答えが返されます。 これは、平均エネルギーが基底状態のエネルギーに収束した状態と、単に広い重ね合わせの中に基底状態が含まれているだけの状態とを区別する、実用的な特徴である。\n", "\n" ] }, { "cell_type": "markdown", "id": "0d6db390-e7a8-4efe-902c-8d9a312170c6", "metadata": {}, "source": [ "\n", "\n", "# 大規模ハードウェアの例\n", "\n" ] }, { "cell_type": "markdown", "id": "ae69c5e0-32b1-4f03-ab13-7b95a9acfd25", "metadata": {}, "source": [ "\n", "\n", "### 手順 1～4 を 1 つのコードブロックにまとめる\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "d58164ca-3b20-441c-9777-728497580cab", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using backend: ibm_boston\n" ] } ], "source": [ "# Selecting a backend using real hardware\n", "service = QiskitRuntimeService()\n", "backend = service.least_busy(\n", " operational=True, simulator=False, min_num_qubits=127\n", ")\n", "print(f\"Using backend: {backend.name}\")" ] }, { "cell_type": "code", "execution_count": 17, "id": "a35f3b21", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Graph: 40 nodes, 60 edges (3-regular)\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Cost operator: 40 qubits, 60 Pauli terms\n", "c* range: [0.000, 1.000] theta range: [1.047, 2.094] rad\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "Transpiled circuit: 2Q depth=86\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# ── Step 1a: Build the 40-node Max-Cut problem ─────────────────────────────\n", "# A 3-regular graph (every node has exactly 3 neighbors) is a standard QAOA\n", "N_LARGE = 40\n", "G_large = nx.random_regular_graph(d=3, n=N_LARGE, seed=0)\n", "edges_large = list(G_large.edges())\n", "print(f\"Graph: {N_LARGE} nodes, {len(edges_large)} edges (3-regular)\")\n", "\n", "# Visualize the graph so it is clear what problem we are solving before any\n", "# quantum work. Nodes in a circular layout; each edge contributes +1 to the\n", "# cut value when its endpoints land in different partitions.\n", "pos_large = nx.circular_layout(G_large)\n", "fig, ax = plt.subplots(figsize=(6, 6))\n", "nx.draw(\n", " G_large,\n", " pos_large,\n", " with_labels=True,\n", " node_color=\"lightblue\",\n", " node_size=400,\n", " font_size=7,\n", " ax=ax,\n", ")\n", "ax.set_title(f\"40-node 3-regular Max-Cut graph ({len(edges_large)} edges)\")\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "\n", "# Same Maxcut → OptimizationProblem → QUBO → Ising pipeline as the small example,\n", "# applied to the 40-node graph.\n", "prob_large = Maxcut(G_large).to_optimization_problem()\n", "converter_large = OptimizationProblemToQubo()\n", "qubo_large = converter_large.convert(prob_large)\n", "cost_op_large, offset_large = to_ising(qubo_large)\n", "n_qubits_large = cost_op_large.num_qubits\n", "print(\n", " f\"Cost operator: {n_qubits_large} qubits, {len(cost_op_large)} Pauli terms\"\n", ")\n", "\n", "# ── Step 1b: QP relaxation (multi-start L-BFGS-B) ─────────────────────────\n", "# Same multi-start approach as the small example. At 40 qubits the relaxed\n", "# landscape has many more local minima, so 200 random starts are essential\n", "# to find a low-energy warm-start point.\n", "Q_large = qubo_large.objective.quadratic.to_array(symmetric=True)\n", "mu_large = qubo_large.objective.linear.to_array()\n", "\n", "\n", "def qp_obj_large(x):\n", " return x @ Q_large @ x + mu_large @ x + qubo_large.objective.constant\n", "\n", "\n", "bounds_large = [(0.0, 1.0)] * n_qubits_large\n", "rng_qp = np.random.default_rng(42)\n", "best_val_large, c_star_large = np.inf, None\n", "\n", "for _ in range(200):\n", " x0 = rng_qp.uniform(0.0, 1.0, n_qubits_large)\n", " res = minimize(qp_obj_large, x0, method=\"L-BFGS-B\", bounds=bounds_large)\n", " if res.fun < best_val_large:\n", " best_val_large, c_star_large = res.fun, res.x\n", "\n", "# Regularize and convert to rotation angles (same formula as small example)\n", "epsilon_large = 0.25\n", "c_clipped_large = np.clip(c_star_large, epsilon_large, 1 - epsilon_large)\n", "thetas_large = 2 * np.arcsin(np.sqrt(c_clipped_large))\n", "print(\n", " f\"c* range: [{c_star_large.min():.3f}, {c_star_large.max():.3f}] \"\n", " f\"theta range: [{thetas_large.min():.3f}, {thetas_large.max():.3f}] rad\"\n", ")\n", "\n", "# Plot the distribution of c* values to see how much structure the relaxation\n", "# extracted. Values near 0/1 mean confident assignments; values near 0.5 mean\n", "# the classical solver was uncertain and quantum exploration is most needed there.\n", "fig, ax = plt.subplots(figsize=(6, 3))\n", "ax.hist(c_star_large, bins=20, color=\"steelblue\", edgecolor=\"white\")\n", "ax.axvline(0.5, color=\"k\", linestyle=\"--\", label=\"Uniform prior (std QAOA)\")\n", "ax.set_xlabel(r\"$c^*_i$\")\n", "ax.set_ylabel(\"Count\")\n", "ax.set_title(r\"Distribution of warm-start values $c^*_i$ (40-node graph)\")\n", "ax.legend()\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "# ── Step 1c: Build WS-QAOA circuit ─────────────────────────────────────────\n", "# Reuse build_ws_qaoa from the small-scale section unchanged; the helper\n", "# scales automatically with n_qubits and the cost operator size.\n", "p_large = 1\n", "ws_qc_large, ws_gammas_large, ws_betas_large = build_ws_qaoa(\n", " cost_op_large, p_large, n_qubits_large, thetas_large\n", ")\n", "ws_qc_large.measure_all()\n", "\n", "# ── Step 2: Transpile to hardware-native gates ──────────────────────────\n", "# generate_preset_pass_manager compiles the abstract circuit to th\n", "# gate set of the backend and inserts SWAP gates wherever the cost Hamiltonian\n", "# couples qubits that are not directly connected on the processor.\n", "pm = generate_preset_pass_manager(optimization_level=3, backend=backend)\n", "ws_isa_large = pm.run(ws_qc_large)\n", "\n", "ecr_count = ws_isa_large.count_ops().get(\"ecr\", 0)\n", "print(\n", " f\"\\nTranspiled circuit: 2Q depth={ws_isa_large.depth(lambda x: x.operation.num_qubits == 2)}\"\n", ")\n", "ws_isa_large.draw(\"mpl\", fold=-1)" ] }, { "cell_type": "code", "execution_count": 18, "id": "00ae1953", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Simulated annealing cut value: 53 (classical reference)\n", " iter 31 = -12.4094\n", "Optimization complete: energy=-13.0256, iterations=31\n", "Most-probable bitstring frequency: 4/8192 (0.0%)\n", "WS-QAOA cut: 53 | SA cut: 53 | Approximation ratio vs SA: 1.0000\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "=== Large Scale Summary ===\n", "Metric Value\n", "--------------------------------------------------\n", "Nodes / Edges 40 / 60 \n", "QAOA layers (p) 1\n", "Transpiled ECR gate count 0\n", "Transpiled circuit depth 276\n", "Optimizer iterations 31\n", "WS-QAOA energy (hardware) -13.0256\n", "Cut value 53\n", "Simulated annealing cut value 53\n", "Approximation ratio (vs SA) 1.0000\n" ] } ], "source": [ "# ── Classical baseline via simulated annealing ────────────────────\n", "# Run SA before any hardware calls to get a strong classical reference cut\n", "# value. SA is fast (seconds), needs no solver license, and reliably finds\n", "# near-optimal solutions on 40-node graphs. We use sa_cut as the denominator\n", "# for the approximation ratio instead of the looser QP upper bound.\n", "#\n", "# At each step we flip a random node and accept the move if it improves the\n", "# cut, or with probability exp(delta/T) otherwise. Temperature T decays\n", "# geometrically, allowing uphill moves early on to escape local minima.\n", "def simulated_annealing_maxcut(\n", " G, seed=0, T0=2.0, T_min=1e-4, alpha=0.995, n_steps=100_000\n", "):\n", " rng_sa = np.random.default_rng(seed)\n", " n = G.number_of_nodes()\n", " x = rng_sa.integers(0, 2, n)\n", " best_x = x.copy()\n", " best_cut = sum(1 for u, v in G.edges() if x[u] != x[v])\n", " T = T0\n", " for _ in range(n_steps):\n", " i = rng_sa.integers(0, n)\n", " delta = sum((-1 if x[i] != x[nb] else 1) for nb in G.neighbors(i))\n", " if delta > 0 or rng_sa.random() < np.exp(delta / T):\n", " x[i] ^= 1\n", " cut = sum(1 for u, v in G.edges() if x[u] != x[v])\n", " if cut > best_cut:\n", " best_cut, best_x = cut, x.copy()\n", " T = max(T * alpha, T_min)\n", " return best_x, best_cut\n", "\n", "\n", "sa_solution, sa_cut = simulated_annealing_maxcut(G_large)\n", "print(f\"Simulated annealing cut value: {sa_cut} (classical reference)\")\n", "\n", "# ── Step 3: Execution on hardware ───────────────────────────\n", "# A Session reserves the backend so the COBYLA iterations and final sampling\n", "# run back-to-back without re-queuing between jobs — important when the\n", "# optimizer submits many short jobs sequentially. All jobs are tagged with\n", "# \"TUT_WSQAOA\" for traceability in the IBM Quantum dashboard.\n", "#\n", "# EstimatorV2 with resilience_level=1 enables twirled readout error extinction\n", "# (TREX), which corrects systematic measurement bit-flip errors without extra\n", "# circuit overhead. 4096 shots per call balances estimation noise vs. job time.\n", "estimator_options = EstimatorOptions()\n", "estimator_options.resilience_level = 1\n", "estimator_options.default_shots = 4096\n", "estimator_options.environment.job_tags = [\"TUT_WSQAOA\"]\n", "\n", "# Align the cost observable with the physical qubit layout chosen by the transpiler\n", "cost_op_isa = cost_op_large.apply_layout(ws_isa_large.layout)\n", "ws_param_order_isa = list(ws_isa_large.parameters)\n", "\n", "ws_history_hw = []\n", "\n", "with Session(backend=backend) as session:\n", " estimator_hw = Estimator(mode=session, options=estimator_options)\n", "\n", " def hw_cost_fn(params):\n", " bound = ws_isa_large.assign_parameters(\n", " dict(zip(ws_param_order_isa, params))\n", " )\n", " energy = (\n", " estimator_hw.run([(bound, cost_op_isa)]).result()[0].data.evs.real\n", " )\n", " ws_history_hw.append(float(energy))\n", " print(\n", " f\" iter {len(ws_history_hw):>3d} = {energy:.4f}\", end=\"\\r\"\n", " )\n", " return float(energy)\n", "\n", " # Warm-start initialization: gamma=0 means the cost unitary is the identity on\n", " # the first call, so COBYLA immediately evaluates the warm-start state itself —\n", " # a much better starting signal than a random point.\n", " ws_params0_hw = np.concatenate(\n", " [np.zeros(p_large), np.full(p_large, np.pi / 4)]\n", " )\n", "\n", " ws_result_hw = minimize(\n", " hw_cost_fn,\n", " ws_params0_hw,\n", " method=\"COBYLA\",\n", " options={\"maxiter\": 150, \"rhobeg\": 0.3},\n", " )\n", " print(\n", " f\"\\nOptimization complete: energy={ws_result_hw.fun:.4f}, \"\n", " f\"iterations={len(ws_history_hw)}\"\n", " )\n", "\n", " # ── Step 3b: Sample the optimized circuit ──────────────────────────────────\n", " # Use 8192 shots for the final sample to get a reliable mode estimate.\n", " sampler_hw = Sampler(\n", " mode=session,\n", " options={\"environment\": {\"job_tags\": [\"TUT_WSQAOA\"]}},\n", " )\n", " ws_bound_hw = ws_isa_large.assign_parameters(\n", " dict(zip(ws_param_order_isa, ws_result_hw.x))\n", " )\n", " counts_hw = (\n", " sampler_hw.run([ws_bound_hw], shots=8192)\n", " .result()[0]\n", " .data.meas.get_counts()\n", " )\n", "\n", "best_bs_hw = max(counts_hw, key=counts_hw.get)\n", "best_count = counts_hw[best_bs_hw]\n", "total_shots = sum(counts_hw.values())\n", "\n", "# Decode: Qiskit returns bitstrings with qubit 0 at the rightmost position,\n", "# so reversing the string maps character index i to variable x_i.\n", "cut_val_hw, s0_hw, s1_hw = evaluate_cut(best_bs_hw[::-1], G_large)\n", "\n", "# Compare against simulated annealing.\n", "# A ratio >= 1.0 means WS-QAOA matched or beat the classical SA solution.\n", "# A ratio close to 1.0 (e.g. > 0.95) shows the quantum result is competitive.\n", "approx_ratio_hw = cut_val_hw / sa_cut\n", "print(\n", " f\"Most-probable bitstring frequency: {best_count}/{total_shots} \"\n", " f\"({100*best_count/total_shots:.1f}%)\"\n", ")\n", "print(\n", " f\"WS-QAOA cut: {cut_val_hw} | SA cut: {sa_cut} \"\n", " f\"| Approximation ratio vs SA: {approx_ratio_hw:.4f}\"\n", ")\n", "\n", "# Visualize both solutions side-by-side on the graph.\n", "# Blue = partition S, orange = partition S-bar.\n", "# Edges crossing between colors are the ones counted in the cut.\n", "fig, axes = plt.subplots(1, 2, figsize=(14, 6))\n", "for ax, assignment, cut, title in [\n", " (\n", " axes[0],\n", " list(sa_solution),\n", " sa_cut,\n", " f\"Simulated Annealing (cut={sa_cut})\",\n", " ),\n", " (\n", " axes[1],\n", " [int(b) for b in best_bs_hw[::-1]],\n", " cut_val_hw,\n", " f\"WS-QAOA hardware (cut={cut_val_hw})\",\n", " ),\n", "]:\n", " colors = [\n", " \"skyblue\" if assignment[i] == 0 else \"salmon\" for i in G_large.nodes()\n", " ]\n", " nx.draw(\n", " G_large,\n", " pos_large,\n", " with_labels=True,\n", " node_color=colors,\n", " node_size=400,\n", " font_size=7,\n", " ax=ax,\n", " )\n", " ax.set_title(title)\n", "plt.suptitle(\"Max-Cut partitions: SA vs WS-QAOA\", fontsize=13)\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "# ── Step 4: Convergence plot and summary ──────────────────────────────────\n", "# On real hardware the trace will be noisy (shot noise + gate errors), but the\n", "# overall downward trend confirms that COBYLA is making progress despite noise.\n", "fig, ax = plt.subplots(figsize=(7, 4))\n", "ax.plot(ws_history_hw, color=\"tab:orange\", label=\"WS-QAOA (hardware)\")\n", "ax.axhline(\n", " ws_result_hw.fun,\n", " color=\"tab:orange\",\n", " linestyle=\":\",\n", " label=f\"Final energy ({ws_result_hw.fun:.3f})\",\n", ")\n", "ax.set_xlabel(\"Optimizer call\")\n", "ax.set_ylabel(r\"$\\langle H_C \\rangle$\")\n", "ax.set_title(f\"WS-QAOA convergence on {backend.name} (40 qubits, p=1)\")\n", "ax.legend()\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "\n", "print(\"\\n=== Large Scale Summary ===\")\n", "print(f\"{'Metric':<38} {'Value':>10}\")\n", "print(\"-\" * 50)\n", "print(f\"{'Nodes / Edges':<38} {N_LARGE:>5} / {len(edges_large):<4}\")\n", "print(f\"{'QAOA layers (p)':<38} {p_large:>10}\")\n", "print(f\"{'Transpiled ECR gate count':<38} {ecr_count:>10}\")\n", "print(f\"{'Transpiled circuit depth':<38} {ws_isa_large.depth():>10}\")\n", "print(f\"{'Optimizer iterations':<38} {len(ws_history_hw):>10}\")\n", "print(f\"{'WS-QAOA energy (hardware)':<38} {ws_result_hw.fun:>10.4f}\")\n", "print(f\"{'Cut value':<38} {cut_val_hw:>10}\")\n", "print(f\"{'Simulated annealing cut value':<38} {sa_cut:>10}\")\n", "print(f\"{'Approximation ratio (vs SA)':<38} {approx_ratio_hw:>10.4f}\")" ] }, { "cell_type": "markdown", "id": "de87f93a", "metadata": {}, "source": [ "\n", "\n", "## 次のステップ\n", "\n", "\n", " この作品に興味を持たれた方は、以下の資料もご参考になるかもしれません：\n", "\n", " * **QAOAの上位層** ：値を増やして `p` 、回路層が増えるにつれて両アルゴリズムがどのように性能を向上させるか、また、層数が少ない場合でもWS-QAOAの優位性が維持されるかどうかを確認してください。\n", " * **Qiskit アドオン「optimization mapper** 」： [ドキュメント](https://qiskit.github.io/qiskit-addon-opt-mapper/)を参照し、さまざまな組み合わせ問題をモデル化したり、連続緩和問題に対して異なるソルバーを試したりしてみてください。\n", "\n", "\n" ] }, { "cell_type": "markdown", "id": "aafc36e2", "metadata": {}, "source": [ "\n", "\n", "## 参照\n", "\n", "\n", "\n", "[\\[1\\]](#Reference1) D. J. エガー、J. マレチェク、および S. Woerner, 「ウォームスタートによる量子最適化」, *『Quantum』*, 第5巻, p. 479, 2021年. [arXiv:2009.10095](https://arxiv.org/abs/2009.10095)\n", "\n", "\n", "\n", "[\\[2\\]](#Reference2) E. ファーヒ、J. ゴールドストーン、および S. Gutmann, 「量子近似最適化アルゴリズム」、『 [arXiv:1411.4028](https://arxiv.org/abs/1411.4028) 』、2014年。\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "id": "a1b8767d", "source": "© IBM Corp., 2017-2026" } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3" }, "hours": 1, "qpuSeconds": 540 }, "nbformat": 4, "nbformat_minor": 5 }