From 332936fb030457921d88b5e89f7451409be3e03f Mon Sep 17 00:00:00 2001 From: Albert Frisch Date: Wed, 13 Mar 2019 16:49:28 +0100 Subject: [PATCH] Add Aqua tutorial solving linear systems of equation with HHL --- .../general/linear_systems_of_equations.ipynb | 559 ++++++++++++++++++ 1 file changed, 559 insertions(+) create mode 100644 qiskit/aqua/general/linear_systems_of_equations.ipynb diff --git a/qiskit/aqua/general/linear_systems_of_equations.ipynb b/qiskit/aqua/general/linear_systems_of_equations.ipynb new file mode 100644 index 000000000..b9ad0233d --- /dev/null +++ b/qiskit/aqua/general/linear_systems_of_equations.ipynb @@ -0,0 +1,559 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# _*Qiskit Aqua: Solving linear systems of equations with the HHL algorithm*_\n", + "\n", + "***\n", + "### Contributors\n", + "David Bucher[1], Jan Mueggenburg[1], Gawel Kus[1], Isabel Haide[1], Shubha Deutschle[1], Harry Barowski[1], Dominik Steenken[1], and Albert Frisch[1]\n", + "### Affiliation\n", + "- [1]IBMQ" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The HHL algorithm (after the author’s surnames Harrow-Hassidim-Lloyd) [1] is a quantum algorithm to solve systems of linear equations $A \\vec{x} = \\vec{b}$. To perform this calculation quantum mechanically, we need in general 4 main steps requiring three qubit registers:\n", + "
    \n", + "
  1. First, we have to express the vector $\\vec{b}$ as a quantum state $|b\\rangle$ on a quantum register.
    2. \n", + "
  2. Now, we have to decompose $\\vec{b}$ into a superposition of eigenvectors of A remembering on the linear combination of the vector $\\vec{b}$. We achieve this using the Quantum Phase Estimation algorithm (Quantum Phase Estimation (QPE)). Since the matrix is hereby diagonalized wherefore $A$ is easily invertible.
    3. \n", + "
  3. The inversion of the eigenvector base of $A$ is achieved by rotating an ancillary qubit by an angle $\\arcsin \\left( \\frac{C}{\\lambda _{\\text{i}}} \\right)$ around the y-axis where $\\lambda_{\\text{i}}$ are the eigenvalues of $A$. Now, we obtain the state $A^{-1}|b\\rangle = |x \\rangle$.
    4. \n", + "
  4. We need to uncompute the register storing the eigenvalues using the inverse QPE. We measure the ancillary qubit whereby the measurement of 1 indicates that the matrix inversion was successful. The inverse QPE leaves the system in a state proportional to the solution vector $|x\\rangle$. In many cases one is not interested in the single vector elements of $|x\\rangle$ but only on certain properties. These are accessible by applying a problem-specific operator $M$ to the state $|x\\rangle$. Another use-case of the HHL algorithm is the implementation in a larger quantum program.
    5. \n", + "
\n", + "Currently only hermitian matrices with a dimension of $2^n$ are supported.\n", + "\n", + "Take into account that in the general case, the entries of $\\vec{x}$ can not be efficiently read out because we would need to know all coefficients describing the quantum state.\n", + "In the following examples, we ignore this constraint and show for our small linear system as a proof of principle that $\\vec{x}$ is calculated correctly.\n", + "\n", + "References:\n", + "- A. W. Harrow, A. Hassidim, and S. Lloyd, Phys. Rev. Lett. 103, 150502 (2009), e-print arXiv 0811.3171\n", + "- S. Barz, I. Kassal, M. Ringbauer, Y. Ole Lipp, B. Dakić, A. Aspuru-Guzik, and P. Walther, Sci Rep. 4: 6115 (2014), e-print arXiv 1302.1210" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from qiskit.aqua import run_algorithm\n", + "from qiskit.aqua.input import LinearSystemInput\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "params = {\n", + " 'problem': {\n", + " 'name': 'linear_system'\n", + " },\n", + " 'algorithm': {\n", + " 'name': 'HHL'\n", + " },\n", + " 'eigs': {\n", + " 'expansion_mode': 'suzuki',\n", + " 'expansion_order': 2,\n", + " 'name': 'EigsQPE',\n", + " 'num_ancillae': 3,\n", + " 'num_time_slices': 50\n", + " },\n", + " 'reciprocal': {\n", + " 'name': 'Lookup'\n", + " },\n", + " 'backend': {\n", + " 'name': 'statevector_simulator'\n", + " }\n", + "}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2x2 diagonal matrix" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, we show an example for running the HHL algorithm with Qiskit Aqua on a diagonal matrix as input\n", + "$$\n", + "A=\n", + "\\begin{bmatrix}\n", + "1 & 0 \\\\\n", + "0 & 2\n", + "\\end{bmatrix}$$ with the vector $$\\vec{b}= \\left( \\begin{array}{c}1 \\\\ 4 \\end{array} \\right)$$\n", + "The `result` dictionary contains several return values. The HHL solution for $\\vec{x}$ is accessible by the key `'solution_hhl'`. For comparison, also the classical solution of the linear system of equations is calculated using standard linear algebra functions in numpy. The fidelity between the HHL solution and the classical solution is also given in the output. Furthermore, the probability is shown with which HHL was running successfully, i.e. the HHL ancillary qubit has been measured to be $|1\\rangle$." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "matrix = [[1, 0], [0, 2]]\n", + "vector = [1, 4]\n", + "params['input'] = {\n", + " 'name': 'LinearSystemInput',\n", + " 'matrix': matrix,\n", + " 'vector': vector\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "hhl solution [1.05859-0.j 1.99245-0.j]\n", + "classical solution [1. 2.]\n", + "fidelity 0.999389\n", + "probability 0.024630\n" + ] + } + ], + "source": [ + "result = run_algorithm(params)\n", + "\n", + "print(\"hhl solution \", np.round(result['solution_hhl'], 5))\n", + "print(\"classical solution \", np.round(result['solution_classical'], 5))\n", + "print(\"fidelity %f\" % result['fidelity_hhl_to_classical'])\n", + "print(\"probability %f\" % result['probability_result'])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The probabilty that HHL runs successfully depends on the constant $C$ (see step 3. in the introduction). In the HHL algorithm, $C$ can be given as the parameter `scale` $\\in [0,1]$. In the above example `scale` is not defined in the `params` dictionary and the HHL algorithm initializes it to the smallest possible eigenvalue before execution. Alternatively, we can set `scale` to 0.5 and see how the results are influenced thereby." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "hhl solution [0.84664-0.j 2.01762-0.j]\n", + "classical solution [1. 2.]\n", + "fidelity 0.995605\n", + "probability 0.361437\n" + ] + } + ], + "source": [ + "params2 = params\n", + "params2['reciprocal'] = { \n", + " 'scale': 0.5\n", + "}\n", + "\n", + "result = run_algorithm(params2)\n", + "\n", + "print(\"hhl solution \", np.round(result['solution_hhl'], 5))\n", + "print(\"classical solution \", np.round(result['solution_classical'], 5))\n", + "print(\"fidelity %f\" % result['fidelity_hhl_to_classical'])\n", + "print(\"probability %f\" % result['probability_result'])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If you want to know how many qubits are required (circuit width) or how large the maximum number of gates applied to a single qubit (circuit depth) is, you can print it out by" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "circuit_depth 12255\n", + "circuit_width 7\n" + ] + } + ], + "source": [ + "print(\"circuit_depth\", result['circuit_depth'])\n", + "print(\"circuit_width\", result['circuit_width'])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2x2 non-diagonal matrix" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here we show an example for running the HHL algorithm with Qiskit Aqua on a non-diagonal matrix as input\n", + "$$\n", + "A=\n", + "\\begin{bmatrix}\n", + "1 & 3 \\\\\n", + "3 & 2\n", + "\\end{bmatrix}$$ with the vector $$\\vec{b}= \\left( \\begin{array}{c}1 \\\\ 1 \\end{array} \\right)$$\n" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "matrix = [[1, 3], [3, 2]]\n", + "vector = [1, 1]\n", + "params['input'] = {\n", + " 'name': 'LinearSystemInput',\n", + " 'matrix': matrix,\n", + " 'vector': vector\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "hhl solution [0.22147+0.j 0.22034-0.j]\n", + "classical solution [0.14286 0.28571]\n", + "fidelity 0.898454\n", + "probability 0.424639\n" + ] + } + ], + "source": [ + "result = run_algorithm(params)\n", + "\n", + "print(\"hhl solution \", np.round(result['solution_hhl'], 5))\n", + "print(\"classical solution \", np.round(result['solution_classical'], 5))\n", + "print(\"fidelity %f\" % result['fidelity_hhl_to_classical'])\n", + "print(\"probability %f\" % result['probability_result'])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Compared to the the first example, the circuit depth is increased approximately by a factor 2,5" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "circuit_depth 30253\n", + "circuit_width 7\n" + ] + } + ], + "source": [ + "print(\"circuit_depth\", result['circuit_depth'])\n", + "print(\"circuit_width\", result['circuit_width'])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 8x8 non-diagonal matrix" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For simplicity, we show a HHL execution of a linear systom of equations defined by the following 8x8 dimensional matrix\n", + "$$\n", + "A=\n", + "\\begin{bmatrix}\n", + "4 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\\\\n", + "0 & 3 & 0 & 0 & 0 & 0 & 0 & 0 \\\\\n", + "0 & 0 & 8 & 0 & 0 & 0 & 0 & 0 \\\\\n", + "0 & 0 & 0 & 5 & 0 & 0 & 0 & 0 \\\\\n", + "0 & 0 & 0 & 0 & 2 & 1 & 0 & 0 \\\\\n", + "0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 \\\\\n", + "0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\\\\n", + "1 & 0 & 0 & 0 & 0 & 0 & 0 & 5\n", + "\\end{bmatrix}$$ and the vector $$\\vec{b}= \\left( \\begin{array}{c}1 \\\\ 0 \\\\ 0 \\\\ 0 \\\\ 0 \\\\ 0 \\\\ 0 \\\\ 1 \\end{array} \\right)$$" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "matrix = [[4, 0, 0, 0, 0, 0, 0, 1],\n", + " [0, 3, 0, 0, 0, 0, 0, 0],\n", + " [0, 0, 8, 0, 0, 0, 0, 0],\n", + " [0, 0, 0, 5, 0, 0, 0, 0],\n", + " [0, 0, 0, 0, 2, 1, 0, 0],\n", + " [0, 0, 0, 0, 1, 1, 0, 0],\n", + " [0, 0, 0, 0, 0, 0, 1, 0],\n", + " [1, 0, 0, 0, 0, 0, 0, 5]]\n", + "vector = [1, 0, 0, 0, 0, 0, 0, 1]\n", + "params['input'] = {\n", + " 'name': 'LinearSystemInput',\n", + " 'matrix': matrix,\n", + " 'vector': vector\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "hhl solution [ 0.18195-0.j 0. -0.j -0. -0.j -0. +0.j 0. +0.j\n", + " 0. +0.j -0. -0.j 0.18041+0.j]\n", + "classical solution [0.21053 0. 0. 0. 0. 0. 0. 0.15789]\n", + "fidelity 0.981173\n", + "probability 0.935566\n" + ] + } + ], + "source": [ + "result = run_algorithm(params)\n", + "\n", + "print(\"hhl solution \", np.round(result['solution_hhl'], 5))\n", + "print(\"classical solution \", np.round(result['solution_classical'], 5))\n", + "print(\"fidelity %f\" % result['fidelity_hhl_to_classical'])\n", + "print(\"probability %f\" % result['probability_result'])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Considering the circuit depth, it is increased approximately by a factor 10 compared to the two dimensional matrices. The circuit width is increased by two additional qubits" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "circuit_depth 315268\n", + "circuit_width 9\n" + ] + } + ], + "source": [ + "print(\"circuit_depth\", result['circuit_depth'])\n", + "print(\"circuit_width\", result['circuit_width'])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 4x4 randomly-generated matrix" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, we show the application of HHL on a randomly-generated 4x4 matrix. We use the function `random_hermitian` to generate a random hermitian matrix and initialize the random seed to achieve reproducibility of the HHL run. Since the matrix can have negative eigenvalues, the `params` dictionary has to be modified by `\"negative_evals\": True` in `\"eigs\"` and `\"reciprocal\"`, respectively. We choose $$\\vec{b}= \\left( \\begin{array}{c}1 \\\\ 2 \\\\ 3 \\\\ 1 \\end{array} \\right)$$" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "from qiskit import Aer\n", + "from qiskit.aqua import QuantumInstance\n", + "from qiskit.aqua.algorithms.single_sample import HHL\n", + "from qiskit.aqua.utils import random_hermitian" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "It is needed for this example to define the \"initial_state\", the \"qft\" and the \"iqft\" additionally:" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "params3 = params\n", + "params3[\"reciprocal\"] = {\n", + " \"name\": \"Lookup\",\n", + " \"negative_evals\": True\n", + "}\n", + "params3[\"eigs\"] = {\n", + " \"expansion_mode\": \"suzuki\",\n", + " \"expansion_order\": 2,\n", + " \"name\": \"EigsQPE\",\n", + " \"negative_evals\": True,\n", + " \"num_ancillae\": 6,\n", + " \"num_time_slices\": 70\n", + "}\n", + "params3[\"initial_state\"] = {\n", + " \"name\": \"CUSTOM\"\n", + "}\n", + "params3[\"iqft\"] = {\n", + " \"name\": \"STANDARD\"\n", + "}\n", + "params3[\"qft\"] = {\n", + " \"name\": \"STANDARD\"\n", + "}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this example, we create an instance of the `HHL` class and run the algorithm with an input that is created programatically. To get the same pseudo-random matrix for every run, we set the random seed by using `np.random.seed(1)`." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "random matrix:\n", + "[[ 0.284+0.j -0.257-0.051j -0.124+0.033j 0.038+0.023j]\n", + " [-0.257+0.051j 0.404+0.j 0.067-0.079j 0.054+0.055j]\n", + " [-0.124-0.033j 0.067+0.079j 0.282-0.j 0.043+0.004j]\n", + " [ 0.038-0.023j 0.054-0.055j 0.043-0.004j 0.206+0.j ]]\n", + "HHL results:\n", + "hhl solution [ 79.9768 +4.52073j 60.28272 +3.09211j 37.51853 -9.5858j\n", + " -35.02324+26.46894j]\n", + "classical solution [ 76.1399 +1.92451j 57.30622 +1.20141j 35.96381-10.07775j\n", + " -32.03837+25.90593j]\n", + "fidelity 0.999946\n", + "probability 0.256771\n" + ] + } + ], + "source": [ + "# set the random seed to get the same pseudo-random matrix for every run\n", + "np.random.seed(1)\n", + "matrix = random_hermitian(4)\n", + "vector = [1, 2, 3, 1]\n", + "\n", + "print(\"random matrix:\")\n", + "m = np.array(matrix)\n", + "print(np.round(m, 3))\n", + "\n", + "algo_input = LinearSystemInput(matrix=matrix, vector=vector)\n", + "hhl = HHL.init_params(params3, algo_input)\n", + "backend = Aer.get_backend('statevector_simulator')\n", + "quantum_instance = QuantumInstance(backend=backend)\n", + "result_hhl = hhl.run(quantum_instance)\n", + "\n", + "print(\"HHL results:\")\n", + "print(\"hhl solution \", np.round(result_hhl['solution_hhl'], 5))\n", + "print(\"classical solution \", np.round(result_hhl['solution_classical'], 5))\n", + "print(\"fidelity %f\" % result_hhl['fidelity_hhl_to_classical'])\n", + "print(\"probability %f\" % result_hhl['probability_result'])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The circuit depth and width are" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "circuit_depth 973532\n", + "circuit_width 12\n" + ] + } + ], + "source": [ + "print(\"circuit_depth\", result_hhl['circuit_depth'])\n", + "print(\"circuit_width\", result_hhl['circuit_width'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}