notebooks/advanced_algorithms/Quantum_Principal_Component_Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "84f5ead3-4910-4b85-b1f6-60187001d382",
   "metadata": {
    "tags": []
   },
   "source": [
    "$$\\newcommand{\\ket}[1]{\\left|{#1}\\right\\rangle}$$\n",
    "$$\\newcommand{\\bra}[1]{\\left\\langle{#1}\\right|}$$\n",
    "$$\\newcommand{\\braket}[2]{\\left\\langle{#1}\\middle|{#2}\\right\\rangle}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "152d5aca-e16e-4341-a714-13c0986d9f83",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Quantum Principal Component Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6667554e-a2ff-4238-91d4-5b005ffeb977",
   "metadata": {},
   "source": [
    "Recent areas of computing such as Machine Learning had to face different obstacles to achieve their successes. An example of these problems is the dimensionality of the data that must be processed in our current digital age. With a highly connected world, processing data with high dimensionality is a difficult task due to the resource demand, even for existing supercomputers. For this reason, different techniques have been created to subtract the most relevant information from data sets with high dimensionality, such as the **Principal Component Analysis** (PCA) technique. In a nutshell, PCA seeks to reduce the dimensionality of a data set, while preserving as much variability (i.e., statistical information) as possible [1]."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6eeb39cd-8541-4d7c-a6f2-b0bed05c76c0",
   "metadata": {},
   "source": [
    "On the other hand, all these techniques to improve data processing and extract the most relevant information have been conceived under the paradigm of classical computing. Classical computing receives its name to differentiate it from the paradigm of quantum computing which leverages the laws of quantum mechanics to take advantage of new properties that are not present in classical computing [2]. Therefore, these and other machine learning techniques are now being explored using this new paradigm (this new branch is known as Quantum Machine Learning [3]) since these quantum devices can potentially store and process exponentially more information than their classical counterparts. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4300e8a1-aa5c-438b-8ad4-1103e238698f",
   "metadata": {},
   "source": [
    "The purpose of this notebook is to implement PCA on a quantum processor by using Amazon Braket. We will use the well-known scenario for PCA of housing prices in the United States (the same one used by [7] from 1996). The idea behind this scenario is to show how a single variable (i.e. price of the house) can be function of many other features such as the number of bedrooms, number of bathrooms, if the kitchen has been recently remodeled, if it has a balcony, square footage, if it has sea view, lot size, amount of years of the construction, location of the house (neighborhood), if it has attic or basement, number of parking spots, etc. Due to the amount of features, it is important to reduce the data to the features that capture the largest variance in the data. Therefore, we are going to use the PCA technique which seeks to determine which features capture the largest variance on the data [1]."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c874eb1-fac6-44d5-ba96-09bbf5071b78",
   "metadata": {},
   "source": [
    "For our case, let's consider the number of bedrooms (first feature) and the square footage (second feature) of several houses for sale in Los Alamos (following the same example as [4]). Here is the raw data, taken from __[Zillow](http://www.zillow.com)__, for 15 houses:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d2b98db-b496-4767-9df9-07c015eb71d8",
   "metadata": {
    "tags": []
   },
   "source": [
    "Number of bedrooms:\n",
    "\\begin{align}\n",
    "X_1 = [4,3,4,4,3,3,3,3,4,4,4,5,4,3,4]\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9e4b171-09bf-4909-a84a-00971b18dba4",
   "metadata": {},
   "source": [
    "Square footage:\n",
    "\n",
    "\\begin{align}\n",
    "X_2 = [3028, 1365, 2726, 2538, 1318, 1693, 1412, 1632, 2875, 3564, 4412, 4444, 4278, 3064, 3857]\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b423903c-a6a4-4918-9770-278bb56a73a4",
   "metadata": {},
   "source": [
    "## Principal Component Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4b0375d-7107-48e0-bac8-7c83d83c816d",
   "metadata": {},
   "source": [
    "PCA has three main steps. The first one is called standarization which means we must adjust the data to an appropiate scale (the idea is to range the continuous initial variables so that each one of them contributes equally to the analysis. To learn more about this point you can see the following [link](https://erdogant.github.io/pca/pages/html/Algorithm.html). Then, we're going to divide each feature ($X_1$,$X_2$) by the standard deviation ($\\sigma$) so we can have both features on the same scale, and then we're going to substract off the mean of both features (we substract the mean because we want the features to have the features of a standard normal distribution with mean of zero): "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "12dfd4b3-4628-4dc8-85d0-e7f97d5f0d40",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "X_1 = [4,3,4,4,3,3,3,3,4,4,4,5,4,3,4]\n",
    "X_2 = [3028,1365,2726,2538,1318,1693,1412,1632,2875,3564,4412,4444,4278,3064,3857]\n",
    "X_1 = (X_1 - np.average(X_1)) / np.std(X_1)\n",
    "X_2 = (X_2 - np.average(X_2)) / np.std(X_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ea084e16-fc60-4791-b609-e05eda4803f8",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The rescaled feature vectors are\n",
      "X_1 =  [ 0.55901699 -1.11803399  0.55901699  0.55901699 -1.11803399 -1.11803399\n",
      " -1.11803399 -1.11803399  0.55901699  0.55901699  0.55901699  2.23606798\n",
      "  0.55901699 -1.11803399  0.55901699]\n",
      "X_2 =  [ 0.19474995 -1.31677383 -0.07974204 -0.25061785 -1.35949278 -1.01865008\n",
      " -1.27405488 -1.07409383  0.05568612  0.68192779  1.45268675  1.481772\n",
      "  1.33089229  0.22747085  0.94823955]\n"
     ]
    }
   ],
   "source": [
    "print('The rescaled feature vectors are')\n",
    "print('X_1 = ', X_1)\n",
    "print('X_2 = ', X_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fb7ac73-582e-4ac2-839e-8509acb1097f",
   "metadata": {},
   "source": [
    "The next step for PCA is calculating the covariance matrix (variance is a measure of dispersion and can be defined as the spread of data from the mean of the given dataset. Covariance is calculated between two variables and is used to measure how the two variables vary together. To learn more about these terms you can see the following [link](https://www.cuemath.com/algebra/covariance-matrix/)) which is defined by (where $\\sigma[X]$ is the standard deviation of $X$, and note the $cov[X,X]=\\sigma^2[X]=var(X)$ which means the covariance of variable with itself is the variance of the variable): \n",
    "\n",
    "\\begin{align}\n",
    "\\Sigma = \n",
    "\\begin{pmatrix}\n",
    "var[X_1] & cov[X_1 * X_2] \\\\\n",
    "cov[X_2 * X_1] & var[X_2]\n",
    "\\end{pmatrix}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18447c23-a870-4f91-97ef-3c007f16a4c8",
   "metadata": {},
   "source": [
    "We're going to use pandas library on Python to calculate the covariance matrix we need:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f483f3f2-8fc4-431c-b53d-eb3cfcd552a5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df = pd.DataFrame(\n",
    "    {'X_1': X_1,\n",
    "    'X_2': X_2}\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "820fa8db-9d4f-41a9-a194-5972eda10251",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          X_1       X_2\n",
      "X_1  1.071429  0.874147\n",
      "X_2  0.874147  1.071429\n"
     ]
    }
   ],
   "source": [
    "sigma = df.cov()\n",
    "print(sigma)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8537dee3-3f4c-4e99-a9c9-1ea4f8d3df46",
   "metadata": {},
   "source": [
    "Then, our covariance matrix is:\n",
    "\n",
    "\\begin{align}\n",
    "\\Sigma = \n",
    "\\begin{pmatrix}\n",
    "1.071429 & 0.874147 \\\\\n",
    "0.874147 & 1.071429\n",
    "\\end{pmatrix}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21794c07-80d5-4f8e-815b-8e423d72a941",
   "metadata": {},
   "source": [
    "Now, the final step is to find the eigenvalues and the eigenvectors of the covariance matrix so we can find the Principal Components of the data. According to the number of the eigenvalues we know how much variance is being captured by each principal component. In our 2-feature case, we would keep the component with the highest eigenvalue. We can calculate using classical computation these eigenvalues with the linalg library on Python of Numpy (remember that we are going to diagonalize the covariance matrix to find the principal components of our technique):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c1add6af-dc9d-499f-b3c1-2bbda377450c",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sigma_eigenvalues:  [0.19728173 1.94557541]\n",
      "sigma_eigenvectors:  [[-0.70710678 -0.70710678]\n",
      " [ 0.70710678 -0.70710678]]\n"
     ]
    }
   ],
   "source": [
    "sigma_eigenvalues, sigma_eigenvectors = np.linalg.eig(sigma)\n",
    "print('sigma_eigenvalues: ', sigma_eigenvalues)\n",
    "print('sigma_eigenvectors: ', sigma_eigenvectors)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d729e2c0-0c65-47f0-b695-cf264c73cc2f",
   "metadata": {},
   "source": [
    "Then, our eigenvalues (calculated with classical computation) are:\n",
    "\n",
    "\\begin{align}\n",
    "e_1 = 1.945575, e_2 = 0.197281\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9690b7a-cb6f-435e-984a-a7a2afec309d",
   "metadata": {},
   "source": [
    "The reason we are calculating the eigenvalues and eigenvectors is because the principal components are eigenvectors of the data's covariance matrix (this can be mathematically shown but it is out of the scope of this notebook. To learn more information related to eigenvectors and eigenvalues on PCA see [this blog](https://medium.com/@dareyadewumi650/understanding-the-role-of-eigenvectors-and-eigenvalues-in-pca-dimensionality-reduction-10186dad0c5c)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00e0a0bc-2597-4bfe-b799-d5961e60fe42",
   "metadata": {},
   "source": [
    "## Quantum algorithm to compute PCA"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55a9cefa-625b-4ee8-84dd-2fa20d9e3fe5",
   "metadata": {},
   "source": [
    "We're going to implement the same quantum algorithm proposed by [4] which has the following steps:\n",
    "\n",
    "    1. Classical pre-processing.\n",
    "    2. State preparation.\n",
    "    3. Purity calculation.\n",
    "    4. Classical post-processing."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27db4a31-9a23-4d86-9094-56f2065d25d8",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 1. Classical pre-processing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b169f04-93a6-4c1a-88fa-d8edaa13da15",
   "metadata": {},
   "source": [
    "Before talking about the classical pre-processing step, it's important to understand what a density matrix is. Density matrices are used to represent quantum states when they're not pure states but mixed states. If we don't know the exact state where our quantum system is but rather we know that it can be in one of the $M$ states, $\\ket{\\psi_i}$, each with a probability of $p_i$. Then, we can define the density matrix for the quantum system to be: \n",
    "\n",
    "\\begin{align}\n",
    "\\rho = \\sum_{i=1}^{M} p_i \\ket{\\psi_i} \\bra{\\psi_i}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ee1de6e-bee8-4872-b726-45460da75070",
   "metadata": {},
   "source": [
    "Note that the previous density matrix definition gives the same result when we know the quantum system is in a specific state with $p=1$. Now, it's important to note also that from the desnity matrix definition it can be seen that: \n",
    "\n",
    "    i) it's positive semi-definite\n",
    "    ii) it has unit trace\n",
    "    \n",
    "In fact, **any matrix that satisfies these two properties can be interpreted as a density matrix** [4] (more details about density matrix interpretation on the reference). This important fact is how we are going to encode the classical information of the covariance matrix into a quantum state. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6c9e80f-2c26-455c-ba57-3a39d78711e5",
   "metadata": {},
   "source": [
    "Then, our density matrix from our covariance matrix in this case of two features would be (after normalization): \n",
    "\n",
    "\\begin{align}\n",
    "\\rho = \\frac{\\Sigma}{Tr(\\Sigma)}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "b50ca88e-0571-43a3-ba9d-5ee49bf686de",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.5        0.40793519]\n",
      " [0.40793519 0.5       ]]\n"
     ]
    }
   ],
   "source": [
    "rho = sigma / np.trace(sigma)\n",
    "print(np.array(rho))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bcda7a0-463a-41cc-ae9e-de168a25f6fb",
   "metadata": {},
   "source": [
    "Our density matrix $\\rho$ would be:\n",
    "\n",
    "\\begin{align}\n",
    "\\rho = \n",
    "\\begin{pmatrix}\n",
    "0.5 & 0.407935 \\\\\n",
    "0.407935 & 0.5\n",
    "\\end{pmatrix}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe97dae6-f8bf-4ef5-99df-ec87e9a4f6f1",
   "metadata": {},
   "source": [
    "## 2. State preparation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "998f743c-af27-4b32-a427-4f5e71a23f66",
   "metadata": {},
   "source": [
    "Before the state preparation step, we must talk about the purification process. Density matrices are a more useful way to represent quantum systems since we can have both pure and mixed states, whereas state vectors can only represent pure states. However, even a mixed state can be seen as a part of a larger system that is in a pure state. This process of converting a mixed state into a pure state of an enlarged system is called **purification** [4]."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3407992c-548d-4cdc-b310-fc27285551ff",
   "metadata": {
    "tags": []
   },
   "source": [
    "Recalling the definition of the density matrix, and by using the [Schrödinger–HJW theorem](https://en.wikipedia.org/wiki/Schr%C3%B6dinger%E2%80%93HJW_theorem) it can be shown that any density matrix $\\rho$ can be purified, which means it can be seen as the partial trace of a pure state defined in a larger Hilbert space. In other words, it's always possible to find a larger Hilbert space $\\mathcal{H_a}$ with a pure state $\\ket{\\psi_{sa}} \\in \\mathcal{H_s} \\otimes \\mathcal{H_a}$ such that $\\rho = Tr_a(\\ket{\\psi_{sa}}\\bra{\\psi_{sa}})$, and those states satisfy (for some orthonormal basis {$a_i$}): \n",
    "\n",
    "\\begin{align}\n",
    "\\ket{\\psi_{sa}} = \\sum_i \\sqrt{p_i} \\ket{\\phi_i} \\otimes \\ket{a_i}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b50865a1-f0eb-4854-83d4-d99856bedc44",
   "metadata": {},
   "source": [
    "For the state preparation process, it's important to note that in our case of only two features, $\\Sigma$ and $\\rho$, are $2x2$ matrices (one qubit states), then $\\rho$, can be purified to a pure state $\\ket{\\psi}$ on two qubits. Rigorously, we should design a state preparation circuit for this purification process which is made on [4] with details but for this notebook we directly calculate it to focus on the quantum PCA process."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0b6f1b1-789a-47ba-9427-527902881e9d",
   "metadata": {},
   "source": [
    "First, we find the eigenvalues and eigenvector of our density matrix $\\rho$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "0284edbe-2534-4a2c-be01-082299253376",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rho_eig_val:  [0.09206481 0.90793519]\n",
      "rho_eig_vec:  [[-0.70710678 -0.70710678]\n",
      " [ 0.70710678 -0.70710678]]\n"
     ]
    }
   ],
   "source": [
    "rho_eig_val, rho_eig_vec = np.linalg.eig(rho)\n",
    "print('rho_eig_val: ', rho_eig_val)\n",
    "print('rho_eig_vec: ', rho_eig_vec)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b5e9c7e-aac7-4a25-815d-d355e68db378",
   "metadata": {},
   "source": [
    "Then, our eigenvectors (after normalization process) would be:\n",
    "\n",
    "\\begin{align}\n",
    "\\vec{v_1} = \n",
    "\\begin{pmatrix}\n",
    "0.707106 \\\\\n",
    "-0.707106\n",
    "\\end{pmatrix},\n",
    "\\vec{v_2} = \n",
    "\\begin{pmatrix}\n",
    "-0.707106 \\\\\n",
    "-0.707106\n",
    "\\end{pmatrix}\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94e01cc4-ebfc-4e46-9ee8-14801198dbd5",
   "metadata": {},
   "source": [
    "Following the definition of the pure state, we are going to use the computational basis which is an orthonormal basis (i.e. {$a=\\ket0,\\ket1$}):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "2c50a6f7-b7e7-450d-b997-0778d1843ab3",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.30342183 0.95285633]\n"
     ]
    }
   ],
   "source": [
    "sqrt_eig_val = np.sqrt(rho_eig_val)\n",
    "print(sqrt_eig_val)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee62ac24-dccb-4ad5-9b53-6e5995e2f4ea",
   "metadata": {},
   "source": [
    "Our eigenvalues would be:\n",
    "\n",
    "\\begin{align}\n",
    "\\sqrt{\\lambda_1} = 0.303421, \\sqrt{\\lambda_2} = 0.952856\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3eb7a844-9699-44d2-b2d2-d00a60f0ac05",
   "metadata": {},
   "source": [
    "Using the previous definition, our pure state would be: \n",
    "\n",
    "\\begin{align}\n",
    "\\ket{\\psi} = \\sum_i \\sqrt{\\lambda_i} \\ket{v_i} \\otimes \\ket{a_i}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd6067ce-ac16-482b-89e3-d1fbedb00644",
   "metadata": {
    "tags": []
   },
   "source": [
    "\\begin{align}\n",
    "\\ket{\\psi} = \\sqrt{\\lambda_1} \\ket{v_1} \\otimes \\ket{0} + \\sqrt{\\lambda_2} \\ket{v_2} \\otimes \\ket{1}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1396b29c-19fd-4417-b626-fa874cd10887",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "\\ket{\\psi} = 0.303421\n",
    "\\begin{pmatrix}\n",
    "0.707106 \\\\\n",
    "-0.707106\n",
    "\\end{pmatrix} \\otimes \n",
    "\\begin{pmatrix}\n",
    "1 \\\\\n",
    "0\n",
    "\\end{pmatrix} + 0.952856\n",
    "\\begin{pmatrix}\n",
    "-0.707106 \\\\\n",
    "-0.707106\n",
    "\\end{pmatrix} \\otimes \n",
    "\\begin{pmatrix}\n",
    "0 \\\\\n",
    "1\n",
    "\\end{pmatrix}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bec2187e-f8e9-474b-a5ec-973ce6c00196",
   "metadata": {
    "tags": []
   },
   "source": [
    "\\begin{align}\n",
    "\\ket{\\psi} = 0.303421\n",
    "\\begin{pmatrix}\n",
    "0.707106 \\\\\n",
    "0 \\\\\n",
    "-0.707106  \\\\\n",
    "0\n",
    "\\end{pmatrix} + 0.952856\n",
    "\\begin{pmatrix}\n",
    "0 \\\\\n",
    "-0.707106 \\\\\n",
    "0  \\\\\n",
    "−0.707106\n",
    "\\end{pmatrix}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eaf4515c-733b-40c2-8ee0-e12dc23b898b",
   "metadata": {},
   "source": [
    "We must use the *flip* function of Numpy because the eigenvectors provided by the linalg library are actually the columns and not the rows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "c43fb834-0d9d-41c9-a024-37b7a395784b",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.70710678  0.          0.70710678  0.        ]\n",
      "[ 0.         -0.70710678  0.         -0.70710678]\n"
     ]
    }
   ],
   "source": [
    "tensor_product1 = np.vstack((np.flip(rho_eig_vec[1]),np.zeros(2))).ravel('F')\n",
    "tensor_product2 = np.vstack((np.zeros(2),np.flip(rho_eig_vec[0]))).ravel('F')\n",
    "print(tensor_product1)\n",
    "print(tensor_product2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "8fe46133-e7fc-49bc-9774-47966966484e",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.21455164 -0.67377117  0.21455164 -0.67377117]\n"
     ]
    }
   ],
   "source": [
    "psi = sqrt_eig_val[0]*tensor_product1 + sqrt_eig_val[1]*tensor_product2\n",
    "print(psi)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6cfbe95-8483-4318-9297-ad9080f40120",
   "metadata": {},
   "source": [
    "Our purified state would be: \n",
    "\n",
    "\\begin{align}\n",
    "\\ket{\\psi} = \n",
    "\\begin{pmatrix}\n",
    "-0.214551 & -0.673771 \\\\\n",
    "0.214551 & -0.673771 \\\\\n",
    "\\end{pmatrix}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc2e9659-35ca-4544-9871-79bb2d860467",
   "metadata": {},
   "source": [
    "The way to confirm the state above is purified is checking if we can recover the original density matrix $\\rho$. \n",
    "\n",
    "\\begin{align}\n",
    "\\rho = Tr_a(\\ket{\\psi}\\bra{\\psi}) \n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "d43c48b7-9995-483f-9c20-b4486e149f7a",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.0460324   0.14455871 -0.0460324   0.14455871]\n",
      " [ 0.14455871  0.4539676  -0.14455871  0.4539676 ]\n",
      " [-0.0460324  -0.14455871  0.0460324  -0.14455871]\n",
      " [ 0.14455871  0.4539676  -0.14455871  0.4539676 ]]\n"
     ]
    }
   ],
   "source": [
    "rho_partial_trace = np.dot(psi.reshape((4,1)),psi.reshape((4,1)).transpose())\n",
    "print(rho_partial_trace)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bf8dc6d-1385-449a-a5d4-57c3bf9379c2",
   "metadata": {},
   "source": [
    "As the reader can verify, taking the partial trace of the previous matrix lead us to the original mixed state with density matrix $\\rho$ which proves that $\\ket{\\psi}$ is a purified state. \n",
    "\n",
    "\\begin{align}\n",
    "Tr_a(\\ket{\\psi}\\bra{\\psi}) = \\begin{pmatrix}\n",
    "0.5 & 0.407935 \\\\\n",
    "0.407935 & 0.5\n",
    "\\end{pmatrix}\n",
    "= \\rho\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e967fbf7-1bf1-4002-8203-c90635000d49",
   "metadata": {},
   "source": [
    "## 3. Purity calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae992d3a-7c79-4e2a-9f70-0f55b3ea6ceb",
   "metadata": {},
   "source": [
    "Before getting into the details of calculating the purity, we must discuss what the purity is and how to find it in our scenario. The purity is a measure which states how much a quantum state is mixed and is defined by: \n",
    "\n",
    "\\begin{align}\n",
    "P = Tr(\\rho^2) = \\sum_{i} p_i^2\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdde1f56-cd3b-46e8-86c6-34820185f79d",
   "metadata": {},
   "source": [
    "Using the useful definitions and theory provided by [5], we can see that the purity can be expressed in terms of the lenght $r = || \\vec{r} ||$ of the **Bloch vector**: \n",
    "\n",
    "\\begin{align}\n",
    "P = Tr(\\rho^2) = \\frac{1 + r^2}{2}\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "2P - 1 = r^2\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "\\sqrt{2P - 1} = || \\vec{r} ||\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "319c59e8-0f7d-4979-be0b-f6cee64a269d",
   "metadata": {},
   "source": [
    "Using the definition of the eigenvalues of the density matrix $\\rho$ expressed in terms of the lenght $r$ of the Bloch vector:\n",
    "\n",
    "\\begin{align}\n",
    "\\lambda_{1,2} =  \\frac{(1 \\pm \\sqrt{(x^2 + y^2 + z^2)})}{2} = \\frac{(1 \\pm ||\\vec{r} ||)}{2} = \\frac{(1 \\pm \\sqrt{2P-1})}{2}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af93da6a-6d6c-49c0-b2dc-c54a7653716a",
   "metadata": {},
   "source": [
    "The previous equations states that we can calculate the eigenvalues of our covariance matrix (i.e. the eigenvalues for the Principal Components) by measuring the purity of our quantum state. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32d7dc6d-783a-43d7-a05e-12432e3aa6b3",
   "metadata": {},
   "source": [
    "Now, following the same procedure presented by [6], the purity is equal to the expected value of the SWAP gate between two purification copies (using also the previous definitions):"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6ff8507-f63f-4ce3-9870-19025ee9f249",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "\\rho = \\sum_{i} p_i \\ket{v_i} \\bra{v_i}\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "\\ket{\\psi} = \\sum_{i} \\sqrt{p_i} \\ket{v_i} \\otimes \\ket{w_i} = \\sum_{i} \\sqrt{p_i} \\ket{v_i} \\ket{w_i}\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "P = Tr(\\rho^2) = \\sum_{i} p_i^2\n",
    "\\end{align}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "099b3159-f048-4c5c-8729-51a752213986",
   "metadata": {},
   "source": [
    "Now, let's see what the SWAP operation makes to the purified state:\n",
    "\n",
    "\\begin{align}\n",
    "\\bra{\\psi} \\bra{\\psi} SWAP_{a} \\ket{\\psi} \\ket{\\psi} = \n",
    "\\sum_{ij} \\bra{w_{j}} \\bra{v_{j}} \\bra{w_{i}}  \\bra{v_{i}} \\sqrt{p_i} \\sqrt{p_j} \\ket{v_{j}} \\ket{w_{i}} \\ket{v_{i}} \\ket{w_{j}} = \n",
    "\\sum_{i} p_i^2\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30fd3ea9-f285-401b-b8b3-bf583cbc2c66",
   "metadata": {},
   "source": [
    "Then, \n",
    "\n",
    "\\begin{align}\n",
    "P = Tr(\\rho^2) = \\bra{\\psi} \\bra{\\psi} SWAP_{a} \\ket{\\psi} \\ket{\\psi}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03e6cd75-b1c6-404e-af5f-fee77fb65ff4",
   "metadata": {},
   "source": [
    "The previous equation states that purity can be measured using the Hadamard test on a controlled-SWAP gate. We're going to use the Amazon Braket service to find the purity of our quantum system to finally convert that data into the eigenvalues we're looking for our PCA scenario.\n",
    "\n",
    "To do this, we're going to follow the algorithm proposed by [4] which involves two processes: the state preparation and the purity calculation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "1a640cd5-2dd9-4a57-a6ae-0ec4e4a8755b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# AWS imports: Import Braket SDK modules\n",
    "from braket.circuits import Circuit\n",
    "from braket.devices import LocalSimulator\n",
    "from braket.circuits import Gate\n",
    "from braket.circuits import Observable\n",
    "import matplotlib.pyplot as plt\n",
    "from braket.aws import AwsDevice"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c7acacf-df03-4090-87bb-2da52c2e346e",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Implementation on Amazon Braket"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "161d7903-b604-4235-8ffb-fec5f7fefb18",
   "metadata": {},
   "source": [
    "We just need to create our quantum circuit with 5 qubits (two copies for the quantum state which means we need 4 qubits and a fifth qubit to interact with the answer which is known as an ancilla qubit). Each pair of qubits must be initialized to start with the quantum state that we already found on the previous part of this notebook (i.e. $\\ket{\\psi}$ and the last qubit just needs to be initialized as the $\\ket{0}$ state. Then, we need to apply a Hadamard gate to the ancilla qubit so we can use the superposition property before applying the controlled-SWAP gate between the ancilla qubit and the pair of quantum states (here the basis state is the target qubit and the control qubits are the two copies of the quantum state). Finally, another Hadamard gate to recover the answer we're looking for and the measure over the ancilla qubit. \n",
    "\n",
    "Then, we need to implement each step using the Amazon Braket SDK because not all quantum providers offer the same gates, that's why we need to decompose complex gates (e.g. Controlled-SWAP gate) into simple gates that any quantum hardware provider can run. That is the reason why we are using the Amazon Braket SDK since it allows us to write one circuit that can be compiled down to run on various devices."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c84310a-1f2f-45c9-9220-0e0977094411",
   "metadata": {},
   "source": [
    "Let's start by defining the quantum gate that we'll use many times during our quantum circuit: the unitary gate. This unitary gate is a generic single-qubit rotation gate with 3 Euler angles defined by:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d407730b-90cd-48eb-a94b-4cb7fe451a89",
   "metadata": {
    "tags": []
   },
   "source": [
    "\\begin{align}\n",
    "U(\\theta, \\phi, \\lambda) = \n",
    "\\begin{pmatrix}\n",
    "cos(\\frac{\\theta}{2}) & -e^{i\\lambda}sin(\\frac{\\theta}{2}) \\\\\n",
    "e^{i\\phi}sin(\\frac{\\theta}{2}) & e^{i(\\phi+\\lambda)}cos(\\frac{\\theta}{2}) \\\\\n",
    "\\end{pmatrix}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "3c00b87b-abae-43ab-b483-c74a96c21c48",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def rotation_matrix(value):\n",
    "    return np.array([[np.cos(value/2), -np.exp(1j*value)*np.sin(value/2)],[np.exp(1j*value)*np.sin(value/2), np.exp(2*1j*value)*np.cos(value/2)]])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e8b1028-9572-4fdd-90c2-e2082368b9c6",
   "metadata": {},
   "source": [
    "Considering our special case has only 2 features, for the PCA process we are doing the following values can be calculated using the above equations and data. First, we need to take the basis states to the quantum state $\\ket{\\psi}$ we found earlier (this is the state preparation process using unitary matrices, hadamard and cnot gates). Then, for the purity calculation process we need to decompose the controlled-swap operation into simple gates (the way of doing this is beyond the scope of this notebook, but if the reader wants further details about how to do decomposition using the Toffoli gate please see the references provided by [4] on this topic). Finally, we'd use the same strategy of the ancilla qubit to measure the quantum circuit to get the eigenvalues we're looking for."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "49f2f954-03af-4f6e-bc08-059f52d43edc",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T  : |0|1|2|3|4|5|6|7 |8|9|10|11|12|13|14|15|16|17|18|19|20|21|22|23|            Result Types            |\n",
      "                                                                                                          \n",
      "q0 : ---X-U-----------------------------------------------------------------------------------------------\n",
      "        |                                                                                                 \n",
      "q1 : -U-C-U-H-X---X-Ti-X---X--T--X--Ti-X-----X--T-----------------X---------------------------------------\n",
      "              |   |    |   |     |     |     |                    |                                       \n",
      "q2 : -U-C-U-H-C-H-C----C-X-C--X--C-----C--X--C--X--T--X--Ti-X--H--C---------------------------------------\n",
      "        |                |    |           |     |     |     |                                             \n",
      "q3 : -H-|----------------C----C-----------C-----C-----C--T--C--H-----Expectation(Z)-Sample(Z)-Probability-\n",
      "        |                                                                                                 \n",
      "q4 : ---X-U-----------------------------------------------------------------------------------------------\n",
      "\n",
      "T  : |0|1|2|3|4|5|6|7 |8|9|10|11|12|13|14|15|16|17|18|19|20|21|22|23|            Result Types            |\n"
     ]
    }
   ],
   "source": [
    "circuit = Circuit()\n",
    "\n",
    "for i in [1,2]:\n",
    "    circuit.unitary(matrix=rotation_matrix(0.465), targets=[i])\n",
    "\n",
    "circuit.h(3)\n",
    "\n",
    "for (i, j) in [[1,0], [2,4]]:\n",
    "    circuit.cnot(i, j)\n",
    "    \n",
    "for i in [0,4]:\n",
    "    circuit.unitary(matrix=rotation_matrix(1.570), targets=[i])\n",
    "\n",
    "for i in [1,2]:\n",
    "    circuit.unitary(matrix=rotation_matrix(1.950), targets=[i])\n",
    "\n",
    "for i in [1,2]:\n",
    "    circuit.h(i)\n",
    "\n",
    "circuit.cnot(2,1)\n",
    "\n",
    "circuit.h(2)\n",
    "\n",
    "circuit.cnot(2,1)\n",
    "\n",
    "circuit.ti(1)\n",
    "\n",
    "for (i, j) in [[2,1], [3,2], [2,1]]:\n",
    "    circuit.cnot(i, j)\n",
    "\n",
    "circuit.t(1)\n",
    "\n",
    "for (i, j) in [[3,2], [2,1]]:\n",
    "    circuit.cnot(i, j)\n",
    "\n",
    "circuit.ti(1)\n",
    "\n",
    "for (i, j) in [[2,1], [3,2], [2,1]]:\n",
    "    circuit.cnot(i, j)\n",
    "\n",
    "for (i,j,k) in [[1,3,2],[2,3,2]]:\n",
    "    circuit.t(i)\n",
    "    circuit.cnot(j,k)\n",
    "\n",
    "circuit.ti(2)\n",
    "circuit.t(3)\n",
    "\n",
    "circuit.cnot(3,2)\n",
    "\n",
    "for i in [2,3]:\n",
    "    circuit.h(i)\n",
    "\n",
    "circuit.cnot(2,1)\n",
    "\n",
    "#measurement part of the quantum circuit\n",
    "circuit.expectation(Observable.Z(), target=[3])\n",
    "circuit.sample(observable=Observable.Z(), target=3)\n",
    "circuit.probability(target=3)\n",
    "\n",
    "print(circuit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f4b0678-6825-4468-a470-d0aa78f23050",
   "metadata": {},
   "source": [
    "## Amazon Braket local simulator"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1969209f-13aa-4d06-992e-8d6fff48a123",
   "metadata": {},
   "source": [
    "For the full article, please go to the following [link](https://aws.amazon.com/es/blogs/quantum-computing/simulating-quantum-circuits-with-amazon-braket/).\n",
    "\n",
    "Today’s quantum computers still have limited capacity, qubit counts, and accuracy. Therefore, quantum circuit simulators, which emulate the behavior of quantum computers using classical hardware, are often the tool of choice to study quantum algorithms for quantum computing researchers and enthusiasts. Amazon Braket provides a suite of quantum simulators for different use cases to help customers prototype and develop quantum algorithms [8].\n",
    "\n",
    "The most common use case for simulators is to prototype, validate, and debug quantum programs. Whether you want to understand the foundation of Simon’s algorithm, or validate that your variational algorithm runs without a bug before running it on a quantum computer, it is useful to have a quick way to run small-scale circuits. For this reason, the Amazon Braket SDK comes pre-installed with a local simulator that runs wherever you use the SDK, for example, on Amazon Braket notebooks, or your own laptop [8]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "922e6d47-5bd4-49d3-8bf1-52da18f7f027",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Counter({'01100': 1585, '11100': 1039, '01101': 1009, '11101': 683, '10001': 597, '10101': 589, '11001': 588, '11000': 511, '01001': 493, '00101': 491, '10100': 479, '01000': 419, '00100': 415, '10000': 263, '00001': 226, '10110': 139, '00111': 122, '01011': 118, '11010': 118, '00000': 116})\n"
     ]
    }
   ],
   "source": [
    "local_device = LocalSimulator()\n",
    "\n",
    "local_result = local_device.run(circuit, shots=10000).result()\n",
    "local_counts = local_result.measurement_counts\n",
    "print(local_counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "b2776f8c-a66d-4407-91b9-e31056e8231c",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Probability: [0.9503 0.0497]\n"
     ]
    }
   ],
   "source": [
    "print(\"Probability:\", local_result.values[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "bbe77e55-f74e-4678-afcc-b1a9a3626ad8",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9006000000000001\n"
     ]
    }
   ],
   "source": [
    "purity_sdk = local_result.values[2][0] - local_result.values[2][1]\n",
    "print(purity_sdk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "d234b2f3-3831-4809-9ca1-d17c554b057c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The first eigenvalue obtained by the quantum PCA using Amazon Braket SDK is: \n",
      " 2.03046188574001\n"
     ]
    }
   ],
   "source": [
    "v_1 = (1 + np.sqrt(2 * purity_sdk - 1)) / 2 * np.trace(sigma)\n",
    "print('The first eigenvalue obtained by the quantum PCA using Amazon Braket SDK is: \\n', v_1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3d4480f-483e-4c31-9456-1e3e19c0ae7c",
   "metadata": {},
   "source": [
    "Recall the eigenvalues we found using linear algebra at the beginning of this notebook were:\n",
    "\n",
    "\\begin{align}\n",
    "e_1 = 1.945575, e_2 = 0.197281\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "773a359a-e186-4b47-a155-c43edfd0685d",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "percent error first eigenvalue (%):  4.363074450484305\n"
     ]
    }
   ],
   "source": [
    "perc_v1 = abs((v_1-1.945575)/1.945575)*100\n",
    "print('percent error first eigenvalue (%): ', perc_v1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8b015aa-5725-4454-bb9d-bcfe1e5292a3",
   "metadata": {},
   "source": [
    "## Amazon Braket SV1 device"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea352ce9-ff56-4c4b-9fb4-f5bd7c88eb38",
   "metadata": {},
   "source": [
    "State vector simulators are the workhorse of circuit simulation. They keep a precise representation of the quantum state (the state vector) at every point in the simulation, and iteratively change the state under the action of the different gates of the circuit, one by one. Each gate that is applied corresponds to a matrix-vector multiplication resulting in predictable runtimes of SV1 with little variation [8]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "9d123c53-bb72-4134-8acf-b7e5f3f5e0f4",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Counter({'01100': 1562, '11100': 1044, '01101': 1021, '10101': 641, '11001': 630, '11101': 619, '10001': 608, '00101': 497, '10100': 489, '01001': 489, '11000': 450, '00100': 445, '01000': 426, '00001': 257, '10000': 249, '01011': 133, '10110': 129, '11010': 121, '00000': 99, '00111': 91})\n"
     ]
    }
   ],
   "source": [
    "sv1_device = AwsDevice(\"arn:aws:braket:::device/quantum-simulator/amazon/sv1\")\n",
    "\n",
    "sv1_result = sv1_device.run(circuit, shots=10000).result()\n",
    "sv1_counts = sv1_result.measurement_counts\n",
    "print(sv1_counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "d739d9b7-b3ba-479a-829b-be69b2216ec5",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Expectation value for: [0.9526 0.0474]\n"
     ]
    }
   ],
   "source": [
    "print(\"Expectation value for:\", sv1_result.values[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "fc79fad8-d2f4-4838-93cb-9e2dc4348cc2",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9052\n"
     ]
    }
   ],
   "source": [
    "purity_sv1 = sv1_result.values[2][0] - sv1_result.values[2][1]\n",
    "print(purity_sv1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "5e5608e5-981c-42af-abfd-3774fc09e5df",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The first eigenvalue obtained by the quantum PCA using Amazon Braket SV1 is: \n",
      " 2.0359523515651623\n"
     ]
    }
   ],
   "source": [
    "a_1 = (1 + np.sqrt(2 * purity_sv1 - 1)) / 2 * np.trace(sigma)\n",
    "print('The first eigenvalue obtained by the quantum PCA using Amazon Braket SV1 is: \\n', a_1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af2efca9-c263-4519-83d5-fc8eb29eaee1",
   "metadata": {},
   "source": [
    "Recall the eigenvalues we found using linear algebra at the beginning of this notebook were:\n",
    "\n",
    "\\begin{align}\n",
    "e_1 = 1.945575, e_2 = 0.197281\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "4f6a91eb-0592-4889-8c47-37f87a7a1a36",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "percent error first eigenvalue (%):  4.6452771836172975\n"
     ]
    }
   ],
   "source": [
    "perc_sv1 = abs((a_1-1.945575)/1.945575)*100\n",
    "print('percent error first eigenvalue (%): ', perc_sv1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6eab3eb-3985-4f8a-80fa-5bc021ecd444",
   "metadata": {},
   "source": [
    "# Conclusions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d6e9c55-28ff-4f12-bf05-5f09948e5d02",
   "metadata": {},
   "source": [
    "Summarizing what we've done, along the entire notebook we implemented the Principal Component Analysis technique under the quantum computation paradigm. We implemented the PCA technique for the famous problem of house pricing in the US where our final goal is to find the eigenvalues of our correlation matrix of our data (i.e. the principal components). The way how we implemented the PCA technique was divided in four stages: i) Classical pre-processing, ii) State preparation, iii) Purity calculation, iv) Classical post-processing. After explaining each step with the mathematics and physics details needed, we implemented the quantum algorithm using the Amazon Braket SDK with two devices: the local simulator and the Amazon Braket SV1. \n",
    "\n",
    "Our results for the eigenvalue we were looking for is the same we predicted using classical computation with a percent error less than 5%. Finally, we saw that for this particular scenario, PCA technique was correctly implemented using quantum algorithms, circuits and devices. However, we know the technology is still in an early stage, but proving that we can make these small steps brings hope and optimism for the future of quantum computing and its applications on machine learning. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb255b4e-9dc9-422b-9519-5299fe077c02",
   "metadata": {
    "tags": []
   },
   "source": [
    "# References\n",
    "\n",
    "1. Ian T. Jolliffe and Jorge Cadima. Principal component analysis: A review and recent developments. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 374, April 2016.\n",
    "\n",
    "2. C. He, J. Li, W. Liu, J. Peng, and Z. J. Wang. A low-complexity quantum principal component analysis algorithm,. IEEE Transactions on Quantum Engineering, 3(3):1–13, 2022.\n",
    "\n",
    "3. Lin Jie et al. An improved quantum principal component analysis algorithm based on the quantum singular threshold method. Physics Letters A, 383:2862–68, aug 2019.\n",
    "\n",
    "4. J. , Abhijith, et al. «Quantum Algorithm Implementations for Beginners». ACM Transactions on Quantum Computing, vol. 3, n.o 4, dic 2022, pp. 1-92. arXiv.org, https://doi.org/10.1145/3517340.\n",
    "\n",
    "5. Schmied, Roman. «Quantum State Tomography of a Single Qubit: Comparison of Methods». Journal of Modern Optics, vol. 63, n.o 18, oct 2016, pp. 1744-58. DOI.org (Crossref), https://doi.org/10.1080/09500340.2016.1142018.\n",
    "\n",
    "6. An special thanks to the Github repository of Haokai Zhang available at: https://github.com/Haokai-Zhang/ExampleQPCA/blob/master/5qubit-qPCA.ipynb\n",
    "\n",
    "7. Kelley Pace, R., y Ronald Barry. «Sparse spatial autoregressions». Statistics & Probability Letters, vol. 33, n.o 3, may 1997, pp. 291-97. ScienceDirect, https://doi.org/10.1016/S0167-7152(96)00140-X.\n",
    "\n",
    "8. Simulating Quantum Circuits with Amazon Braket | AWS Quantum Technologies Blog. sept 2021, https://aws.amazon.com/blogs/quantum-computing/simulating-quantum-circuits-with-amazon-braket/.\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "conda_braket",
   "language": "python",
   "name": "conda_braket"
  },
  "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}