12-towards measurement mitigation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d2d70b4f-799e-4ae0-96b6-ecb3132c6b67",
   "metadata": {},
   "source": [
    "# Towards Quantum Measurement Error Mitigation\n",
    "\n",
    "## How to connect CDR with quantum state tomography\n",
    "\n",
    "Quantum computers are astonishing devices. Yet, they are error-prone, too. Therefore, we need to implement quantum error mitigation methods to reduce the negative effect of errors on our computation results. \n",
    "\n",
    "In a series of previous posts, we [learned the Clifford Data Regression method](https://pyqml.medium.com/mitigating-quantum-errors-using-clifford-data-regression-98ab663bf4c6) and mitigated errors in a [simulated environment](https://towardsdatascience.com/how-to-implement-quantum-error-mitigation-with-qiskit-and-mitiq-e2f6a933619c) and on a [real quantum computer](https://towardsdatascience.com/practical-error-mitigation-on-a-real-quantum-computer-41a99dddf740).\n",
    "\n",
    "The results are encouraging. Yet, an unexpected obstacle appeared when I tried to use it to participate in IBM's Quantum Open Science Prize.\n",
    "\n",
    "IBM asks us to simulate a three-particle Heisenberg Hamiltonian using Trotterization. No, that's not the problem. The problem is that they assess any submission through quantum state tomography. This is an approach to recreate a quantum state through measurements. More specifically, the problem is that they use the `StateTomographyFitter` inside Qiksit. This implementation builds upon experiment counts. But the CDR method works with expectation values.\n",
    "\n",
    "Let me illustrate the problem a little bit. The following figure depicts the simple case of a 1-qubit quantum circuit.\n",
    "\n",
    "![](assets/ev.png)\n",
    "\n",
    "Whenever we look at a qubit, it is either 0 or 1. That's it. Which of both depends on chance and the internal quantum state. Say the qubit is in state $|+\\rangle$. In this state, the measuring 0 and 1 have the same probability. In this state, measuring zero and one have the same probability. But, we can't conclude on the probabilities when running the quantum circuit only once.\n",
    "\n",
    "But when we run it repeatedly, say a thousand times, we would see zero and one 500 times each, except for a slight statistical variance.\n",
    "\n",
    "Unfortunately, current quantum computers are noisy devices. We sometimes measure a 0 for a 1 and vice versa. So, the results become blurred. For instance, we would measure 0 412 times. \n",
    "\n",
    "The raw result of running a quantum circuit is the measurement. Since we almost always execute circuits multiple times, we count the measurements. \n",
    "\n",
    "So, let's look at the source code in Qiskit of such a circuit.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "681eac5c-e565-4543-8af6-a7b1f16cda03",
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit import QuantumCircuit, execute, Aer\n",
    "\n",
    "# Create a quantum circuit with one qubit\n",
    "qc = QuantumCircuit(1)  \n",
    "qc.h(0)\n",
    "qc.measure_all()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b5dd85f-d9eb-4214-965b-26ec3887fc1c",
   "metadata": {},
   "source": [
    "We define a quantum circuit with a single qubit and apply the Hadamard gate to it. This transformation gate puts the qubit (from the initial state $|0\\rangle$) into state $|+\\rangle$. At the end of the circuit, we measure the qubit. The following image shows the circuit diagram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f2cf0c98-1d3a-4e0d-b352-4dde3a111eab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"word-wrap: normal;white-space: pre;background: #fff0;line-height: 1.1;font-family: &quot;Courier New&quot;,Courier,monospace\">        ┌───┐ ░ ┌─┐\n",
       "     q: ┤ H ├─░─┤M├\n",
       "        └───┘ ░ └╥┘\n",
       "meas: 1/═════════╩═\n",
       "                 0 </pre>"
      ],
      "text/plain": [
       "        ┌───┐ ░ ┌─┐\n",
       "     q: ┤ H ├─░─┤M├\n",
       "        └───┘ ░ └╥┘\n",
       "meas: 1/═════════╩═\n",
       "                 0 "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qc.draw(output='text')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11c0892d-d63e-4a1d-97d7-9264c4d7b650",
   "metadata": {},
   "source": [
    "Let's run this circuit 1,000 times. We do this by getting a `backend` (here, a noise-free statistical simulator) and calling the `execute` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "16d1ce16-9ab3-4596-9251-8981b6f4e84e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Tell Qiskit how to simulate our circuit\n",
    "backend = Aer.get_backend('qasm_simulator') \n",
    "\n",
    "# Do the simulation, returning the result\n",
    "result = execute(qc, backend, shots=1000).result()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dd15531-e0f4-48a7-a0ee-9cc2908f9b43",
   "metadata": {},
   "source": [
    "We get back a Qiskit result object ([Qiskit reference](https://qiskit.org/documentation/stubs/qiskit.result.Result.html)). The essential function of this object is `get_counts` because it tells us what we saw when looking at the qubit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c4e8e689-7fee-45b7-bcfd-2752b0b2142c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'1': 499, '0': 501}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result.get_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d078d0d6-16b1-41a5-b564-422ab3df62a6",
   "metadata": {},
   "source": [
    "It is a simple Python dictionary with the measurement result as the key and the number of times we observed this result as the value. It is now up to us to interpret these counts and do something meaningful. From now on, these results are as good as any other statistical data. We can use them to calculate further values, such as the expectation value. This is the probabilistic expected value of the measurement of an experiment. It is similar to a classical expectation value. For example, consider tossing a fair coin that lands on heads and tails equally likely. If you assign the value 1 to heads and 0 to tails, the expectation value is $0.5*1+0.5*0=0.5$. It is the same for a qubit in state $|+\\rangle$.\n",
    "\n",
    "Usually, when you see the calculation of expectation values in quantum computing, it looks like this $\\langle \\psi|Z|\\psi\\rangle$. The letter \"psi\" ($\\psi$) denotes the quantum state, and the Z in the middle symbolizes the observable. In this case, it is the Z-axis. \n",
    "\n",
    "The important thing here is that this notation introduces the concept of an observable. When we measured our qubit in Qiskit before, we implicitly chose the Z-observable because this is the default measurement basis for a qubit in Qiskit. So, we are essentially talking about the same concept. In [this post](https://towardsdatascience.com/how-to-implement-quantum-error-mitigation-with-qiskit-and-mitiq-e2f6a933619c), we looked at observables in more detail. The one thing that we need to know is that there is not a single observable, but many. Just think about a sphere like the earth. The observable is the specific point from which you look at the globe. The world looks different depending on the perspective, yet, it is the same.\n",
    "\n",
    "![](assets/earth.png)\n",
    "\n",
    "Essentially, the counts and the expectation value are closely tight together. They both have their uses. While the counts contain more information about the different measurements, the expectation value is better to work with because it is a single number.\n",
    "\n",
    "This is the point of struggle I described in my previous post. While the quantum error mitigation method CDR uses the simplicity of the expectation value, the quantum state tomography that IBM uses to evaluate the performance of the error mitigation works with the counts.\n",
    "\n",
    "To participate in IBM's challenge, it is now our job to integrate the two.\n",
    "\n",
    "I chose to adapt the CDR to work with an unchanged state tomography because the latter is IBM's evaluation method. I believe messing around with their assessment tool may disqualify myself from the contest.\n",
    "\n",
    "So, we need to change the CDR method to change a counts dictionary instead of a single number.\n",
    "\n",
    "Let's revisit the CDR briefly.\n",
    "\n",
    "![](assets/cdr_steps.png)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf7c837c-262f-4b10-b7e3-5894ecfde15a",
   "metadata": {},
   "source": [
    "The CDR method has three steps. First, we generate training data. Then, in this step, we run our quantum circuit twice. Once on a classical computer to obtain the exact value of an observable's expectation value of interest. And once on a real quantum computer to produce the noisy value.\n",
    "\n",
    "In the second step, we create a linear model of the relationship between the noisy and the exact values. Building a linear model from a set of data points is known as regression. Therefore the name CDR--Clifford Data Regression.\n",
    "\n",
    "Finally, we use the model to transform a noisy expectation value into a mitigated one by predicting the noise-free value."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "147c5548-216e-41fe-ae65-cd29b909a830",
   "metadata": {},
   "source": [
    "All these steps need to work with the Qiskit [experiment result](https://qiskit.org/documentation/stubs/qiskit.result.Result.html). However, the problem is that this is an object that the Qiskit `execute` function creates. It stores most of its data in a read-only manner that we can't change any more.\n",
    "\n",
    "But, we can apply a little trick. We write our own `Result` class that allows us to change the counts afterward.\n",
    "\n",
    "Conceptually, we create a new Python class that serves the one function the state tomography uses. This is the `get_counts` function. So, when the state tomography function queries the counts, it gets a response. But since we implement this new class, we can also provide a function that overwrites the counts.\n",
    "\n",
    "The following listing depicts the source code of our `OwnResult` class.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "c7025efd-30f7-4a69-831a-fdc27cb26ca2",
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.result import Result\n",
    "\n",
    "from qiskit.circuit.quantumcircuit import QuantumCircuit\n",
    "from qiskit.exceptions import QiskitError\n",
    "from qiskit.result.counts import Counts\n",
    "\n",
    "class OwnResult(Result):\n",
    "    \n",
    "    def __init__(self, result):\n",
    "        self._result = result\n",
    "        self._counts = {}\n",
    "            \n",
    "        \n",
    "    def get_counts(self, experiment=None):\n",
    "\n",
    "        if experiment is None:\n",
    "            exp_keys = range(len(self._result.results))\n",
    "        else:\n",
    "            exp_keys = [experiment]\n",
    "        \n",
    "\n",
    "        dict_list = []\n",
    "        for key in exp_keys:\n",
    "            exp = self._result._get_experiment(key)\n",
    "            try:\n",
    "                header = exp.header.to_dict()\n",
    "            except (AttributeError, QiskitError):  # header is not available\n",
    "                header = None\n",
    "\n",
    "            if \"counts\" in self._result.data(key).keys():\n",
    "                if header:\n",
    "                    counts_header = {\n",
    "                        k: v\n",
    "                        for k, v in header.items()\n",
    "                        if k in {\"time_taken\", \"creg_sizes\", \"memory_slots\"}\n",
    "                    }\n",
    "                else:\n",
    "                    counts_header = {}\n",
    "                    \n",
    "                    \n",
    "                # CUSTOM CODE STARTS HERE #######################\n",
    "                dict_list.append(Counts(\n",
    "                    self._counts[str(key)] if str(key) in map(lambda k: str(k), self._counts.keys()) else self._result.data(key)[\"counts\"]\n",
    "                    , **counts_header))\n",
    "                # CUSTOM CODE ENDS HERE #######################\n",
    "                \n",
    "            elif \"statevector\" in self._result.data(key).keys():\n",
    "                vec = postprocess.format_statevector(self._result.data(key)[\"statevector\"])\n",
    "                dict_list.append(statevector.Statevector(vec).probabilities_dict(decimals=15))\n",
    "            else:\n",
    "                raise QiskitError(f'No counts for experiment \"{repr(key)}\"')\n",
    "\n",
    "        # Return first item of dict_list if size is 1\n",
    "        if len(dict_list) == 1:\n",
    "            return dict_list[0]\n",
    "        else:\n",
    "            return dict_list\n",
    "        \n",
    "        \n",
    "    def set_counts(self, counts, experiment=None):\n",
    "        self._counts[str(experiment) if experiment is not None else \"0\"] = counts"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4e498c2-3b76-4c8d-9758-d8694fa3eca2",
   "metadata": {},
   "source": [
    "The class takes the existing Qiskit result as an initialization parameter. Further, we specify `_counts` as a member variable that we initialize with an empty dictionary. This will hold our changed counts. The code of the `get_count` function is copied from the original source code except for two little thing. First, whenever we refer to the a property of the result, such as `data` we need to look into `self._result.data` instead of `self.data`. Second, at lines 40-44, we look into the custom member function for the actual counts. If they exist (`if str(key) in self._counts.keys()if str(key) in self._counts.keys()`) we return the changed counts (`self._counts[str(key)]self._counts[str(key)]`). If they don't exist, we return the original counts (`self._result.data(key)[\"counts\"]`)\n",
    "\n",
    "Let's have a look at how that works."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "400c9784-c9c9-4b87-b3b7-85c98cbc543e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "original result:  {'1': 499, '0': 501}\n",
      "Changed counts:   {'0': 100, '1': 900}\n"
     ]
    }
   ],
   "source": [
    "# create our own Result object\n",
    "own = OwnResult(result)\n",
    "print(\"original result: \", own.get_counts())\n",
    "\n",
    "# set new counts\n",
    "own.set_counts({\"0\": 100, \"1\": 900})\n",
    "print(\"Changed counts:  \", own.get_counts())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18986ad0-9dfb-4357-bde9-522043d6dfbc",
   "metadata": {},
   "source": [
    "So, let's see whether we can use our `OwnResult` inside the state tomography.\n",
    "\n",
    "In the following snippet, we create the same simple circuit we used before. The only difference is that we omit the measurement because we need to create the `state_tomography_circuits` from it. Then, we run these circuits on a noise-free `Qasm-Simulator` and store the result of it.\n",
    "\n",
    "So, now we're ready for the interesting part. We loop through the circuits in the list of experiments (`st_qcs`). For each of the circuits, we set the counts to a fixed dictionary with arbitrary values. We don't care about the values for now because we only want to verify whether the `StateTomographyFitter` works with our `OwnResult`.\n",
    "\n",
    "Finally, we compute the fidelity based on the original and the changed results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "e19ea5b5-20c6-4323-add8-7886525fba55",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "original fidelity: 0.49800015998080377\n",
      "changed fidelity:  0.7886751345948128\n"
     ]
    }
   ],
   "source": [
    "# Import state tomography modules\n",
    "from qiskit.ignis.verification.tomography import state_tomography_circuits, StateTomographyFitter\n",
    "from qiskit.quantum_info import state_fidelity\n",
    "from qiskit.opflow import Zero, One\n",
    "from qiskit.providers.aer import QasmSimulator\n",
    "from qiskit import execute, transpile, assemble\n",
    "\n",
    "# create a circuit without measurement\n",
    "qc = QuantumCircuit(1)  \n",
    "qc.h(0)\n",
    "\n",
    "# The expected final state; necessary to determine state tomography fidelity\n",
    "target_state = (One).to_matrix()  # |1>\n",
    "\n",
    "# create the three state tomography circuits\n",
    "st_qcs = state_tomography_circuits(qc, [0])\n",
    "\n",
    "# Noiseless simulated backend\n",
    "sim = QasmSimulator()\n",
    "\n",
    "job = execute(st_qcs, sim, shots=1000)\n",
    "\n",
    "# get the result\n",
    "result = job.result()\n",
    "\n",
    "# put the result in our own class\n",
    "own_res = OwnResult(result)\n",
    "\n",
    "# loop through the experiments\n",
    "for experiment in st_qcs:\n",
    "    exp_keys = [experiment]\n",
    "    \n",
    "    for key in exp_keys:\n",
    "        exp = result._get_experiment(key)        \n",
    "        own_res.set_counts({\"0\": 100, \"1\": 900}, key)\n",
    "\n",
    "\n",
    "# Compute fidelity with original result\n",
    "orig_tomo_fitter = StateTomographyFitter(result, st_qcs)\n",
    "orig_rho_fit = orig_tomo_fitter.fit(method='lstsq')\n",
    "orig_fid = state_fidelity(orig_rho_fit, target_state)\n",
    "print(\"original fidelity:\", orig_fid)\n",
    "\n",
    "\n",
    "# Compute fidelity with chnaged result\n",
    "tomo_fitter = StateTomographyFitter(own_res, st_qcs)\n",
    "rho_fit = tomo_fitter.fit(method='lstsq')\n",
    "fid = state_fidelity(rho_fit, target_state)\n",
    "print(\"changed fidelity: \", fid)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c57c765d-ed37-4ac1-b5da-0a9c527311d4",
   "metadata": {},
   "source": [
    "The output shows that the fidelity we calculated based on the changed counts differs significantly from the one we calculated based on the original counts. Apparently, the `StateTomographyFitter` works with our custom counts. This creates the prerequisite for mitigating errors in the next step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d4680ff-6a2c-4a52-b7f8-9bce4f454fe3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}