VQE_submission.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Variational Quantum Eigensolver for 2-qubit system\n",
    "\n",
    "This is an attempt at solving task 4 of the screening tasks i.e. finding the lowest eigenvalue for the given matrix using VQE like circuits.\n",
    "\n",
    "## Introduction\n",
    "Variational Quantum Eigensolver is an algorithm that helps us find an upper bound on the lowest eigen value of a given Hamiltonian. This notebook will not go into much depth on the theoretical aspects of the algorithm such as the variational principle and what a Hamiltonian. This article<sup>[[1]]</sup> by Michał Stęchły and the original paper on VQE<sup>[[2]]</sup> do a great job of explaining them. Here, the focus will be on the implementation of the algorithm.  \n",
    "<br>\n",
    "\n",
    "<div align=\"center\">\n",
    "    <img alt=\"Variational Quantum Algorithm\" src=\"./VQE_algorithm.png\"><br>\n",
    "    Variational Quantum Algorithm<sup>[2]</sup>\n",
    "</div>\n",
    "\n",
    "<br><br>\n",
    "The steps involved in the algorithm can be summarized as follows:  \n",
    "\n",
    "1. Design a parameterized circuit (ansatz) for trial state preparation. Design the quantum modules based on the tensor products of Pauli matrices obtained by decomposing the Hamiltonian.  \n",
    "\n",
    "2. Initialize the parameters of the circuit and prepare a trial state.  \n",
    "\n",
    "3. Pass the trial state to the quantum modules to perform measurements.  \n",
    "\n",
    "4. Calculate the expectation value $\\left\\langle H_{1} \\right\\rangle$ to $\\left\\langle H_{n} \\right\\rangle$ based on the counts of measurements obtained from quantum modules.  \n",
    "\n",
    "5. Add all the expectation values using classical adder to obtain expectation value $\\left\\langle H \\right\\rangle$.  \n",
    "\n",
    "6. Obtain new input parameters based on the expectation value using a classical optimizer.  \n",
    "\n",
    "7. Repeat steps 2 to 7 until the expectation value is optimised. This expectation value is the upper bound on the lowest eigenvalue of the given Hamiltonian.\n",
    "\n",
    "\n",
    "[1]: (https://www.mustythoughts.com/variational-quantum-eigensolver-explained)\n",
    "[2]: (https://arxiv.org/abs/1304.3061)\n"
   ]
  },
  {
   "source": [
    "Before we proceed with the code, let's make sure that all the required dependencies are installed on your system. It is recommended that these dependencies be installed in a virtual environment so that the global ones are not messed up. \n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "!pip install -U -r requirements.txt\n",
    "\n",
    "from IPython.display import clear_output\n",
    "clear_output()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "''' Importing the required modules and defining global constants '''\n",
    "import numpy as np\n",
    "pi = np.pi\n",
    "from qiskit import *\n",
    "from scipy import optimize\n",
    "import pandas as pd\n",
    "from matplotlib.pyplot import savefig\n",
    "from qiskit.providers.ibmq import *"
   ]
  },
  {
   "source": [
    "## Decomposing given Hamiltonian into Pauli products\n",
    "\n",
    "The first step in the algorithm requires us to decompose the given Hamiltonian matrix into a linear combination of Pauli matrices and tensor products of Pauli matrices (henceforth called \"Pauli product\"). This is done as a quantum computer can efficiently evaluate the expectation value of Pauli products<sup>[[2]]</sup>.\n",
    "\n",
    "After digging around the internet,  these resources<sup>[[3]][[4]][[5]][[8]]</sup> were found that leverage the properties of Pauli matrices, Hermitian matrices and the Hilbert-Schmidt inner product to calculate the coefficients of the decomposition. If we represent the decomposed 4 x 4 Hamiltonian matrix H as\n",
    "\n",
    "$$H = \\ \\sum_{i,j = I,X,Y,Z}^{}{a_{\\text{ij}}\\left( \\sigma_{i} \\otimes \\sigma_{j} \\right)}$$\n",
    "\n",
    "then we can use the relation\n",
    "\n",
    "$$a_{\\text{ij}} = \\frac{1}{4}tr(\\left( \\sigma_{i} \\otimes \\sigma_{j} \\right)H)\\ where\\ i,\\ j = I,\\ X,\\ Y,\\ Z$$\n",
    "\n",
    "to calculate the coefficient of each Pauli product.\n",
    "\n",
    "[2]: https://arxiv.org/abs/1304.3061\n",
    "[3]: https://quantumcomputing.stackexchange.com/questions/11899/example-of-hamiltonian-decomposition-into-pauli-matrices\n",
    "[4]: https://quantumcomputing.stackexchange.com/questions/8725/can-arbitrary-matrices-be-decomposed-using-the-pauli-basis\n",
    "[5]: https://quantumcomputing.stackexchange.com/questions/6882/decomposition-of-a-matrix-in-the-pauli-basis\n",
    "[8]: https://michaelgoerz.net/notes/decomposing-two-qubit-hamiltonians-into-pauli-matrices.html\n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def decompose_matrix(matrix):\n",
    "    \"\"\" \n",
    "    This function uses the formula describe above to calculate the coefficient of each Pauli product. \n",
    "  \n",
    "    It essentially decomposes the 4 x 4 hamiltonian matrix into a linear combination of Pauli products.  \n",
    "  \n",
    "    Parameters: \n",
    "    matrix (np.array): A 4 x 4 hamiltonian matrix\n",
    "  \n",
    "    Returns: \n",
    "    dict: Dictionary of coefficients of each Pauli product\n",
    "    \"\"\"\n",
    "\n",
    "    pauli_I = np.array([[1,  0], \n",
    "                        [0,  1]], dtype=complex)\n",
    "    pauli_X = np.array([[0,  1], \n",
    "                        [1,  0]], dtype=complex)\n",
    "    pauli_Y = np.array([[0, -1j], \n",
    "                        [1j,  0]], dtype=complex)\n",
    "    pauli_Z = np.array([[1,  0], \n",
    "                        [0, -1]], dtype=complex) \n",
    "    pauli_matrices = [[\"I\", pauli_I], [\"X\", pauli_X], [\"Y\", pauli_Y], [\"Z\", pauli_Z]]\n",
    "\n",
    "    coefficient_dict = {}\n",
    "\n",
    "    for pauli_matrix_1 in pauli_matrices:\n",
    "        for pauli_matrix_2 in pauli_matrices:\n",
    "            tensor_product = np.kron(pauli_matrix_1[1], pauli_matrix_2[1])\n",
    "            coefficient_dict[f\"{pauli_matrix_1[0]}{pauli_matrix_2[0]}\"] = 0.25 * np.trace(np.matmul(tensor_product, given_matrix))\n",
    "\n",
    "    return coefficient_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "Coefficient of each tensor product in the decomposition of given matrix: \n\n {'II': (0.5+0j), 'IX': 0j, 'IY': 0j, 'IZ': 0j, 'XI': 0j, 'XX': (-0.5+0j), 'XY': 0j, 'XZ': 0j, 'YI': 0j, 'YX': 0j, 'YY': (-0.5+0j), 'YZ': 0j, 'ZI': 0j, 'ZX': 0j, 'ZY': 0j, 'ZZ': (0.5+0j)}\n"
    }
   ],
   "source": [
    "given_matrix = np.array([[1,  0,  0,  0], \n",
    "                         [0,  0, -1,  0],\n",
    "                         [0, -1,  0,  0],\n",
    "                         [0,  0,  0,  1]], dtype=np.float64)\n",
    "print(\"Coefficient of each tensor product in the decomposition of given matrix: \\n\\n\", decompose_matrix(given_matrix))"
   ]
  },
  {
   "source": [
    "It is seen from the above output that the given matrix is decomposed into a linear combination of $\\sigma_{I} \\otimes \\sigma_{I}$, $\\sigma_{X} \\otimes \\sigma_{X}$, $\\sigma_{Y} \\otimes \\sigma_{Y}$ and $\\sigma_{Z} \\otimes \\sigma_{Z}$ i.e.\n",
    "\n",
    "$$H = \\frac{1}{2}\\left( \\sigma_{I} \\otimes \\sigma_{I} \\right) - \\frac{1}{2}\\ \\left( \\sigma_{X} \\otimes \\sigma_{X} \\right) - \\frac{1}{2}{(\\sigma}_{Y} \\otimes \\sigma_{Y}) + \\frac{1}{2}\\left( \\sigma_{Z} \\otimes \\sigma_{Z} \\right)$$"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "## Creating Trial States using Parametrized Ansatz\n",
    "\n",
    "The next step is to design the circuit for trial state preparation. The goal is to create a state that is exactly the eigen state or very close to it. We could do this by iterating over all possible states in the Hilbert Space until we find the correct one but this would be computationally intensive. We need a better way to do so. \n",
    "\n",
    "An ansatz, in general sense, is an educated guess or an additional assumption made to help solve a problem, and which is later verified to be part of the solution by its results<sup>[[6]]</sup>. In this case, an ansatz is a set of parametrized gates that gives us access to a portion of the Hilbert space. Since it is parameterized, the parameters can be varied iteratively to access the set of states represented by it. Choosing a good ansatz is key as it should represent a sufficient part of the Hilbert space and be shallow (to be less computationally intensive) and should not have too many parameters (optimizing it would become difficult)<sup>[[1]]</sup>. \n",
    "\n",
    "After experimenting with different ansatze, it was found that the ansatz in the hint i.e. ${RX(\\theta)}_{1}{(CX)H}_{1}$ was the best ansatz for the given hamiltonian as it contains the eigenvector, corresponding to the lowest eigenvalue, as a state in its subspace. \n",
    "\n",
    "[6]: https://en.wikipedia.org/wiki/Ansatz\n",
    "[1]: https://www.mustythoughts.com/variational-quantum-eigensolver-explained\n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def create_trial_state_circuit(parameters):\n",
    "    \"\"\" \n",
    "    Creates a parameterized circuit (ansatz) that prepares the trial state based \n",
    "    on the parameters received as input.\n",
    "  \n",
    "    Parameters: \n",
    "    parameters (np.array): List of angles that act as parameters for the circuit\n",
    "  \n",
    "    Returns: \n",
    "    QuantumCircuit(2, 2): Quantum circuit that prepares the trial state\n",
    "    \"\"\"\n",
    "    trial_state_circuit = QuantumCircuit(2, 2)\n",
    "\n",
    "    trial_state_circuit.h(0)\n",
    "    trial_state_circuit.cx(0, 1)\n",
    "    trial_state_circuit.rx(parameters[0], 0)\n",
    "\n",
    "    trial_state_circuit.barrier()\n",
    "    # trial_state_circuit.draw(output='mpl').savefig('./circuits/AnsatzHint.png') # This statement was used during the experimentation phase to store the circuits\n",
    "    return trial_state_circuit"
   ]
  },
  {
   "source": [
    "## Creating quantum modules to perform measurements\n",
    "\n",
    "It is impossible to know the exact state of a qubit as any external interaction collapses the qubit into one of the basis states. Sp, to get an approximate idea of what the trial state would have been, the same circuit is repeatedly prepared and measurements are performed measurements to get the counts of each output state. These can be used to calculate the probability of each output state and which can in turn be used to calculate the expectation value.\n",
    "\n",
    "The problem that arises here is that separate circuits (i.e. quantum modules) are needed for each Pauli Product. The reason being that to calculate the expectation value for $\\sigma_{X}$ and $\\sigma_{y}$, measurements have to be performed and counts have to be obtained in the X and Y basis respectively. Since we cannot directly measure in an arbitrary basis, transformations have to be performed on the trial state to convert it into Z basis (i.e. the $\\left| 0 \\right\\rangle\\ and\\ |1\\rangle$ basis). For $\\sigma_{X}$, a Hadamard (or $H$) gate is applied and for $\\sigma_{y}$, $S^{\\dagger}$ followed by $H$ gate are applied for the transformation.\n",
    "\n",
    "$$S^{\\dagger} = \\ \\begin{bmatrix} 1 & 0 \\\\ 0 & - i \\\\ \\end{bmatrix}\\ \\ \\ \\ \\ \\ \\ \\ \\ \\ H = \\frac{1}{\\sqrt{2}}\\begin{bmatrix} 1 & 1 \\\\ 1 & - 1 \\\\ \\end{bmatrix}$$\n",
    "\n",
    "For $\\sigma_{Z}$, there is no need for a transformation as its expectation value requires counts in the Z basis. For $\\sigma_{I}$, there is no need to create a circuit as its expectation value is always 1. \n",
    "\n",
    "For a 2-qubit system, the expectation value for tensor products has to be calculated. This can be done by doing these transformations on individual qubits according to the Pauli matrix and then doing a measurement. E.g. For calculating the expectation value of the $\\sigma_{X} \\otimes \\sigma_{Y}$, a $H$ gate is applied on the first qubit and a $S^{\\dagger}$ followed by $H$ gate are applied on the second qubit and measurements are performed. "
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quantum_module_simulator(trial_state_circuit, pauli_product, number_of_shots):\n",
    "    \"\"\" \n",
    "    This is a generalized function that adds transformations and performs measurements on the trial\n",
    "    state based on the pauli matrices. The measurements are performed repeatedly to obtain counts of \n",
    "    each output state. The measurements are performed on an ideal simulator.\n",
    "  \n",
    "    Parameters: \n",
    "    trial_state_circuit (QuantumCircuit): Circuit that prepares the trial state\n",
    "    pauli_product (str): String representation of tensor product of pauli matrices\n",
    "    number_of_shots (int): Number of times measurements should be performed on the trial state\n",
    "    \n",
    "    Returns: \n",
    "    counts (dict): Dictionary of counts of each output state\n",
    "    \"\"\"\n",
    "\n",
    "    measurement_circuit = trial_state_circuit.copy()\n",
    "\n",
    "    qubit = 0\n",
    "    cbit = 0\n",
    "\n",
    "    for pauli_matrix in pauli_product:\n",
    "        if pauli_matrix == 'X':\n",
    "            measurement_circuit.h(qubit)\n",
    "        elif pauli_matrix == 'Y':\n",
    "            measurement_circuit.sdg(qubit)\n",
    "            measurement_circuit.h(qubit)\n",
    "        elif pauli_matrix != 'I' and pauli_matrix != 'Z':\n",
    "            raise ValueError(\"Pauli product should consist only of I, X, Y or Z matrices\")\n",
    "        measurement_circuit.measure(qubit, cbit)\n",
    "        qubit += 1\n",
    "        cbit += 1\n",
    "\n",
    "    backend = Aer.get_backend('qasm_simulator')\n",
    "    job = execute(measurement_circuit, backend, shots = number_of_shots)\n",
    "    result = job.result()\n",
    "    counts = result.get_counts()\n",
    "\n",
    "    return counts"
   ]
  },
  {
   "source": [
    "## Calculating expectation values on the basis of counts\n",
    "\n",
    "The quantum modules give us a set of counts for each basis state i.e.$\\ \\left| 00 \\right\\rangle,\\ \\left| 01 \\right\\rangle,$$\\ \\left| 10 \\right\\rangle\\ and\\ |11\\rangle$. The probability of each state is calculated by dividing the count of each state by the total number of counts.\n",
    "\n",
    "The expectation value is the sum of product of probabilities of each state and their associated eigen value. As we know, the eigen value of state $\\left| 0 \\right\\rangle,\\left| + \\right\\rangle$, and $\\left| i \\right\\rangle$ is $+ 1$ and of state $\\left| 1 \\right\\rangle,\\left| - \\right\\rangle$, and $\\left| - i \\right\\rangle$ is $- 1$. The eigen value of tensor products of these states would then be a product of their individual eigen values. Since every state has been transformed into the Z basis, the eigenvalues of only 4 states need to tbe considered. It comes out to be $+ 1$ for $|00\\rangle$ and $|11\\rangle$ and $- 1\\ $for $|01\\rangle$ and $|10\\rangle$.\n",
    "\n",
    "The formula used for the calculation depends on the Pauli product used but can be generalised to three cases:\n",
    "\n",
    "1.  When the Pauli product is $\\sigma_{I} \\otimes \\sigma_{X}$, $\\sigma_{I} \\otimes \\sigma_{Y}$ or $\\sigma_{I} \\otimes \\sigma_{Z}$:\n",
    "\n",
    "The expectation value in these cases depends only on the second qubit as the first Pauli matrix is $I$. There is no need to create a totally new circuit in these cases and the probabilities obtained from the quantum modules can be used. The eigenvalue of the state is considered as the eigenvalue of state of qubit 2 i.e. states having 1 on the second qubit have $- 1$ and for the states having on the first qubit have $+ 1$. Hence, the expectation value for this case is:\n",
    "\n",
    "$$\\left\\langle \\sigma_{I} \\otimes \\sigma_{i} \\right\\rangle\\  = P_{00}\\left( + 1 \\right) + P_{01}\\left( - 1 \\right) + P_{10}\\left( + 1 \\right) + P_{11}( - 1)\\ where\\ i = X,\\ Y\\ and\\ Z\\ $$\n",
    "\n",
    "2.  When the Pauli product is $\\sigma_{X} \\otimes \\sigma_{I}$, $\\sigma_{Y} \\otimes \\sigma_{I}$ or $\\sigma_{Z} \\otimes \\sigma_{I}$:\n",
    "\n",
    "Similar to the above case, the expectation value would depend only on the first qubit as the second Pauli matrix is $I$. Similar to the above case, the eigenvalue for states having 1 on the first qubit is considered as $- 1$ and for the states having 0 on the first qubit is considered as $+ 1$. Therefore, the expectation value for this case is:\n",
    "\n",
    "$$\\left\\langle \\sigma_{i} \\otimes \\sigma_{I} \\right\\rangle\\  = P_{00}\\left( + 1 \\right) + P_{01}\\left( + 1 \\right) + P_{10}\\left( - 1 \\right) + P_{11}( - 1)\\ where\\ i = X,\\ Y\\ and\\ Z\\ $$\n",
    "\n",
    "3.  When the Pauli product is of the form $\\sigma_{i} \\otimes \\sigma_{j}$ where $i,\\ j = X,\\ Y\\ and\\ Z$:\n",
    "\n",
    "In this case, the eigen value of the entire state is considered. The eigen value of the 4 states was defined initially and is used here. Therefore, the expectation value for this case is:\n",
    "\n",
    "$$\\left\\langle \\sigma_{i} \\otimes \\sigma_{j} \\right\\rangle\\  = P_{00}\\left( + 1 \\right) + P_{01}\\left( - 1 \\right) + P_{10}\\left( - 1 \\right) + P_{11}\\left( + 1 \\right)\\ where\\ i,j = X,\\ Y,\\ Z$$\n",
    "\n",
    "This tutorial<sup>[[7]](https://github.com/DavitKhach/quantum-algorithms-tutorials/blob/master/variational_quantum_eigensolver.ipynb)</sup> by David Khachatryan goes into much depth on the significance of adding extra gates before measurement and calculating expectation values and can be referred for further details."
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_expectation_value(counts, pauli_product):\n",
    "    \"\"\" \n",
    "    Calculates the expectation value of the Pauli product based on the counts of each state and the formula defined above\n",
    "  \n",
    "    Parameters: \n",
    "    counts (dict): Dictionary of counts of each output state\n",
    "    pauli_product (str): String representation of tensor product of pauli matrices\n",
    "    \n",
    "    Returns: \n",
    "    expectation_value (int): Expectation value of the Pauli product based on the counts \n",
    "    \"\"\"\n",
    "\n",
    "    if pauli_product == 'II':\n",
    "        return 1\n",
    "    \n",
    "    if '00' not in counts:\n",
    "        counts['00'] = 0\n",
    "    if '01' not in counts:\n",
    "        counts['01'] = 0\n",
    "    if '10' not in counts:\n",
    "        counts['10'] = 0\n",
    "    if '11' not in counts:\n",
    "        counts['11'] = 0\n",
    "\n",
    "    total_count = counts['00'] + counts['01'] + counts['10'] + counts['11']\n",
    "\n",
    "\n",
    "    # The formulae vary slightly from above as Qiskit has reverse ordering of qubits (little-endian)\n",
    "    # i.e. qubit 1 is qubit 2 and qubit 2 is qubit 1\n",
    "\n",
    "    if pauli_product == 'IX' or pauli_product == 'IY' or pauli_product == 'IZ': \n",
    "        expectation_value = (counts['00'] + counts['01'] - counts['10'] - counts['11']) / total_count\n",
    "    elif pauli_product == 'XI' or pauli_product == 'YI' or pauli_product == 'ZI':\n",
    "        expectation_value = (counts['00'] - counts['01'] + counts['10'] - counts['11']) / total_count\n",
    "    else:\n",
    "        expectation_value = (counts['00'] - counts['01'] - counts['10'] + counts['11']) / total_count\n",
    "\n",
    "    return expectation_value"
   ]
  },
  {
   "source": [
    "## Combining everything to calculate the expectation value\n",
    "\n",
    "All the parts needed to calculate the expectation value of the given Hamiltonian are defined above. We use the coefficients of pauli products (obtained by decomposing the matrix) and set of parameters for the ansatz to call the respective functions and calculate expectation values for each Pauli product. These expectation values are multiplied with their respective coefficients and added to give the expectation value of the given Hamiltonian with respect to the current trial state."
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "def calculate_expectation_value_of_hamiltonian(parameters, coefficient_dict):\n",
    "    \"\"\" \n",
    "    Calculates the expectation value of the hamiltonian using the parameters for trial state circuit\n",
    "    and coefficients of pauli products\n",
    "  \n",
    "    Parameters: \n",
    "    parameters (np.array): List of angles that act as parameters for the trial state circuit \n",
    "    coefficient_dict (dict): Coeffiecients of pauli products obtained by decomposing the hamiltonian\n",
    "    \n",
    "    Returns: \n",
    "    expectation_value_of_hamiltonian (float): Expectation value of the hamiltonian\n",
    "    \"\"\"\n",
    "\n",
    "    trial_state_circuit = create_trial_state_circuit(parameters)\n",
    "\n",
    "    expectation_value = 0\n",
    "\n",
    "    for pauli_product in coefficient_dict:\n",
    "        if abs(coefficient_dict[pauli_product]) > 0:\n",
    "            counts = quantum_module_simulator(trial_state_circuit, pauli_product, 8192)\n",
    "            expectation_value += np.real(coefficient_dict[pauli_product]) * calculate_expectation_value(counts, pauli_product)\n",
    "    \n",
    "    return expectation_value            "
   ],
   "cell_type": "code",
   "metadata": {},
   "execution_count": 9,
   "outputs": []
  },
  {
   "source": [
    "## Optimizing the expectation value of the Hamiltonian\n",
    "\n",
    "Using the above methods, the expectation value with respect to the current trial state is obtained. It needs to be optimised as it may not be the lowest upper bound on the eigenvalue. Classical optimisation methods such as gradient descent or Nelder-Mead Simplex method can be used to optimize the above function and obtain the lowest possible value. In this notebook, Powell's method is used for optimization. I have also shown the lowest eigen value that is calculated using classical methods for comparison.\n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "Classically calculated eigenvalues of given hamiltonian are:  [ 1. -1.  1.  1.]\nLowest eigenvalue calculated using classical methods:  -1.0\nUpper bound on lowest eigenvalue calculated using VQE:  -1.0\n"
    }
   ],
   "source": [
    "given_matrix = np.array([[1,  0,  0,  0], \n",
    "                         [0,  0, -1,  0],\n",
    "                         [0, -1,  0,  0],\n",
    "                         [0,  0,  0,  1]], dtype=np.float64) \n",
    "                         \n",
    "coefficient_dict = decompose_matrix(given_matrix)\n",
    "\n",
    "theta1 = 0\n",
    "parameters = [theta1]\n",
    "tolerance = 1e-5\n",
    "\n",
    "''' Running classical algorithm on the given matrix to find exact value '''\n",
    "eigenvalues = np.linalg.eigvals(given_matrix)\n",
    "lowest_eigenvalue = np.min(eigenvalues)\n",
    "\n",
    "print(\"Classically calculated eigenvalues of given hamiltonian are: \", eigenvalues)\n",
    "print(\"Lowest eigenvalue calculated using classical methods: \", lowest_eigenvalue)\n",
    "\n",
    "result = optimize.minimize(fun=calculate_expectation_value_of_hamiltonian, x0=parameters, \n",
    "                           args=(coefficient_dict), method='Powell', tol=tolerance)\n",
    "\n",
    "print(\"Upper bound on lowest eigenvalue calculated using VQE: \", result.fun)"
   ]
  },
  {
   "source": [
    "As it is seen above, the VQE algorithm is able to return the exact lowest eigenvalue for the given Hamiltonian for the selected ansatz (the classical method may give a small error which can be neglected). "
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "## Observations and Results\n",
    "\n",
    "### Selection of Ansatz\n",
    "As I mentioned earlier, I tried several different ansatze and kept a track of the expectation value of the Hamiltonian obtained using the ansatze. My initial approach was to try several ansatze randomly and observe the effect they have on the expectation value. Some of my ansatze were based on the ones shown in lectures 25 to 27 of the course Introduction to Quantum Computing and Quantum Hardware<sup>[[9]]</sup> by IBM. This paper<sup>[[10]]</sup> by Sukin Sim et. al. helped me visualize the expressibility of different circuits and guided me in choosing a good ansatz.\n",
    "\n",
    "I obtained the expectation value for each ansatz thrice to show the variation in values due to random seed of the noiseless simulator and the convergence of the optimiser. The data for each ansatz is stored in a file name experimental_data.csv and the corresponding circuit diagram (parameters of the circuit have values for the last iteration) for each ansatz is stored in the circuits directory. Due to lack of space, I do not display the description and plot of each ansatz below. In all cases, the number of shots (i.e. total count) for each circuit is fixed at 8192.\n",
    "\n",
    "[9]: https://qiskit.org/learn/intro-qc-qh/\n",
    "[10]: https://arxiv.org/abs/1905.10876\n",
    "\n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "experimental_data = pd.read_csv(\"experimental_data.csv\")\n",
    "experimental_data_without_description = experimental_data.loc[:, 'ID':]\n",
    "\n",
    "experimental_data_without_description.head(9)"
   ],
   "cell_type": "code",
   "metadata": {
    "tags": []
   },
   "execution_count": 11,
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "           ID  Number of gates  Number of parameters  Expectation Value 1  \\\n0     Ansatz1                4                     2            -0.002686   \n1     Ansatz2                4                     4            -0.002563   \n2     Ansatz3                4                     2            -0.014160   \n3     Ansatz4                4                     4            -0.010742   \n4     Ansatz5                9                     2            -0.261841   \n5     Ansatz6               14                     2            -0.061035   \n6     Ansatz7                9                     4            -0.505371   \n7     Ansatz8               14                     4            -0.947998   \n8  AnsatzHint                3                     1            -1.000000   \n\n   Expectation Value 2  Expectation Value 3  \n0            -0.011108            -0.006226  \n1             0.000244            -0.009277  \n2            -0.013062            -0.014771  \n3            -0.015503             0.000122  \n4            -0.269043            -0.266235  \n5            -0.062500            -0.057007  \n6            -0.511108            -0.499512  \n7            -0.509644            -0.945435  \n8            -1.000000            -1.000000  ",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>ID</th>\n      <th>Number of gates</th>\n      <th>Number of parameters</th>\n      <th>Expectation Value 1</th>\n      <th>Expectation Value 2</th>\n      <th>Expectation Value 3</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>0</th>\n      <td>Ansatz1</td>\n      <td>4</td>\n      <td>2</td>\n      <td>-0.002686</td>\n      <td>-0.011108</td>\n      <td>-0.006226</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>Ansatz2</td>\n      <td>4</td>\n      <td>4</td>\n      <td>-0.002563</td>\n      <td>0.000244</td>\n      <td>-0.009277</td>\n    </tr>\n    <tr>\n      <th>2</th>\n      <td>Ansatz3</td>\n      <td>4</td>\n      <td>2</td>\n      <td>-0.014160</td>\n      <td>-0.013062</td>\n      <td>-0.014771</td>\n    </tr>\n    <tr>\n      <th>3</th>\n      <td>Ansatz4</td>\n      <td>4</td>\n      <td>4</td>\n      <td>-0.010742</td>\n      <td>-0.015503</td>\n      <td>0.000122</td>\n    </tr>\n    <tr>\n      <th>4</th>\n      <td>Ansatz5</td>\n      <td>9</td>\n      <td>2</td>\n      <td>-0.261841</td>\n      <td>-0.269043</td>\n      <td>-0.266235</td>\n    </tr>\n    <tr>\n      <th>5</th>\n      <td>Ansatz6</td>\n      <td>14</td>\n      <td>2</td>\n      <td>-0.061035</td>\n      <td>-0.062500</td>\n      <td>-0.057007</td>\n    </tr>\n    <tr>\n      <th>6</th>\n      <td>Ansatz7</td>\n      <td>9</td>\n      <td>4</td>\n      <td>-0.505371</td>\n      <td>-0.511108</td>\n      <td>-0.499512</td>\n    </tr>\n    <tr>\n      <th>7</th>\n      <td>Ansatz8</td>\n      <td>14</td>\n      <td>4</td>\n      <td>-0.947998</td>\n      <td>-0.509644</td>\n      <td>-0.945435</td>\n    </tr>\n    <tr>\n      <th>8</th>\n      <td>AnsatzHint</td>\n      <td>3</td>\n      <td>1</td>\n      <td>-1.000000</td>\n      <td>-1.000000</td>\n      <td>-1.000000</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "metadata": {},
     "execution_count": 11
    }
   ]
  },
  {
   "source": [
    "The following observations and results are made based on the above experimental data:\n",
    "\n",
    "1. In case of Ansatz1, Ansatz2, Ansatz3 and Ansatz4, the expectation values that are obtained are 10<sup>4</sup> to 10<sup>2</sup> orders of magnitude away from the lowest eigenvalue. This may be due to 2 reasons: Firstly, the ansatze do not have many number of gates and secondly, they do not have entangling capability and hence, are limited in the number of states they represent. So, there is a good possibility that the subspace they represent does not have any vector near the required eigenvector.\n",
    "\n",
    "Note 1: Entanglement was introduced in all circuits from Ansatz5, which increases the expressibility of the ansatz.\n",
    "\n",
    "2. In case of Ansatz 5, the number of gates have been increased and entanglement has been introduced. So, the circuit is more expressible. But, the subspace represented by the ansatz still does not contain the exact state represented by the required eigen vector. Hence, we are still an order of magnitude away from the solution.\n",
    "\n",
    "3. In case of Ansatz 6, the number of gates have been further increased but this worsens the expectation value. This might be due the fact that we have only 2 parameters and adding more gates makes the task of the optimiser difficult without giving a better solution.\n",
    "\n",
    "Note 2: Adding more gates and parameters increases the time taken to calculate the eigenvalue. This delay in execution is visible while obtaining values for Ansatze 5 through 8.\n",
    "\n",
    "4. In case of Ansatz 7, we add more parameters (4) to the circuit while keeping the number of gates constant (9). This increases the expressibility further and we get a much better upper bound. \n",
    "\n",
    "5. In case of Ansatz 8, we increase the gates again but this time we have 4 parameters instead of 2. This is the best circuit I made with the closest upper bound. Expectation value 2 for this ansatz is an anomaly when compared to Expectation value 1 and Expectation value 3. This is the problem that plagues larger circuits. As the sub space represented by them increases, there may be local minima that occur which cause the optimizer to converge at such values. Repeatedly performing the experiments at different initial points is one of the solutions to this problem. \n",
    "\n",
    "6. I had a look at the ansatz given in the hint and applied it here. It had only 3 gates and 1 parameter so I was suspicious. But surprisingly, it gave the exact eigenvalue of -1 each time I executed it without taking long to execute. Later, I deduced that the subspace represented by AnsatzHint contains the exact eigenvector that we need and varying the parameter would always lead us to the exact solution. Since, the circuit is small, it is not expressible but it is one of the best ansatz for this case.\n",
    "\n",
    "After analysing all these observations, it is seen that we can get a good (or exact) upper bound on the lowest eigen value if we know the subspace in which the corresponding eigenvector lies. We can design very shallow circuits which sweep that subspace to get the eigenvalue. This approach cannot be generalised for large and random hamiltonians as we would not know the subspace in which the eigenvector lies. In such cases, we need expressible circuits and generalised circuits like Ansatz8 which would take more time but give us a good upper bound for most cases."
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "### Introducing noise into the circuit\n",
    "\n",
    "I ran the circuit having AnsatzHint which gave an exact eigen value of -1 on the noiseless simulator on a noisy IBM quantum computer (ibmq_15_melbourne). Since, there is queue for executing on real devices, I executed the circuit once with a parameter of pi (This parameter is the one that gives exact eigen value)."
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "provider = IBMQ.load_account() # Load credentials for IBMQ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quantum_module_hardware(trial_state_circuit, pauli_product, number_of_shots):\n",
    "    \"\"\" \n",
    "    This is a generalized function that adds transformations and performs measurements on the trial\n",
    "    state based on the pauli matrices. The measurements are performed repeatedly to obtain counts of \n",
    "    each output state. The measurements are performed on an real hardware.\n",
    "  \n",
    "    Parameters: \n",
    "    trial_state_circuit (QuantumCircuit): Circuit that prepares the trial state\n",
    "    pauli_product (str): String representation of tensor product of pauli matrices\n",
    "    number_of_shots (int): Number of times measurements should be performed on the trial state\n",
    "    \n",
    "    Returns: \n",
    "    counts (dict): Dictionary of counts of each output state\n",
    "    \"\"\"\n",
    "\n",
    "    measurement_circuit = trial_state_circuit.copy()\n",
    "\n",
    "    qubit = 0\n",
    "    cbit = 0\n",
    "\n",
    "    for pauli_matrix in pauli_product:\n",
    "        if pauli_matrix == 'X':\n",
    "            measurement_circuit.h(qubit)\n",
    "        elif pauli_matrix == 'Y':\n",
    "            measurement_circuit.sdg(qubit)\n",
    "            measurement_circuit.h(qubit)\n",
    "        elif pauli_matrix != 'I' and pauli_matrix != 'Z':\n",
    "            raise ValueError(\"Pauli product should consist only of I, X, Y or Z matrices\")\n",
    "        measurement_circuit.measure(qubit, cbit)\n",
    "        qubit += 1\n",
    "        cbit += 1\n",
    "\n",
    "    backend = provider.backends.ibmq_16_melbourne\n",
    "    job = execute(measurement_circuit, backend, shots = number_of_shots)\n",
    "    result = job.result()\n",
    "    counts = result.get_counts()\n",
    "\n",
    "    return counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_expectation_value_of_hamiltonian_hardware(parameters, coefficient_dict):\n",
    "    \"\"\" \n",
    "    Calculates the expectation value of the hamiltonian using the parameters for trial state circuit\n",
    "    and coefficients of pauli products using real hardware\n",
    "  \n",
    "    Parameters: \n",
    "    parameters (np.array): List of angles that act as parameters for the trial state circuit \n",
    "    coefficient_dict (dict): Coeffiecients of pauli products obtained by decomposing the hamiltonian\n",
    "    \n",
    "    Returns: \n",
    "    expectation_value_of_hamiltonian (float): Expectation value of the hamiltonian\n",
    "    \"\"\"\n",
    "\n",
    "    trial_state_circuit = create_trial_state_circuit(parameters)\n",
    "\n",
    "    expectation_value = 0\n",
    "\n",
    "    for pauli_product in coefficient_dict:\n",
    "        if abs(coefficient_dict[pauli_product]) > 0:\n",
    "            counts = quantum_module_hardware(trial_state_circuit, pauli_product, 8192)\n",
    "            expectation_value += np.real(coefficient_dict[pauli_product]) * calculate_expectation_value(counts, pauli_product)\n",
    "    \n",
    "    return expectation_value "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "Classically calculated eigenvalues of given hamiltonian are:  [ 1. -1.  1.  1.]\nLowest eigenvalue calculated using classical methods:  -1.0\n"
    },
    {
     "output_type": "error",
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-14-23220d33082c>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     17\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"Lowest eigenvalue calculated using classical methods: \"\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlowest_eigenvalue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     18\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 19\u001b[1;33m \u001b[0mexpectation_value\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mcalculate_expectation_value_of_hamiltonian_hardware\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mparameters\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcoefficient_dict\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     20\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     21\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"Upper bound on lowest eigenvalue calculated using VQE: \"\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mexpectation_value\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-13-6a31c3decd59>\u001b[0m in \u001b[0;36mcalculate_expectation_value_of_hamiltonian_hardware\u001b[1;34m(parameters, coefficient_dict)\u001b[0m\n\u001b[0;32m     18\u001b[0m     \u001b[1;32mfor\u001b[0m \u001b[0mpauli_product\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mcoefficient_dict\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     19\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mabs\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcoefficient_dict\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mpauli_product\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m>\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 20\u001b[1;33m             \u001b[0mcounts\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mquantum_module_hardware\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtrial_state_circuit\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mpauli_product\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m8192\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     21\u001b[0m             \u001b[0mexpectation_value\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mreal\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcoefficient_dict\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mpauli_product\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mcalculate_expectation_value\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcounts\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mpauli_product\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     22\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-12-d5008e52518f>\u001b[0m in \u001b[0;36mquantum_module_hardware\u001b[1;34m(trial_state_circuit, pauli_product, number_of_shots)\u001b[0m\n\u001b[0;32m     33\u001b[0m     \u001b[0mbackend\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mprovider\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbackends\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mibmq_16_melbourne\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     34\u001b[0m     \u001b[0mjob\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mexecute\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmeasurement_circuit\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbackend\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mshots\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnumber_of_shots\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 35\u001b[1;33m     \u001b[0mresult\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mjob\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mresult\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     36\u001b[0m     \u001b[0mcounts\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mresult\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_counts\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     37\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mc:\\Users\\Rohit Pai\\Quantum Computing\\Variational Quantum Eigensolver\\vqe\\lib\\site-packages\\qiskit\\providers\\ibmq\\job\\ibmqjob.py\u001b[0m in \u001b[0;36mresult\u001b[1;34m(self, timeout, wait, partial, refresh)\u001b[0m\n\u001b[0;32m    263\u001b[0m         \"\"\"\n\u001b[0;32m    264\u001b[0m         \u001b[1;31m# pylint: disable=arguments-differ\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 265\u001b[1;33m         if not self._wait_for_completion(timeout=timeout, wait=wait,\n\u001b[0m\u001b[0;32m    266\u001b[0m                                          required_status=(JobStatus.DONE,)):\n\u001b[0;32m    267\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_status\u001b[0m \u001b[1;32mis\u001b[0m \u001b[0mJobStatus\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mCANCELLED\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mc:\\Users\\Rohit Pai\\Quantum Computing\\Variational Quantum Eigensolver\\vqe\\lib\\site-packages\\qiskit\\providers\\ibmq\\job\\ibmqjob.py\u001b[0m in \u001b[0;36m_wait_for_completion\u001b[1;34m(self, timeout, wait, required_status, status_queue)\u001b[0m\n\u001b[0;32m    814\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    815\u001b[0m         \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 816\u001b[1;33m             status_response = self._api.job_final_status(\n\u001b[0m\u001b[0;32m    817\u001b[0m                 self.job_id(), timeout=timeout, wait=wait, status_queue=status_queue)\n\u001b[0;32m    818\u001b[0m         \u001b[1;32mexcept\u001b[0m \u001b[0mUserTimeoutExceededError\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mc:\\Users\\Rohit Pai\\Quantum Computing\\Variational Quantum Eigensolver\\vqe\\lib\\site-packages\\qiskit\\providers\\ibmq\\api\\clients\\account.py\u001b[0m in \u001b[0;36mjob_final_status\u001b[1;34m(self, job_id, timeout, wait, status_queue)\u001b[0m\n\u001b[0;32m    330\u001b[0m             \u001b[0mstart_time\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtime\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtime\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    331\u001b[0m             \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 332\u001b[1;33m                 status_response = self._job_final_status_websocket(\n\u001b[0m\u001b[0;32m    333\u001b[0m                     job_id=job_id, timeout=timeout, status_queue=status_queue)\n\u001b[0;32m    334\u001b[0m             \u001b[1;32mexcept\u001b[0m \u001b[0mWebsocketTimeoutError\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0mex\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mc:\\Users\\Rohit Pai\\Quantum Computing\\Variational Quantum Eigensolver\\vqe\\lib\\site-packages\\qiskit\\providers\\ibmq\\api\\clients\\account.py\u001b[0m in \u001b[0;36m_job_final_status_websocket\u001b[1;34m(self, job_id, timeout, status_queue)\u001b[0m\n\u001b[0;32m    382\u001b[0m             \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    383\u001b[0m                 \u001b[1;32mraise\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 384\u001b[1;33m         return loop.run_until_complete(\n\u001b[0m\u001b[0;32m    385\u001b[0m             self.client_ws.get_job_status(job_id, timeout=timeout, status_queue=status_queue))\n\u001b[0;32m    386\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mc:\\Users\\Rohit Pai\\Quantum Computing\\Variational Quantum Eigensolver\\vqe\\lib\\site-packages\\nest_asyncio.py\u001b[0m in \u001b[0;36mrun_until_complete\u001b[1;34m(self, future)\u001b[0m\n\u001b[0;32m     87\u001b[0m                 \u001b[0mf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_log_destroy_pending\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mFalse\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     88\u001b[0m             \u001b[1;32mwhile\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[0mf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdone\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 89\u001b[1;33m                 \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_run_once\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     90\u001b[0m                 \u001b[1;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_stopping\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     91\u001b[0m                     \u001b[1;32mbreak\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mc:\\Users\\Rohit Pai\\Quantum Computing\\Variational Quantum Eigensolver\\vqe\\lib\\site-packages\\nest_asyncio.py\u001b[0m in \u001b[0;36m_run_once\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m    111\u001b[0m             \u001b[1;32melse\u001b[0m \u001b[0mmin\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmax\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mscheduled\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_when\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mnow\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m10\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mscheduled\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m\\\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    112\u001b[0m             \u001b[1;32melse\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 113\u001b[1;33m         \u001b[0mevent_list\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_selector\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mselect\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    114\u001b[0m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_process_events\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mevent_list\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    115\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mc:\\users\\rohit pai\\appdata\\local\\programs\\python\\python38\\lib\\selectors.py\u001b[0m in \u001b[0;36mselect\u001b[1;34m(self, timeout)\u001b[0m\n\u001b[0;32m    321\u001b[0m         \u001b[0mready\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    322\u001b[0m         \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 323\u001b[1;33m             \u001b[0mr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mw\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0m_\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_select\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_readers\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_writers\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtimeout\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    324\u001b[0m         \u001b[1;32mexcept\u001b[0m \u001b[0mInterruptedError\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    325\u001b[0m             \u001b[1;32mreturn\u001b[0m \u001b[0mready\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mc:\\users\\rohit pai\\appdata\\local\\programs\\python\\python38\\lib\\selectors.py\u001b[0m in \u001b[0;36m_select\u001b[1;34m(self, r, w, _, timeout)\u001b[0m\n\u001b[0;32m    312\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[0msys\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplatform\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;34m'win32'\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    313\u001b[0m         \u001b[1;32mdef\u001b[0m \u001b[0m_select\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mw\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0m_\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtimeout\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mNone\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 314\u001b[1;33m             \u001b[0mr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mw\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mselect\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mselect\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mw\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mw\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtimeout\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    315\u001b[0m             \u001b[1;32mreturn\u001b[0m \u001b[0mr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mw\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    316\u001b[0m     \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "given_matrix = np.array([[1,  0,  0,  0], \n",
    "                         [0,  0, -1,  0],\n",
    "                         [0, -1,  0,  0],\n",
    "                         [0,  0,  0,  1]], dtype=np.float64) \n",
    "                         \n",
    "coefficient_dict = decompose_matrix(given_matrix)\n",
    "\n",
    "theta1 = pi\n",
    "parameters = [theta1]\n",
    "tolerance = 1e-5\n",
    "\n",
    "''' Running classical algorithm on the given matrix to find exact value '''\n",
    "eigenvalues = np.linalg.eigvals(given_matrix)\n",
    "lowest_eigenvalue = np.min(eigenvalues)\n",
    "\n",
    "print(\"Classically calculated eigenvalues of given hamiltonian are: \", eigenvalues)\n",
    "print(\"Lowest eigenvalue calculated using classical methods: \", lowest_eigenvalue)\n",
    "\n",
    "expectation_value = calculate_expectation_value_of_hamiltonian_hardware(parameters, coefficient_dict)\n",
    "\n",
    "print(\"Upper bound on lowest eigenvalue calculated using VQE: \", expectation_value)"
   ]
  },
  {
   "source": [
    "I wasn't able to execute the circuit within time but previously I have achieved an expectation value of -0.9023 with theta1 as 0. Noise played a big role in the deviation and lead to errors. This can be fixed only when we have better devices that are resistant to noise and errors. There is constant research being done for this and we have a bright future ahead of us. "
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "## Conclusion\n",
    "\n",
    "I want to conclude by saying that it was a great experience working on NISQ algorithm and learning it in depth. The main outcomes for me were that choosing an ansatz for such circuits is an art and it would take a lot of time to develop an intuition in it. I also learnt that there are several factors that can affect the outcome of a variational circuit such as cost of the circuit, expressibility of the circuit, number of parameters, noise of the device, etc. and striking the right balance is important to get a good result. "
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "## References\n",
    "\n",
    "1. https://www.mustythoughts.com/variational-quantum-eigensolver-explained\n",
    "2. https://arxiv.org/abs/1304.3061\n",
    "3. https://quantumcomputing.stackexchange.com/questions/11899/example-of-hamiltonian-decomposition-into-pauli-matrices\n",
    "4. https://quantumcomputing.stackexchange.com/questions/8725/can-arbitrary-matrices-be-decomposed-using-the-pauli-basis\n",
    "5. https://quantumcomputing.stackexchange.com/questions/6882/decomposition-of-a-matrix-in-the-pauli-basis\n",
    "6. https://en.wikipedia.org/wiki/Ansatz\n",
    "7. https://github.com/DavitKhach/quantum-algorithms-tutorials/blob/master/variational_quantum_eigensolver.ipynb\n",
    "8. https://michaelgoerz.net/notes/decomposing-two-qubit-hamiltonians-into-pauli-matrices.html\n",
    "9. [Introduction to Quantum Computing and Quantum Hardware - Lecture 25 to 27](https://qiskit.org/learn/intro-qc-qh/)\n",
    "10. Sim, Sukin, Peter D. Johnson, and Alán Aspuru‐Guzik. “Expressibility and Entangling Capability of Parameterized Quantum Circuits for Hybrid Quantum‐Classical Algorithms.” Advanced Quantum Technologies 2.12 (2019): 1900070. Crossref. Web.\n",
    "\n",
    "## Additional Resources used for studying Variational Quantum Eigensolver\n",
    "    \n",
    "1. [Framework agnostic VQE tutorial by Alexander Soare](https://github.com/alexander-soare/framework-agnostic-vqe-tutorial)"
   ],
   "cell_type": "markdown",
   "metadata": {}
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "python_defaultSpec_1600878381426",
   "display_name": "Python 3.8.5 64-bit ('vqe': venv)"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}