diff --git a/docs/examples/algorithms_and_protocols/pytket_qaoa_maxcut_example.ipynb b/docs/examples/algorithms_and_protocols/pytket_qaoa_maxcut_example.ipynb new file mode 100644 index 00000000..2b274a51 --- /dev/null +++ b/docs/examples/algorithms_and_protocols/pytket_qaoa_maxcut_example.ipynb @@ -0,0 +1,719 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "a3ba3449", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Solving Maxcut with QAOA\n", + "\n", + "**Download this notebook - {nb-download}`pytket_qaoa_maxcut_example.ipynb`**" + ] + }, + { + "cell_type": "markdown", + "id": "45668632", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "## The Max-Cut problem" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9456fbeb", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [], + "source": [ + "import networkx as nx\n", + "import matplotlib.pyplot as plt\n", + "\n", + "G = nx.Graph()\n", + "G.add_edges_from([(0, 1), (1, 2), (2, 0)])\n", + "plt.figure(figsize=(2, 2))\n", + "nx.draw(G, node_color=[\"red\", \"blue\", \"red\"])\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "cad41481", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "There are $2^3$ possible assignments of colour to nodes. In general there are $2^n$. The Max-cut problem can then be stated as that of finding the colour assignment which maximises the number of edges between vertices of a different colour." + ] + }, + { + "cell_type": "markdown", + "id": "7389bbe5", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "## Quantum Approximate Optimization Algorithm (QAOA)\n", + "\n", + "Introduced in 'A Quantum Approximate Optimization Algorithm' (found at https://arxiv.org/abs/1411.4028). The idea is to prepare a quantum state which encodes a solution to the Max-cut problem.\n", + "\n", + "\n", + "This is a variational algorithm, which is to say that a paramaterised state is prepared, with the parameters varied to improve the solution. We will have $2p$ parameters where p is our number of layers. In particular, the state prepared has the form \n", + "\n", + "\n", + "$$\n", + "| \\psi ( \\beta, \\gamma ) \\rangle = U ( \\beta_m ) U ( \\gamma_m ) ... U (\\beta_0) U ( \\gamma_0 ) | \\psi_0 \\rangle\n", + "$$\n", + "\n", + "where\n", + "\n", + "$$\n", + "U( \\beta_i ) = e^{i \\beta_i H_B} \\quad \\& \\quad U ( \\gamma_i) = e^{i \\gamma_i H_P}\n", + "$$\n", + "\n", + "with $H_P$ depending on the problem instance. " + ] + }, + { + "cell_type": "markdown", + "id": "d720b387", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "For the previous 3 vertex graph the *problem Hamiltonian* is\n", + "\n", + "$$\n", + "H_P = \\frac{3}{2} I - \\frac{1}{2} \\Big[ ( Z \\otimes Z \\otimes I ) + ( Z \\otimes I \\otimes Z ) + ( I \\otimes Z \\otimes Z ) \\Big]\n", + "$$\n", + "\n", + "where you will notice that there is a $ Z \\otimes Z$ acting between each vertex which is connected by an edge." + ] + }, + { + "cell_type": "markdown", + "id": "8ca2e1b9", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "The *mixer Hamiltonian* has the form \n", + "\n", + "$$\n", + "H_B = ( X \\otimes I \\otimes I ) + ( I \\otimes X \\otimes I ) + ( I \\otimes I \\otimes X )\n", + "$$\n", + "\n", + "where you will notice that there is an $X$ operator acting on each vertex." + ] + }, + { + "cell_type": "markdown", + "id": "de6b9e03", + "metadata": {}, + "source": [ + "## The cost function for Maxcut\n", + "\n", + "A solution to maxcut can be found by maximising the following cost function $C$ .\n", + "\n", + "\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "C= \\frac{1}{2}\\sum_{(i,j)} (1-z_i\\,z_j)\n", + "\\end{equation}\n", + "$$\n", + "\n", + "Here $z_i$ and $z_j$ are the the \"colours\" of each vertex. \n", + "\n", + "$$\n", + "\\begin{equation}\n", + "z_i,z_j \\in \\{0,1\\}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "61d4e798", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "The Pauli Z operator can be used to distinguish between the $|0\\rangle$ and $|1\\rangle$ basis states as these are eigenstates with eigenvalues $\\pm 1$ .\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "H_P = \\frac{1}{2}\\sum_{(i, \\,j)} (I-Z_i \\,Z_j)\n", + "\\end{equation}\n", + "$$\n", + "\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "H_B = \\sum_i X_i\n", + "\\end{equation}\n", + "$$\n", + "\n", + "Here we use the the convention that $X_i$ means a Pauli X operator will be applied to the \"ith\" qubit and the identity operator will be applied to all other qubits in the circuit." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "aeb6abd9", + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "import warnings\n", + "\n", + "warnings.filterwarnings(\"ignore\")" + ] + }, + { + "cell_type": "markdown", + "id": "3cbb783a", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Circuit Construction for QAOA" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "688f1332", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [], + "source": [ + "import networkx as nx\n", + "\n", + "max_cut_graph_edges = [(0, 1), (1, 2), (1, 3), (3, 4), (4, 5), (4, 6)]\n", + "n_nodes = 7\n", + "\n", + "max_cut_graph = nx.Graph()\n", + "max_cut_graph.add_edges_from(max_cut_graph_edges)\n", + "nx.draw(max_cut_graph, labels={node: node for node in max_cut_graph.nodes()})\n", + "\n", + "expected_results = [(0, 1, 0, 0, 1, 0, 0), (1, 0, 1, 1, 0, 1, 1)]" + ] + }, + { + "cell_type": "markdown", + "id": "18a5bd16", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Define Cost Hamiltonian: $\\gamma H$" + ] + }, + { + "cell_type": "markdown", + "id": "543f87ca", + "metadata": {}, + "source": [ + "$$\n", + "\\begin{equation}\n", + "H_P = \\frac{1}{2}\\sum_{(i, \\, j)} (I -Z_i \\,Z_j)\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "6da499ac", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "$$\n", + "\\begin{equation}\n", + "H_P = 3 I^{\\otimes 7} -0.5 \\big[ Z_0 Z_1 + Z_1 Z_2 +Z_1 Z_3 +Z_3 Z_4 +Z_4 Z_5 +Z_4 Z_6 \\big]\n", + "\\end{equation}\n", + "$$\n", + "\n", + "Using the same index convention as above" + ] + }, + { + "cell_type": "markdown", + "id": "785ff56c", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Hamiltonian Circuit" + ] + }, + { + "cell_type": "markdown", + "id": "4690b787", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Construction of the Mixer Hamiltonian: $\\beta B$" + ] + }, + { + "cell_type": "markdown", + "id": "4d128a70", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Define the Initial State" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0a9db628", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [], + "source": [ + "from pytket import Circuit\n", + "from pytket.circuit.display import render_circuit_jupyter as draw\n", + "\n", + "\n", + "def qaoa_initial_circuit(n_qubits: int) -> Circuit:\n", + " c = Circuit(n_qubits)\n", + " for i in range(n_qubits):\n", + " c.H(i)\n", + " return c\n", + "\n", + "\n", + "superposition_circuit = qaoa_initial_circuit(n_nodes)\n", + "\n", + "draw(superposition_circuit)" + ] + }, + { + "cell_type": "markdown", + "id": "da759b59", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Construct QAOA Circuit" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "57c9ed7e", + "metadata": {}, + "outputs": [], + "source": [ + "def build_cost_layer(graph: nx.Graph, gamma_val: float) -> Circuit:\n", + " circ = Circuit(graph.number_of_nodes())\n", + "\n", + " for i, j in graph.edges:\n", + " circ.ZZPhase(-gamma_val / 2, i, j)\n", + "\n", + " return circ" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "5f43f97a", + "metadata": {}, + "outputs": [], + "source": [ + "def build_mixer_layer(n_vertices: int, beta_val: float) -> Circuit:\n", + " circ = Circuit(n_vertices)\n", + "\n", + " for qubit in circ.qubits:\n", + " circ.Rx(beta_val, qubit)\n", + "\n", + " return circ" + ] + }, + { + "cell_type": "markdown", + "id": "359a1a0f-e92e-40ae-bbe6-ce960b118f49", + "metadata": {}, + "source": [ + "Now lets define a function to create our entire QAOA circuit. For $p$ QAOA layers we expect that our circuit will require $2p$ parameters. Here we will pass and cost mixer parameters in as a list where the length of the list defines the number of layers." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "23f8910a", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [], + "source": [ + "def build_qaoa_circuit(\n", + " graph: nx.Graph, mixer_angles: list[float], cost_angles: list[float]\n", + ") -> Circuit:\n", + "\n", + " n_nodes = graph.number_of_nodes()\n", + "\n", + " assert len(mixer_angles) == len(cost_angles)\n", + "\n", + " qaoa_circuit = qaoa_initial_circuit(n_nodes)\n", + "\n", + " # Iteratively append \"p\" cost and mixer layers\n", + " for cost_angle, mixer_angle in zip(cost_angles, mixer_angles):\n", + " qaoa_circuit.append(build_cost_layer(graph, cost_angle))\n", + " qaoa_circuit.add_barrier(list(range(n_nodes)))\n", + " qaoa_circuit.append(build_mixer_layer(n_nodes, mixer_angle))\n", + "\n", + " qaoa_circuit.measure_all()\n", + " return qaoa_circuit" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f97de321", + "metadata": {}, + "outputs": [], + "source": [ + "from networkx import path_graph\n", + "\n", + "three_vertex_path = path_graph(3)\n", + "\n", + "draw(build_qaoa_circuit(three_vertex_path, [0.8, 0.1], [0.75, 0.6]))" + ] + }, + { + "cell_type": "markdown", + "id": "bc2f8939-41b7-476b-a5a8-09de07211079", + "metadata": {}, + "source": [ + "We also need to extract our energy expectation values from a `BackendResult` object after our circuit is processed by the device/simulator. We do this with the `energy_from_result` function below. Note that the fact that the maxcut Hamiltonian contains only commuting terms means that we do not need to calculate our energy expectation using multiple measurement circuits. This may not the the case for a different problem Hamiltonian." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "df387eea-4198-428e-9b92-4f3bceb12f0e", + "metadata": {}, + "outputs": [], + "source": [ + "from pytket.backends.backendresult import BackendResult\n", + "\n", + "\n", + "def energy_from_result(graph: nx.Graph, results: BackendResult) -> float:\n", + " energy = 0.0\n", + " dist = results.get_distribution()\n", + " for i, j in graph.edges:\n", + " energy += sum((meas[i] ^ meas[j]) * prob for meas, prob in dist.items())\n", + "\n", + " return energy" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "e5abad7b-e989-4156-9708-3d8c97d8ca2a", + "metadata": {}, + "outputs": [], + "source": [ + "from pytket.backends.backend import Backend\n", + "from pytket.passes import RemoveBarriers\n", + "from typing import Callable\n", + "import numpy as np\n", + "\n", + "\n", + "def eval_qaoa_energy(\n", + " backend: Backend,\n", + " guess_mixer_angles: np.array,\n", + " guess_cost_angles: np.array,\n", + " graph: nx.Graph,\n", + " seed: int,\n", + " shots: int = 5000,\n", + " compiler_pass: Callable[[Circuit], bool] | None = None,\n", + ") -> tuple[float, BackendResult]:\n", + "\n", + " # Build Circuit\n", + " qaoa_circuit = build_qaoa_circuit(\n", + " graph=graph, cost_angles=guess_cost_angles, mixer_angles=guess_mixer_angles\n", + " )\n", + "\n", + " # Compile\n", + " RemoveBarriers().apply(qaoa_circuit)\n", + " if compiler_pass:\n", + " compiler_pass(qaoa_circuit)\n", + "\n", + " # Execute\n", + " result: BackendResult = backend.run_circuit(qaoa_circuit, shots, seed=seed)\n", + "\n", + " return energy_from_result(graph, result), result" + ] + }, + { + "cell_type": "markdown", + "id": "2c01c28b", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Optimise Energy by Guessing Parameters" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "0a44bed8", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [], + "source": [ + "def solve_maxcut_instance(\n", + " graph: nx.Graph,\n", + " backend: Backend,\n", + " iterations: int = 100,\n", + " p_value: int = 3,\n", + " n_shots: int = 5000,\n", + " seed: int = 12345,\n", + " compiler_pass: Callable[[Circuit], bool]| None = None,\n", + ") -> tuple[BackendResult, np.array, np.array]:\n", + "\n", + " highest_energy = 0\n", + " best_guess_mixer_angles = [0 for _ in range(p_value)]\n", + " best_guess_cost_angles = [0 for _ in range(p_value)]\n", + "\n", + " rng = np.random.default_rng(seed)\n", + "\n", + " for _ in range(iterations):\n", + " guess_mixer_angles = rng.uniform(0, 1, p_value)\n", + " guess_cost_angles = rng.uniform(0, 1, p_value)\n", + " qaoa_energy, result = eval_qaoa_energy(\n", + " backend,\n", + " guess_mixer_angles,\n", + " guess_cost_angles,\n", + " seed=seed,\n", + " shots=n_shots,\n", + " graph=graph,\n", + " compiler_pass=None,\n", + " )\n", + "\n", + " if qaoa_energy > highest_energy:\n", + " print(\"new highest energy found: \", qaoa_energy)\n", + "\n", + " best_guess_mixer_angles = np.round(guess_mixer_angles, 3)\n", + " best_guess_cost_angles = np.round(guess_cost_angles, 3)\n", + " highest_energy = qaoa_energy\n", + " best_result: BackendResult = result\n", + "\n", + " print(\"highest energy: \", highest_energy)\n", + " print(\"best guess mixer angles: \", best_guess_mixer_angles)\n", + " print(\"best guess cost angles: \", best_guess_cost_angles)\n", + " best_outputs = tuple([best_result, best_guess_cost_angles, best_guess_mixer_angles])\n", + " return best_outputs" + ] + }, + { + "cell_type": "markdown", + "id": "d22226cc", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Calculate the State for the final Parameters" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "e7afb38e", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [], + "source": [ + "from pytket.extensions.qiskit import AerBackend\n", + "\n", + "backend = AerBackend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aaea7e2f", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [], + "source": [ + "%%time\n", + "qaoa_result, cost_angles, mixer_angles = solve_maxcut_instance(\n", + " backend=backend,\n", + " compiler_pass=None,\n", + " n_shots=5000,\n", + " iterations=100,\n", + " graph=max_cut_graph,\n", + " seed=12345,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3b86301e-7645-4553-be38-3ebf89eedd37", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "\n", + "def plot_maxcut_results(result: BackendResult, n_strings: int) -> None:\n", + " \"\"\"\n", + " Plots Maxcut results in a barchart with the two most common bitstrings highlighted in green.\n", + " \"\"\"\n", + " counts_dict = result.get_counts()\n", + " sorted_shots = counts_dict.most_common()\n", + " n_shots = sum(counts_dict.values())\n", + "\n", + " n_most_common_strings = sorted_shots[:n_strings]\n", + " x_axis_values = [str(entry[0]) for entry in n_most_common_strings] # basis states\n", + " y_axis_values = [entry[1] for entry in n_most_common_strings] # counts\n", + " num_successful_shots = sum(y_axis_values[:2])\n", + " print(f\"Success ratio {num_successful_shots/n_shots} \")\n", + "\n", + " fig = plt.figure()\n", + " ax = fig.add_axes([0, 0, 1.5, 1])\n", + " color_list = [\"green\"] * 2 + ([\"orange\"] * (len(x_axis_values) - 2))\n", + " ax.bar(\n", + " x=x_axis_values,\n", + " height=y_axis_values,\n", + " color=color_list,\n", + " )\n", + " ax.set_title(label=\"Maxcut Results\")\n", + " plt.ylim([0, 0.25 * n_shots])\n", + " plt.xlabel(\"Basis State\")\n", + " plt.ylabel(\"Number of Shots\")\n", + " plt.show()\n", + "\n", + "\n", + "plot_maxcut_results(qaoa_result, 6)" + ] + }, + { + "cell_type": "markdown", + "id": "6e36c4fb-a118-4ab8-be01-77a674f273e3", + "metadata": {}, + "source": [ + "Here the binary strings in the results correspond to the two optimal colourings of our graph." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ffce2a97-902c-44d8-8d64-b25498752907", + "metadata": {}, + "outputs": [], + "source": [ + "G = nx.Graph()\n", + "G.add_edges_from(max_cut_graph_edges)\n", + "\n", + "H = nx.Graph()\n", + "H.add_edges_from(max_cut_graph_edges)\n", + "\n", + "plt.figure(1)\n", + "nx.draw(\n", + " G,\n", + " labels={node: node for node in max_cut_graph.nodes()},\n", + " node_color=[\"red\", \"blue\", \"red\", \"red\", \"blue\", \"red\", \"red\"],\n", + ")\n", + "plt.figure(2)\n", + "nx.draw(\n", + " H,\n", + " labels={node: node for node in max_cut_graph.nodes()},\n", + " node_color=[\"blue\", \"red\", \"blue\", \"blue\", \"red\", \"blue\", \"blue\"],\n", + ")\n", + "\n", + "plt.show()" + ] + } + ], + "metadata": { + "celltoolbar": "Slideshow", + "kernelspec": { + "display_name": ".venv", + "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.11.1" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}