{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"tags": [
"remove_cell"
]
},
"source": [
"# Solving linear systems of equations using HHL and its Qiskit implementation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this tutorial, we introduce the HHL algorithm, derive the circuit, and implement it using Qiskit. We show how to run the HHL on a simulator and on a five qubit device."
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"contents"
]
},
"source": [
"## Contents\n",
"1. [Introduction](#introduction)\n",
"2. [The HHL algorithm](#hhlalg)\n",
" 1. [Some mathematical background](#mathbackground)\n",
" 2. [Description of the HHL](#hhldescription)\n",
" 3. [Quantum Phase Estimation (QPE) within HHL](#qpe)\n",
" 4. [Non-exact QPE](#qpe2)\n",
"3. [Example 1: 4-qubit HHL](#example1)\n",
"4. [Qiskit Implementation](#implementation)\n",
" 1. [Running HHL on a simulator: general method](#implementationsim)\n",
" 2. [Running HHL on a real quantum device: optimised example](#implementationdev)\n",
"5. [Problems](#problems)\n",
"6. [References](#references)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1. Introduction \n",
"\n",
"Systems of linear equations arise naturally in many real-life applications in a wide range of areas, such as in the solution of Partial Differential Equations, the calibration of financial models, fluid simulation or numerical field calculation. The problem can be defined as, given a matrix $A\\in\\mathbb{C}^{N\\times N}$ and a vector $\\vec{b}\\in\\mathbb{C}^{N}$, find $\\vec{x}\\in\\mathbb{C}^{N}$ satisfying $A\\vec{x}=\\vec{b}$\n",
"\n",
"For example, take $N=2$, \n",
"\n",
"$$A = \\begin{pmatrix}1 & -1/3\\\\-1/3 & 1 \\end{pmatrix},\\quad \\vec{x}=\\begin{pmatrix} x_{1}\\\\ x_{2}\\end{pmatrix}\\quad \\text{and} \\quad \\vec{b}=\\begin{pmatrix}1 \\\\ 0\\end{pmatrix}$$\n",
"\n",
"Then the problem can also be written as find $x_{1}, x_{2}\\in\\mathbb{C}$ such that\n",
"$$\\begin{cases}x_{1} - \\frac{x_{2}}{3} = 1 \\\\ -\\frac{x_{1}}{3} + x_{2} = 0\\end{cases} $$\n",
"\n",
"A system of linear equations is called $s$-sparse if $A$ has at most $s$ non-zero entries per row or column. Solving an $s$-sparse system of size $N$ with a classical computer requires $\\mathcal{ O }(Ns\\kappa\\log(1/\\epsilon))$ running time using the conjugate gradient method ^{[1](#conjgrad)}. Here $\\kappa$ denotes the condition number of the system and $\\epsilon$ the accuracy of the approximation.\n",
"\n",
"The HHL is a quantum algorithm to estimate a function of the solution with running time complexity of $\\mathcal{ O }(\\log(N)s^{2}\\kappa^{2}/\\epsilon)$^{[2](#hhl)} when $A$ is a Hermitian matrix under the assumptions of efficient oracles for loading the data, Hamiltonian simulation and computing a function of the solution. This is an exponential speed up in the size of the system, however one crucial remark to keep in mind is that the classical algorithm returns the full solution, while the HHL can only approximate functions of the solution vector."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2. The HHL algorithm\n",
"\n",
"### A. Some mathematical background\n",
"The first step towards solving a system of linear equations with a quantum computer is to encode the problem in the quantum language. By rescaling the system, we can assume $\\vec{b}$ and $\\vec{x}$ to be normalised and map them to the respective quantum states $|b\\rangle$ and $|x\\rangle$. Usually the mapping used is such that $i^{th}$ component of $\\vec{b}$ (resp. $\\vec{x}$) corresponds to the amplitude of the $i^{th}$ basis state of the quantum state $|b\\rangle$ (resp. $|x\\rangle$). From now on, we will focus on the rescaled problem\n",
"\n",
"$$ A|x\\rangle=|b\\rangle$$\n",
"\n",
"Since $A$ is Hermitian, it has a spectral decomposition\n",
"$$\n",
"A=\\sum_{j=0}^{N-1}\\lambda_{j}|u_{j}\\rangle\\langle u_{j}|,\\quad \\lambda_{j}\\in\\mathbb{ R }\n",
"$$\n",
"where $|u_{j}\\rangle$ is the $j^{th}$ eigenvector of $A$ with respective eigenvalue $\\lambda_{j}$. Then,\n",
"$$\n",
"A^{-1}=\\sum_{j=0}^{N-1}\\lambda_{j}^{-1}|u_{j}\\rangle\\langle u_{j}|\n",
"$$\n",
"and the right hand side of the system can be written in the eigenbasis of $A$ as\n",
"$$\n",
"|b\\rangle=\\sum_{j=0}^{N-1}b_{j}|u_{j}\\rangle,\\quad b_{j}\\in\\mathbb{ C }\n",
"$$\n",
"It is useful to keep in mind that the goal of the HHL is to exit the algorithm with the readout register in the state\n",
"$$\n",
"|x\\rangle=A^{-1}|b\\rangle=\\sum_{j=0}^{N-1}\\lambda_{j}^{-1}b_{j}|u_{j}\\rangle\n",
"$$\n",
"Note that here we already have an implicit normalisation constant since we are talking about a quantum state."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### B. Description of the HHL algorithm \n",
"\n",
"The algorithm uses three quantum registers, all of them set to $|0\\rangle $ at the beginning of the algorithm. One register, which we will denote with the subindex $n_{l}$, is used to store a binary representation of the eigenvalues of $A$. A second register, denoted by $n_{b}$, contains the vector solution, and from now on $N=2^{n_{b}}$. There is an extra register, for the ancilla qubits. These are qubits used as intermediate steps in the individual computations but will be ignored in the following description since they are set to $|0\\rangle $ at the beginning of each computation and restored back to the $|0\\rangle $ state at the end of the individual operation.\n",
"\n",
"The following is an outline of the HHL algorithm with a high-level drawing of the corresponding circuit. For simplicity all computations are assumed to be exact in the ensuing description, and a more detailed explanation of the non-exact case is given in Section [2.D.](#qpe2).\n",
"\n",
"\n",
"\n",
"1. Load the data $|b\\rangle\\in\\mathbb{ C }^{N}$. That is, perform the transformation\n",
" $$ |0\\rangle _{n_{b}} \\mapsto |b\\rangle _{n_{b}} $$\n",
"2. Apply Quantum Phase Estimation (QPE) with\n",
"\t$$\n",
"\tU = e ^ { i A t } := \\sum _{j=0}^{N-1}e ^ { i \\lambda _ { j } t } |u_{j}\\rangle\\langle u_{j}|\n",
"\t$$\n",
"\tThe quantum state of the register expressed in the eigenbasis of $A$ is now\n",
"\t$$\n",
"\t\\sum_{j=0}^{N-1} b _ { j } |\\lambda _ {j }\\rangle_{n_{l}} |u_{j}\\rangle_{n_{b}}\n",
"\t$$\n",
" where $|\\lambda _ {j }\\rangle_{n_{l}}$ is the $n_{l}$-bit binary representation of $\\lambda _ {j }$.\n",
" \n",
"3. Add an ancilla qubit and apply a rotation conditioned on $|\\lambda_{ j }\\rangle$,\n",
"\t$$\n",
"\t\\sum_{j=0}^{N-1} b _ { j } |\\lambda _ { j }\\rangle_{n_{l}}|u_{j}\\rangle_{n_{b}} \\left( \\sqrt { 1 - \\frac { C^{2} } { \\lambda _ { j } ^ { 2 } } } |0\\rangle + \\frac { C } { \\lambda _ { j } } |1\\rangle \\right)\n",
"\t$$\n",
"\twhere $C$ is a normalisation constant.\n",
" \n",
"4. Apply QPE$^{\\dagger}$. Ignoring possible errors from QPE, this results in\n",
"\t$$\n",
"\t\\sum_{j=0}^{N-1} b _ { j } |0\\rangle_{n_{l}}|u_{j}\\rangle_{n_{b}} \\left( \\sqrt { 1 - \\frac {C^{2} } { \\lambda _ { j } ^ { 2 } } } |0\\rangle + \\frac { C } { \\lambda _ { j } } |1\\rangle \\right)\n",
"\t$$\n",
" \n",
"5. Measure the ancilla qubit in the computational basis. If the outcome is $1$, the register is in the post-measurement state\n",
"\t$$\n",
"\t\\left( \\sqrt { \\frac { 1 } { \\sum_{j=0}^{N-1} \\left| b _ { j } \\right| ^ { 2 } / \\left| \\lambda _ { j } \\right| ^ { 2 } } } \\right) \\sum _{j=0}^{N-1} \\frac{b _ { j }}{\\lambda _ { j }} |0\\rangle_{n_{l}}|u_{j}\\rangle_{n_{b}}\n",
"\t$$\n",
"\twhich up to a normalisation factor corresponds to the solution.\n",
"\n",
"6. Apply an observable $M$ to calculate $F(x):=\\langle x|M|x\\rangle$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### C. Quantum Phase Estimation (QPE) within HHL \n",
"\n",
"Quantum Phase Estimation is described in more detail in Chapter 3. However, since this quantum procedure is at the core of the HHL algorithm, we recall here the definition. Roughly speaking, it is a quantum algorithm which, given a unitary $U$ with eigenvector $|\\psi\\rangle_{m}$ and eigenvalue $e^{2\\pi i\\theta}$, finds $\\theta$. We can formally define this as follows.\n",
"\n",
"**Definition:** Let $U\\in\\mathbb{ C }^{2^{m}\\times 2^{m}}$ be unitary and let $|\\psi\\rangle_{m}\\in\\mathbb{ C }^{2^{m}}$ be one of its eigenvectors with respective eigenvalue $e^{2\\pi i\\theta}$. The **Quantum Phase Estimation** algorithm, abbreviated **QPE**, takes as inputs the unitary gate for $U$ and the state $|0\\rangle_{n}|\\psi\\rangle_{m}$ and returns the state $|\\tilde{\\theta}\\rangle_{n}|\\psi\\rangle_{m}$. Here $\\tilde{\\theta}$ denotes a binary approximation to $2^{n}\\theta$ and the $n$ subscript denotes it has been truncated to $n$ digits.\t\n",
"$$\n",
"\\operatorname { QPE } ( U , |0\\rangle_{n}|\\psi\\rangle_{m} ) = |\\tilde{\\theta}\\rangle_{n}|\\psi\\rangle_{m}\n",
"$$\n",
"\n",
"For the HHL we will use QPE with $U = e ^ { i A t }$, where $A$ is the matrix associated to the system we want to solve. In this case, \n",
"$$\n",
"e ^ { i A t } = \\sum_{j=0}^{N-1}e^{i\\lambda_{j}t}|u_{j}\\rangle\\langle u_{j}|\n",
"$$\n",
"Then, for the eigenvector $|u_{j}\\rangle_{n_{b}}$, which has eigenvalue $e ^ { i \\lambda _ { j } t }$, QPE will output $|\\tilde{\\lambda }_ { j }\\rangle_{n_{l}}|u_{j}\\rangle_{n_{b}}$. Where $\\tilde{\\lambda }_ { j }$ represents an $n_{l}$-bit binary approximation to $2^{n_l}\\frac{\\lambda_ { j }t}{2\\pi}$. Therefore, if each $\\lambda_{j}$ can be exactly represented with $n_{l}$ bits,\n",
"$$\n",
"\\operatorname { QPE } ( e ^ { i A t } , \\sum_{j=0}^{N-1}b_{j}|0\\rangle_{n_{l}}|u_{j}\\rangle_{n_{b}} ) = \\sum_{j=0}^{N-1}b_{j}|\\lambda_{j}\\rangle_{n_{l}}|u_{j}\\rangle_{n_{b}}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### D. Non-exact QPE \n",
"\n",
"In reality, the quantum state of the register after applying QPE to the initial state is\n",
"$$\n",
"\\sum _ { j=0 }^{N-1} b _ { j } \\left( \\sum _ { l = 0 } ^ { 2 ^ { n_{l} } - 1 } \\alpha _ { l | j } |l\\rangle_{n_{l}} \\right)|u_{j}\\rangle_{n_{b}}\n",
"$$\n",
"where\n",
"$$\n",
"\\alpha _ { l | j } = \\frac { 1 } { 2 ^ { n_{l} } } \\sum _ { k = 0 } ^ { 2^{n_{l}}- 1 } \\left( e ^ { 2 \\pi i \\left( \\frac { \\lambda _ { j } t } { 2 \\pi } - \\frac { l } { 2 ^ { n_{l} } } \\right) } \\right) ^ { k }\n",
"$$\n",
"\n",
"Denote by $\\tilde{\\lambda_{j}}$ the best $n_{l}$-bit approximation to $\\lambda_{j}$, $1\\leq j\\leq N$. Then we can relabel the $n_{l}$-register so that $\\alpha _ { l | j }$ denotes the amplitude of $|l + \\tilde { \\lambda } _ { j } \\rangle_{n_{l}}$. So now,\n",
"$$\n",
"\\alpha _ { l | j } : = \\frac { 1 } { 2 ^ { n_{l}} } \\sum _ { k = 0 } ^ { 2 ^ { n_{l} } - 1 } \\left( e ^ { 2 \\pi i \\left( \\frac { \\lambda _ { j } t } { 2 \\pi } - \\frac { l + \\tilde { \\lambda } _ { j } } { 2 ^ { n_{l} } } \\right) } \\right) ^ { k }\n",
"$$\n",
"If each $\\frac { \\lambda _ { j } t } { 2 \\pi }$ can be represented exactly with $n_{l}$ binary bits, then $\\frac { \\lambda _ { j } t } { 2 \\pi }=\\frac { \\tilde { \\lambda } _ { j } } { 2 ^ { n_{l} } }$ $\\forall j$. Therefore in this case $\\forall j$, $1\\leq j \\leq N$, it holds that $\\alpha _ { 0 | j } = 1$ and $\\alpha _ { l | j } = 0 \\quad \\forall l \\neq 0$. Only in this case we can write that the state of the register after QPE is \n",
"$$\n",
"\t\\sum_{j=0}^{N-1} b _ { j } |\\lambda _ {j }\\rangle_{n_{l}} |u_{j}\\rangle_{n_{b}}\n",
"$$\n",
"Otherwise, $|\\alpha _ { l | j }|$ is large if and only if $\\frac { \\lambda _ { j } t } { 2 \\pi } \\approx \\frac { l + \\tilde { \\lambda } _ { j } } { 2 ^ { n_{l} } }$ and the state of the register is\n",
"$$\n",
"\\sum _ { j=0 }^{N-1} \\sum _ { l = 0 } ^ { 2 ^ { n_{l} } - 1 } \\alpha _ { l | j } b _ { j }|l\\rangle_{n_{l}} |u_{j}\\rangle_{n_{b}}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3. Example: 4-qubit HHL\n",
"\n",
"Let's take the small example from the introduction to illustrate the algorithm. That is,\n",
"$$A = \\begin{pmatrix}1 & -1/3\\\\-1/3 & 1 \\end{pmatrix}\\quad \\text{and} \\quad |b\\rangle=\\begin{pmatrix}1 \\\\ 0\\end{pmatrix}$$\n",
"\n",
"We will use $n_{b}=1$ qubit to represent $|b\\rangle$, and later the solution $|x\\rangle$, $n_{l}=2$ qubits to store the binary representation of the eigenvalues and $1$ ancilla qubit to store whether the conditioned rotation, hence the algorithm, was successful.\n",
"\n",
"For the purpose of illustrating the algorithm, we will cheat a bit and calculate the eigenvalues of $A$ to be able to choose $t$ to obtain an exact binary representation of the rescaled eigenvalues in the $n_{l}$-register. However, keep in mind that for the HHL algorithm implementation one does not need previous knowledge of the eigenvalues. Having said that, a short calculation will give\n",
"$$\\lambda_{1} = 2/3\\quad\\text{and}\\quad\\lambda_{2}=4/3$$\n",
"\n",
"Recall from the previous section that the QPE will output an $n_{l}$-bit ($2$-bit in this case) binary approximation to $\\frac{\\lambda_ { j }t}{2\\pi}$. Therefore, if we set \n",
"$$t=2\\pi\\cdot \\frac{3}{8}$$\n",
"the QPE will give a $2$-bit binary approximation to\n",
"$$\\frac{\\lambda_ { 1 }t}{2\\pi} = 1/4\\quad\\text{and}\\quad\\frac{\\lambda_ { 2 }t}{2\\pi}=1/2$$\n",
"which is, respectively,\n",
"$$|01\\rangle_{n_{l}}\\quad\\text{and}\\quad|10\\rangle_{n_{l}}$$\n",
"\n",
"The eigenvectors are, respectively,\n",
"$$|u_{1}\\rangle=\\begin{pmatrix}1 \\\\ -1\\end{pmatrix}\\quad\\text{and}\\quad|u_{2}\\rangle=\\begin{pmatrix}1 \\\\ 1\\end{pmatrix}$$\n",
"Again, keep in mind that one does not need to compute the eigenvectors for the HHL implementation. In fact, a general Hermitian matrix $A$ of dimension $N$ can have up to $N$ different eigenvalues, therefore calculating them would take $\\mathcal{O}(N)$ time and the quantum advantage would be lost.\n",
"\n",
"We can then write $|b\\rangle$ in the eigenbasis of $A$ as\n",
"$$|b\\rangle _{n_{b}}=\\sum_{j=1}^{2}\\frac{1}{\\sqrt{2}}|u_{j}\\rangle _{n_{b}}$$\n",
"\n",
"Now we are ready to go through the different steps of the HHL algorithm. \n",
"\n",
"1. State preparation in this example is trivial since $|b\\rangle=|0\\rangle$.\n",
"2. Applying QPE will yield\n",
"$$\n",
"\\frac{1}{\\sqrt{2}}|01\\rangle|u_{1}\\rangle + \\frac{1}{\\sqrt{2}}|10\\rangle|u_{2}\\rangle\n",
"$$\n",
"3. Conditioned rotation with $C=3/8$ to compensate from having rescaled the eigenvalues gives\n",
"$$\\frac{1}{\\sqrt{2}}|01\\rangle|u_{1}\\rangle\\left( \\sqrt { 1 - \\frac { (3/8)^{2} } {(1/4)^{2} } } |0\\rangle + \\frac { 3/8 } { 1/4 } |1\\rangle \\right) + \\frac{1}{\\sqrt{2}}|10\\rangle|u_{2}\\rangle\\left( \\sqrt { 1 - \\frac { (3/8)^{2} } {(1/2)^{2} } } |0\\rangle + \\frac { 3/8 } { 1/2 } |1\\rangle \\right)\n",
"$$\n",
"$$\n",
"=\\frac{1}{\\sqrt{2}}|01\\rangle|u_{1}\\rangle\\left( \\sqrt { 1 - \\frac { 9 } {4 } } |0\\rangle + \\frac { 3 } { 2 } |1\\rangle \\right) + \\frac{1}{\\sqrt{2}}|10\\rangle|u_{2}\\rangle\\left( \\sqrt { 1 - \\frac { 9 } {16 } } |0\\rangle + \\frac { 3 } { 4 } |1\\rangle \\right)\n",
"$$\n",
"4. After applying QPE$^{\\dagger}$ the quantum computer is in the state\n",
"$$\n",
"\\frac{1}{\\sqrt{2}}|00\\rangle|u_{1}\\rangle\\left( \\sqrt { 1 - \\frac { 9 } {4 } } |0\\rangle + \\frac { 3 } { 2 } |1\\rangle \\right) + \\frac{1}{\\sqrt{2}}|00\\rangle|u_{2}\\rangle\\left( \\sqrt { 1 - \\frac { 9 } {16 } } |0\\rangle + \\frac { 3 } { 4 } |1\\rangle \\right)\n",
"$$\n",
"5. On outcome $1$ when measuring the ancilla qubit, the state is \n",
"$$\n",
"\\frac{\\frac{1}{\\sqrt{2}}|00\\rangle|u_{1}\\rangle\\frac { 3 } { 2 } |1\\rangle + \\frac{1}{\\sqrt{2}}|00\\rangle|u_{2}\\rangle\\frac { 3 } { 4 } |1\\rangle}{\\sqrt{45/32}}\n",
"$$\n",
"A quick calculation shows that\n",
"$$\n",
"\\frac{\\frac{3}{2\\sqrt{2}}|u_{1}\\rangle+ \\frac{3}{4\\sqrt{2}}|u_{2}\\rangle}{\\sqrt{45/32}} = \\frac{|x\\rangle}{||x||}\n",
"$$\n",
"6. Without using extra gates, we can compute the norm of $|x\\rangle$: it is the probability of measuring $1$ in the ancilla qubit from the previous step.\n",
"$$\n",
"P(|1\\rangle) = \\left(\\frac{3}{2\\sqrt{2}}\\right)^{2} + \\left(\\frac{3}{4\\sqrt{2}}\\right)^{2} = \\frac{45}{32} = ||x||^{2}\n",
"$$\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 4. Qiskit Implementation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that we have analytically solved the problem from the example we are going to use it to illustrate how to run the HHL on a quantum simulator and on the real hardware. For the quantum simulator, Qiskit Aqua already provides an implementation of the HHL algorithm requiring the matrix $A$ and $|b\\rangle$ as basic inputs. The main advantage is that it can take a general Hermitian matrix and an arbitrary initial state as inputs. This means that the algorithm is designed for a general purpose and does not optimise the circuit for a particular problem, which is problematic if the goal is to run the circuit on the existing real hardware. At the time of writing, the existing quantum computers are noisy and can only run small circuits. Therefore, in Section [4.B.](#implementationdev) we will see an optimised circuit that can be used for a class of problems to which our example belongs and mention the existing procedures to deal with noise in quantum computers."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## A. Running HHL on a simulator: general method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To run the HHL algorithm provided by Qiskit Aqua we just need to import the right modules and set the parameters as follows. In the worked out example we set the time of the Hamiltonian simulation to $t=2\\pi\\cdot \\frac{3}{8}$, however we will run the simulation without setting this parameter to show that knowledge of the eigenvalues is not required. Nonetheless, if the matrix has some structure it might be possible to obtain information about the eigenvalues and use it to choose a suitable $t$ and improve the accuracy of the solution returned by the HHL. As an exercise to see this, run the algorithm setting the time to $t=2\\pi\\cdot \\frac{3}{8}$. If done correctly, the fidelity of the solution should be $1$."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from qiskit import Aer, transpile, assemble\n",
"from qiskit.circuit.library import QFT\n",
"from qiskit.aqua import QuantumInstance, aqua_globals\n",
"from qiskit.quantum_info import state_fidelity\n",
"from qiskit.aqua.algorithms import HHL, NumPyLSsolver\n",
"from qiskit.aqua.components.eigs import EigsQPE\n",
"from qiskit.aqua.components.reciprocals import LookupRotation\n",
"from qiskit.aqua.operators import MatrixOperator\n",
"from qiskit.aqua.components.initial_states import Custom\n",
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def create_eigs(matrix, num_ancillae, num_time_slices, negative_evals):\n",
" ne_qfts = [None, None]\n",
" if negative_evals:\n",
" num_ancillae += 1\n",
" ne_qfts = [QFT(num_ancillae - 1), QFT(num_ancillae - 1).inverse()]\n",
"\n",
" return EigsQPE(MatrixOperator(matrix=matrix),\n",
" QFT(num_ancillae).inverse(),\n",
" num_time_slices=num_time_slices,\n",
" num_ancillae=num_ancillae,\n",
" expansion_mode='suzuki',\n",
" expansion_order=2,\n",
" evo_time=None, # This is t, can set to: np.pi*3/4\n",
" negative_evals=negative_evals,\n",
" ne_qfts=ne_qfts)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following function will be used to calculate the fidelity of solution returned by the HHL algorithm."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def fidelity(hhl, ref):\n",
" solution_hhl_normed = hhl / np.linalg.norm(hhl)\n",
" solution_ref_normed = ref / np.linalg.norm(ref)\n",
" fidelity = state_fidelity(solution_hhl_normed, solution_ref_normed)\n",
" print(\"Fidelity:\\t\\t %f\" % fidelity)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"matrix = [[1, -1/3], [-1/3, 1]]\n",
"vector = [1, 0]"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"orig_size = len(vector)\n",
"matrix, vector, truncate_powerdim, truncate_hermitian = HHL.matrix_resize(matrix, vector)\n",
"\n",
"# Initialize eigenvalue finding module\n",
"eigs = create_eigs(matrix, 3, 50, False)\n",
"num_q, num_a = eigs.get_register_sizes()\n",
"\n",
"# Initialize initial state module\n",
"init_state = Custom(num_q, state_vector=vector)\n",
"\n",
"# Initialize reciprocal rotation module\n",
"reci = LookupRotation(negative_evals=eigs._negative_evals, evo_time=eigs._evo_time)\n",
"\n",
"algo = HHL(matrix, vector, truncate_powerdim, truncate_hermitian, eigs,\n",
" init_state, reci, num_q, num_a, orig_size)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The reason to choose $t=2\\pi\\cdot \\frac{3}{8}$ was so that the rescaled eigenvalues could be represented exactly with $2$ binary digits. Since now this is not the case, the representation will be approximate, hence QPE not exact and the returned solution will be an approximation."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Solution:\t\t [1.13586-0.j 0.40896-0.j]\n",
"Classical Solution:\t [1.125 0.375]\n",
"Probability:\t\t 0.056291\n",
"Fidelity:\t\t 0.999432\n"
]
}
],
"source": [
"result = algo.run(QuantumInstance(Aer.get_backend('statevector_simulator')))\n",
"print(\"Solution:\\t\\t\", np.round(result['solution'], 5))\n",
"\n",
"result_ref = NumPyLSsolver(matrix, vector).run()\n",
"print(\"Classical Solution:\\t\", np.round(result_ref['solution'], 5))\n",
"\n",
"print(\"Probability:\\t\\t %f\" % result['probability_result'])\n",
"fidelity(result['solution'], result_ref['solution'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can print the resources used by the algorithm. The depth is the maximum number of gates applied to a single qubit, while the width is defined as the number of qubits required. We will also print the number of CNOTs since this number together with the width gives a good idea of whether running the circuit on current real hardware is feasible."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"circuit_width:\t 7\n",
"circuit_depth:\t 101\n",
"CNOT gates:\t 54\n"
]
}
],
"source": [
"print(\"circuit_width:\\t\", result['circuit_info']['width'])\n",
"print(\"circuit_depth:\\t\", result['circuit_info']['depth'])\n",
"print(\"CNOT gates:\\t\", result['circuit_info']['operations']['cx'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## B. Running HHL on a real quantum device: optimised example"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the previous section we ran the standard algorithm provided in Qiskit and saw that it uses $7$ qubits, has a depth of $102$ gates and requires a total of $54$ CNOT gates. These numbers are not feasible for the current available hardware, therefore we need to decrease these quantities. In particular, the goal will be to reduce the number of CNOTs by a factor of $5$ since they have worse fidelity than single-qubit gates. Furthermore, we can reduce the number of qubits to $4$ as was the original statement of the problem: the Qiskit method was written for a general problem and that is why it requires $3$ additional ancilla qubits.\n",
"\n",
"However, solely decreasing the number of gates and qubits will not give a good approximation to the solution on real hardware. This is because there are two sources of errors: those that occur during the run of the circuit and readout errors. \n",
"\n",
"Qiskit provides a module to mitigate the readout errors by individually preparing and measuring all basis states, a detailed treatment on the topic can be found in the paper by Dewes et al.^{[3](#readouterr)} To deal with the errors occurring during the run of the circuit, Richardson extrapolation can be used to calculate the error to the zero limit by running the circuit three times, each replacing each CNOT gate by $1$, $3$ and $5$ CNOTs respectively^{[4](#richardson)}. The idea is that theoretically the three circuits should produce the same result, but in real hardware adding CNOTs means amplifying the error. Since we know that we have obtained results with an amplified error, and we can estimate by how much the error was amplified in each case, we can recombine the quantities to obtain a new result that is a closer approximation to the analytic solution than any of the previous obtained values.\n",
"\n",
"Below we give the optimised circuit that can be used for any problem of the form\n",
"$$A = \\begin{pmatrix}a & b\\\\b & a \\end{pmatrix}\\quad \\text{and} \\quad |b\\rangle=\\begin{pmatrix}\\cos(\\theta) \\\\ \\sin(\\theta)\\end{pmatrix},\\quad a,b,\\theta\\in\\mathbb{R}$$\n",
"\n",
"The following optimisation was extracted from a work on the HHL for tridiagonal symmetric matrices^{[[5]](#tridi)}, this particular circuit was derived with the aid of the UniversalQCompiler software^{[[6]](#qcompiler)}.\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Depth: 26\n",
"CNOTS: 10\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
"