tutorials/quantum_simulation/QDRIFT_EN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hamiltonian Simulation with qDRIFT\n",
    "<em> Copyright (c) 2022 Institute for Quantum Computing, Baidu Inc. All Rights Reserved. </em>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "In quantum mechanics, the energy of the system is described by the Hamiltonian operator $H$, which determines the evolution of the system. So Hamiltonian simulation has great practical value in modeling complex chemical and physical systems. However, because the degree of freedom of the system increases exponentially with the increase of the system (such as the number of qubits), it is generally impossible to use classical computers to effectively simulate quantum systems. At present, the main technology of using quantum computers to simulate Hamiltonian is to use product formula method to simulate time evolution. This tutorial will introduce some basic theories and methods about product formula, and a random method named quantum stochastic drift protocol (qDRIFT) which is based on product formula. Then we give a code demonstration at the end of the article.\n",
    "\n",
    "\n",
    "## Product formula\n",
    "\n",
    "According to the basic axioms of quantum mechanics, the evolution of the system with Hamiltonian $H$ can be described by the following equation\n",
    "\n",
    "$$\n",
    "i \\hbar \\frac{d}{d t} | \\psi \\rangle = H | \\psi \\rangle,\n",
    "\\tag{1}\n",
    "$$\n",
    "\n",
    "where $\\hbar$ is the reduced Planck constant. Therefore, for a time-independent Hamiltonian, the time evolution equation of the system can be written as\n",
    "\n",
    "$$\n",
    "|\\psi(t) \\rangle = U(t) | \\psi (0) \\rangle, ~ U(t) = e^{- i H t}.\n",
    "\\tag{2}\n",
    "$$\n",
    "\n",
    "Here we take the natural unit $\\hbar=1$, $U (t) $ as the time evolution operator. The core idea of using quantum circuits to simulate the time evolution process is to use the unitary transformation constructed by quantum circuits to simulate and approximate the time evolution operator. Seth Lloyd pointed out in his 1996 article that a whole evolution time of $t$ can be divided into $r$ shorter \"time blocks\" to reduce the error of simulation in time evolution[1]. Consider a general Hamiltonian form $H = \\sum_ {k=1}^{L} H_k$, of which $H_k$ is sub-Hamiltonian acting on a part of the system. We consider each sub-Hamiltonian $H_k$ whose evolution operator is $e^ {-i H_k t}$, and we can get $\\prod_{k=1}^{L} e^{-i H_k t}$ by simulating each sub-Hamiltonian in turn. Through Taylor expansion, it can be found that\n",
    "\n",
    "$$\n",
    "e^{-iHt} = \\prod_{k=1}^{L} e^{-i H_k t} + O(t^2).\n",
    "\\tag{3}\n",
    "$$\n",
    "\n",
    "Then let $\\tau = t/r$ and consider the evolution operator $\\left (e^ {-iH \\tau}\\right) ^r$, we can deduce that\n",
    "\n",
    "$$\n",
    "e^{-iHt} = \\left(e^{-iH \\tau}\\right)^r = \\left(\\prod_{k=1}^{L} e^{-i H_k \\tau} + O(\\tau^2) \\right)^r = \\left(\\prod_{k=1}^{L} e^{-i H_k \\tau} \\right)^r + O\\left(\\frac{t^2}{r}\\right).\n",
    "\\tag{4}\n",
    "$$\n",
    "\n",
    "This formula tells us that as long as the whole evolution time can be divided into enough \"fragments\", we can simulate with any high simulation accuracy. That is the basic idea of the product formula. However, what is given in (4) is only a rough estimate. If we want to estimate the depth of the quantum circuits required to achieve a certain simulation accuracy, we need to calculate its rigorous error upper bound. Specifically, we make $U_ {circuit}$ represent the circuit we construct, $\\Vert \\cdot \\Vert$ is the Schatten-$\\infty$ norm, that is, the [spectral norm](https://en.wikipedia.org/wiki/Schatten_norm). Then the simulation error $\\epsilon$ can be written as\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\epsilon\\left(e^{-iH\\tau}, U_{circuit}\\right)  & = \\Vert e^{-iH\\tau} - U_{circuit}\\Vert .\n",
    "\\end{aligned}\n",
    "\\tag{5}\n",
    "$$\n",
    "\n",
    "Next, we will show a brief calculation process of the upper bound of error. We give two conclusions (6) and (7) without proof, which will be used in proving (8). Interested readers can refer to section F.1 in [2] for details.\n",
    "\n",
    "$$\n",
    "\\left\\Vert \\mathcal{R}_k \\left( \\prod_{k=1}^{L} e^{-i H_k \\tau} \\right) \\right\\Vert\n",
    "\\leq\n",
    "\\mathcal{R}_k \\left( e^{\\vert \\tau \\vert   \\sum_{k=1}^{L} \\Vert H_k \\Vert } \\right),\n",
    "\\tag{6}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\vert \\mathcal{R}_k(e^\\alpha) \\vert \\leq \\frac{\\vert \\alpha \\vert^{k+1}}{(k+1)!}  e^{ \\vert \\alpha \\vert }, ~\n",
    "\\forall \\alpha \\in \\mathbb{C},\n",
    "\\tag{7}\n",
    "$$\n",
    "\n",
    "where $\\mathcal{R}_ k (f) $ is the remainder Taylor expansion to order $k$ of the function $f$, such as $\\mathcal{R}_1 (e^x)=\\mathcal{R}_1 (\\sum_{j=0}^\\infty \\frac{x^n}{n!})=\\sum_{j=2}^\\infty \\frac{x^n}{n!}$. \n",
    "Make $\\Lambda = \\max_ k \\Vert H_ k \\Vert$, considering the complete evolution time $t = r \\cdot \\tau$, the error of simulation in complete time $t$ is:\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\left \\Vert \\left ( e^{-i\\tau \\sum_{k=1}^L H_k  }\\right)^r - \\left (\\prod_{k=1}^{L} e^{-i H_k \\tau} \\right)^r \\right \\Vert \\leq &\n",
    "r \\left \\Vert e^{-i\\tau \\sum_{k=1}^L H_k } - \\prod_{k=1}^{L} e^{-i H_k \\tau } \\right \\Vert  \\\\\n",
    "=& r \\left \\Vert \\mathcal{R}_1 \\left(  e^{-i\\tau \\sum_{k=1}^L H_k} \\right)- \\mathcal{R}_1 \\left( \\prod_{k=1}^{L} e^{-i H_k \\tau } \\right) \\right \\Vert \\\\\n",
    "\\leq& r \\left \\Vert \\mathcal{R}_1 \\left(  e^{-i\\tau \\sum_{k=1}^L H_k} \\right) \\right \\Vert+ r\\left \\Vert \\mathcal{R}_1 \\left( \\prod_{k=1}^{L} e^{-i H_k \\tau } \\right) \\right \\Vert \\\\\n",
    "\\leq& 2r \\left \\Vert \\mathcal{R}_1 \\left(  e^{-i |\\tau | \\sum_{k=1}^L \\Vert H_k \\Vert} \\right) \\right \\Vert \\\\\n",
    "\\leq& 2r \\left \\Vert \\mathcal{R}_1 \\left(  e^{-i |\\tau | L \\Lambda} \\right) \\right \\Vert \\\\\n",
    "\\leq& r (  \\tau L \\Lambda )^2 e^{\\vert \\tau \\vert L \\Lambda } \\\\\n",
    "=&\\frac{(  t L \\Lambda )^2}{r} e^{\\frac{\\vert t \\vert L \\Lambda}{r} }.\n",
    "\\end{aligned}\n",
    "\\tag{8}\n",
    "$$\n",
    "\n",
    "The conclusion of linear accumulation of errors in quantum circuits is used here, that is, $\\Vert U^r - V^r \\Vert \\leq r\\Vert U - V \\Vert$. Readers who are not familiar with this conclusion can refer to section 4.5.3 in [3]; and also use the conclusion when $k=1$ in formula (7). So far, we have calculated the upper bound of the simulation error of the product formula for a complete evolution time $t$, that is, the second-order term $O (t^2/r)$ in equation (4).\n",
    "\n",
    "After obtaining the upper bound of the simulation error, the lower bound of the circuit depth required to reach a certain accuracy $\\epsilon$ can be calculated. From (8), we can find that the formula contains a $L$ term, which means that as the number of Hamiltonian terms increases, the upper bound of simulation error will become larger and larger, which will cause a deeper circuit if we need to control the accuracy. qDRIFT introduced in this tutorial optimizes this problem. qDRIFT focuses on the coefficients of the Hamiltonian itself and models it as a probability distribution. Each unitary gate is sampled from the probability distribution independently and repeated a certain number of times to form a quantum circuit. Finally, under a given accuracy, the depth of the quantum circuits will not explicitly contain the number of Hamiltonian items $L$. Now we will introduce it.\n",
    "\n",
    "\n",
    "## qDRIFT\n",
    "\n",
    "First, we give the form of the target Hamiltonian\n",
    "$$\n",
    "H=\\sum_{j=1}^L h_j H_j,\n",
    "\\tag{9}\n",
    "$$\n",
    "\n",
    "It contains $L$ sub-Hamiltonian $H_j$, note that here $H_j$ has been normalized, that is, $\\Vert H_j \\Vert = 1$, where $\\Vert\\cdot\\Vert$ is the Schatten-$\\infty$ norm. $h_j$ is the coefficient of each sub-Hamiltonian, which is a positive real number. Using these coefficients, we can construct a discrete probability distribution, and take the proportion of a single coefficient in the sum of the Hamiltonian coefficients as the probability of each unitary gate being sampled, which is $p_j =h_j / \\lambda $, where $\\lambda = \\sum_j h_j $ is the sum of the coefficients. Then the sampling will be repeated $ N $ times (to compared with product formula, we let $N=Lr$ here ), we will get an ordered list arranged by $j $ and can construct a unitary gate $U_j = e^{i\\tau H_j}$ according to the arrangement. Assuming $L = 3 $, $r = 2 $, we can sample an ordered list according to the above probability distribution, as shown in\n",
    "\n",
    "$$\n",
    "[ 3, 1, 2 ,3 ,3 ,1 ],\n",
    "$$\n",
    "\n",
    "then we can construct the quantum circuits as  \n",
    "\n",
    "$$\n",
    "U_{circuit} = e^{i\\tau H_1}e^{i\\tau H_3}e^{i\\tau H_3}e^{i\\tau H_2}e^{i\\tau H_1}e^{i\\tau H_3},\n",
    "$$\n",
    "\n",
    "$\\tau = t \\lambda / N$. That is an implementation of qDRIFT to simulate Hamiltonian.\n",
    "\n",
    "The implementation process of qDRIFT is very simple, and its advantage is that the complexity of the number of unitary gates is $O ((\\lambda t) ^ 2 / \\epsilon) $) when the target precision is $\\epsilon $. It can be seen that this is a result without $L$. In other words, the number of unitary gates is not explicitly related to the number of Hamiltonian terms, which can effectively reduce the length of the simulated circuit when the number of Hamiltonian terms is large. Next, we will give a proof.\n",
    "\n",
    "We model the process of sampling from probability distribution as a quantum channel. We use the curlicue letters $\\mathcal{E} $ and $\\mathcal{U}$ to represent the channel established through qDRIFT and the channel to be simulated, and use $\\mathcal{E}_ N $ and $\\mathcal{U}_N$ to represent one of the $n $ actions of their respective channels on the quantum state $\\rho $,\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "&\\mathcal{U}_N (\\rho) = e^{\\frac{it}{N}H} \\rho e^{\\frac{-it}{N}H}= e^{\\frac{t}{N}\\mathcal{L}} (\\rho),\n",
    "\\\\\n",
    "&\\mathcal{E}_N (\\rho)=p_j e^{i\\tau H_j} \\rho e^{-i\\tau H_j}=p_j e^{\\tau \\mathcal{L}_j}(\\rho).\n",
    "\\end{aligned}\n",
    "\\tag{10}\n",
    "$$\n",
    "\n",
    "Here we introduce the Liouvillian representation, that is, for the quantum channel $\\mathcal{P} (\\rho) = e^{iHt} \\rho e^{-iHt}$, there is\n",
    "$$\n",
    "\\mathcal{P}(\\rho)=e^{iHt}\\rho e^{-iHt}=e^{t\\mathcal{L}}(\\rho)=\\sum_{k=0}^\\infty \\frac{t^k \\mathcal{L}^k (\\rho)}{k!},\n",
    "\\tag{11}\n",
    "$$\n",
    "\n",
    "where $\\mathcal{L} (\\rho) =i (H\\rho - \\rho H) $, similarly, $\\mathcal{L}_ j(\\rho)=i(H_j\\rho - \\rho H_j)$. It should be noted that the operation rules of its series follow $\\mathcal{L}^ {n+1} (\\rho) =i (H\\mathcal{L}^n (\\rho) -\\mathcal{L}^n (\\rho) H) $. Specifically, $\\mathcal{U}_N = \\sum_{n=0}^\\infty \\frac{t^n\\mathcal{L}^n}{n!N^n}$, $\\mathcal{E}_N =\\sum_{j}p_j \\sum_{n=0}^\\infty \\frac{\\lambda^n t^n \\mathcal{L}_j^n}{n!N^n}$. Next, how do we measure the distance between two channels? Here we introduce the definition of [diamond norm](https://en.wikipedia.org/wiki/Diamond_norm)\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\Vert \\mathcal{P} \\Vert_\\Diamond :=\\sup_{\\rho ; \\Vert \\rho \\Vert _1 =1}\\Vert (\\mathcal{P} \\otimes \\mathbb{I})(\\rho )\\Vert _1 .\n",
    "\\end{aligned}\n",
    "\\tag{12}\n",
    "$$\n",
    "Here $\\mathbb {I} $ is the identity channel which has the same size as $\\mathcal{P}$ and $\\Vert \\cdot \\Vert$ is Schatten-$1$ norm or called [trace norm](https://en.wikipedia.org/wiki/Schatten_norm). We use the diamond norm to define the distance between two quantum channels\n",
    "$$\n",
    "\\begin{aligned}\n",
    "d_\\Diamond (\\mathcal{E},\\mathcal{U}) &=\\frac{1}{2} \\Vert \\mathcal{E} -\\mathcal{U} \\Vert_\\Diamond\n",
    "\\\\\n",
    "&=\\sup_{\\rho ; \\Vert \\rho \\Vert _1 =1} \\frac{1}{2} \\Vert ((\\mathcal{E}-\\mathcal{U}) \\otimes \\mathbb{I})(\\rho )\\Vert _1 .\n",
    "\\end{aligned}\n",
    "\\tag{13}\n",
    "$$\n",
    "The diamond norm represents the maximum possibility that the two channels can be distinguished in all quantum states. The larger its value, the greater the possibility that the two channels can be distinguished, which means that the two channels are far away and the simulation ability is poor; on the contrary, if its value is small, it means that the simulation ability is good. Next, we can calculate the upper bound of the distance of the channel with a single action of the channel\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    \\Vert \\mathcal{U}_N-\\mathcal{E}_N \\Vert_\\Diamond &= \\left\\Vert \\sum_{n=2}^\\infty \\frac{t^n\\mathcal{L}^n}{n!N^n}-\\sum_{j}\\frac{h_j}{\\lambda} \\sum_{n=2}^\\infty \\frac{\\lambda^n t^n \\mathcal{L}_j^n}{n!N^n} \\right\\Vert_\\Diamond\\\\\n",
    "    &\\leq \\sum_{n=2}^\\infty \\frac{t^n\\Vert \\mathcal{L}^n \\Vert_\\Diamond }{n!N^n} + \\sum_{j}\\frac{h_j}{\\lambda} \\sum_{n=2}^\\infty \\frac{\\lambda^n t^n \\Vert\\mathcal{L}_j^n \\Vert_\\Diamond }{n!N^n}\\\\\n",
    "    &\\leq \\sum_{n=2}^\\infty \\frac{1}{n!}\\left( \\frac{2\\lambda t}{N}\\right)^n+\\sum_{j}\\frac{h_j}{\\lambda} \\sum_{n=2}^\\infty \\frac{1}{n!}\\left( \\frac{2\\lambda t}{N}\\right)^n\\\\\n",
    "    &=2\\sum_{n=2}^\\infty \\frac{1}{n!}\\left( \\frac{2\\lambda t}{N}\\right)^n .\n",
    "\\end{aligned}\n",
    "\\tag{14}\n",
    "$$\n",
    "\n",
    "The conclusion $\\Vert \\mathcal{L} \\Vert_ \\Diamond \\leq 2\\Vert H\\Vert \\leq 2\\lambda$ is used here. Similarly, $\\Vert \\mathcal{L}_ j \\Vert_ \\Diamond \\leq 2\\Vert H_ j\\Vert \\leq 2$ [4]. Then, we can use the conclusion of (7) mentioned above to make $k=1$, $\\alpha=2 \\lambda t /N$, and then we can get\n",
    "\n",
    "$$\n",
    "d_\\Diamond (\\mathcal{U}_N,\\mathcal{E}_N) \\leq \\frac{2\\lambda^2 t^2}{N^2} e^{2\\lambda t/N} .\n",
    "\\tag{15}\n",
    "$$\n",
    "\n",
    "Then using the conclusion of $\\Vert U^r - V^r \\Vert \\leq r\\Vert U - V \\Vert$ again (it should be noted that $U$ and $V$ in the formula are linear operators, but quantum channels $\\mathcal{U}_N$ and $\\mathcal{E}_N$ are still feasible. Refer to Chapter 3.3.2 in [6] to get more details.). With $2\\lambda t \\ll N$ generally, we can deduce\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "d_\\Diamond (\\mathcal{U},\\mathcal{E}) &\\leq N d_\\Diamond (\\mathcal{U}_N, \\mathcal{E}_N)\\\\\n",
    "    &=\\frac{2\\lambda^2 t^2}{N} e^{2\\lambda t/N} \\approx \\frac{2\\lambda^2 t^2}{N}.\n",
    "\\end{aligned}\n",
    "\\tag{16}\n",
    "$$\n",
    "\n",
    "So $N \\sim O ((\\lambda t) ^2 /\\epsilon) $. It can be seen from the above formula that under the condition of $\\lambda \\ll \\Lambda L $ (Hamiltonian is defined as $H=\\sum_{j=1}^L h_j H_j$ in qDRIFT, so $\\Lambda = \\max_k h_k $ here), the distance will not be explicitly related to $L$, which means that when $L$ is large, that is, the situation is complex, the quantum circuit will not be increased, and the number of unitary gates can be effectively controlled. Many systems satisfy the condition $\\lambda \\ll \\Lambda L$ such as carbon-dioxide, ethane. But not all cases satisfy. If $\\lambda = \\Lambda L$ or $\\lambda = \\Lambda \\sqrt{L}$, their upper bounds are $O (L^2 (\\Lambda t) ^2 /\\epsilon) $ and $O (L (\\Lambda t) ^2 /\\epsilon) $ respectively, we can see that they still increase with the number of Hamiltonian terms. Interested readers can refer to [4] for more details.\n",
    "\n",
    "\n",
    "## Code Demonstration\n",
    "We will implement qDRIFT in combination with code. We will first demonstrate the performance of its sampling results, and then calculate the simulation error of its channel. First, we need to import the required packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "\n",
    "import math\n",
    "import numpy as np                    \n",
    "import scipy                                          \n",
    "import paddle_quantum as pq         \n",
    "import paddle\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")   # Hide warnings\n",
    "np.set_printoptions(suppress=True,linewidth=np.nan)        # Enable full display, so that line breaks would not appear\n",
    "                                                           # when viewing the matrix on the terminal print\n",
    "pq.set_backend('density_matrix')    # use density matrix representation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We assume that the system consists of two qubits. We can use the `hamiltonian` module of the  Paddle Quantum to construct a Hamiltonian with $L=4$ satisfied $\\lambda \\ll \\Lambda L$ to demonstrate, which is our target Hamiltonian, as follows\n",
    "$$\n",
    "\\begin{aligned}\n",
    "H&=I \\otimes X + 0.05 * X \\otimes Z + 0.05 * I \\otimes Y+0.05 * X \\otimes X   .\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Target Hamiltonian is :\n",
      " [[ 0.  +0.j    1.  -0.05j  0.05+0.j    0.05+0.j  ]\n",
      " [ 1.  +0.05j  0.  +0.j    0.05+0.j   -0.05+0.j  ]\n",
      " [ 0.05+0.j    0.05+0.j    0.  +0.j    1.  -0.05j]\n",
      " [ 0.05+0.j   -0.05+0.j    1.  +0.05j  0.  +0.j  ]]\n"
     ]
    }
   ],
   "source": [
    "qubits = 2  # Set the number of qubits\n",
    "H_j = [(1.0, 'I0,X1'),  # The Pauli string of target Hanmiltonian\n",
    "       (0.05, 'X0,Z1'),\n",
    "       (0.05, 'I0,Y1'),\n",
    "       (0.05, 'X0,X1'), ]\n",
    "\n",
    "H = pq.hamiltonian.Hamiltonian(H_j)    \n",
    "print(f'Target Hamiltonian is :\\n {H.construct_h_matrix(qubit_num=qubits)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we calculate probability according to $\\lambda = \\sum_ j h_ j$，$ p_ j=h_j/\\lambda $. In this experiment, we assume that our target accuracy $\\epsilon=0.1$, simulation time $t=1$, that is, we need to sample $n=\\lceil \\frac{2\\lambda^2 t^2} {\\epsilon}\\rceil = 27 $ times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "To meet the accuracy of 0.1, there are 27 unitary gates needed.\n"
     ]
    }
   ],
   "source": [
    "h_j = np.array(H.coefficients)  # Get coeffcients\n",
    "lamda = h_j.sum()\n",
    "p_j = h_j/lamda  # Calculate the discrete probability distribution\n",
    "accuracy = 0.1\n",
    "t = 1\n",
    "gate_counts = math.ceil(2 * lamda**2 * t**2 / accuracy)\n",
    "\n",
    "print(f'To meet the accuracy of {accuracy}, there are {gate_counts} unitary gates needed.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we will sample 27 times independently from the probability distribution $p_j$ and construct a unitary circuit according to the sample results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample results:\n",
      " [1 1 1 1 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 4 1]\n",
      "The simulation circuit matrix of qDRIFT is: \n",
      " [[ 0.51998969-0.02531459j  0.03408183+0.85099285j -0.03524884+0.02376227j -0.0353899 +0.02366261j]\n",
      " [-0.03408183+0.85099285j  0.51998969+0.02531459j  0.0353899 +0.02366261j -0.03524884-0.02376227j]\n",
      " [-0.03524884+0.02376227j -0.0353899 +0.02366261j  0.51998969-0.02531459j  0.03408183+0.85099285j]\n",
      " [ 0.0353899 +0.02366261j -0.03524884-0.02376227j -0.03408183+0.85099285j  0.51998969+0.02531459j]] \n",
      "The original circuit matrix is: \n",
      " [[ 0.53752508-0.00075235j  0.04202098+0.83966719j -0.04201839+0.04202098j -0.00075235+0.02697398j]\n",
      " [-0.04202098+0.83966719j  0.53752508+0.00075235j  0.00075235+0.02697398j -0.04201839-0.04202098j]\n",
      " [-0.04201839+0.04202098j -0.00075235+0.02697398j  0.53752508-0.00075235j  0.04202098+0.83966719j]\n",
      " [ 0.00075235+0.02697398j -0.04201839-0.04202098j -0.04202098+0.83966719j  0.53752508+0.00075235j]]\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(666)  # Fix the random seed to demonstrate\n",
    "sample_list = np.random.choice(a=range(1, 5), size=gate_counts, replace=True, p=p_j)\n",
    "print(f'Sample results:\\n {sample_list}')\n",
    "\n",
    "# Calculate the sampled unitary circuit according to the sample result\n",
    "simulation = np.identity(2 ** qubits)  # Generate the identity matrix\n",
    "tau = 1j*lamda*t/gate_counts\n",
    "for i in sample_list:\n",
    "    pauli_str_j = (1.0, H_j[i-1][1])   # Get H_ j. Note that its original coefficient should be discarded\n",
    "    H_i = pq.hamiltonian.Hamiltonian([pauli_str_j]).construct_h_matrix(qubit_num=qubits)\n",
    "    simulation = np.matmul(scipy.linalg.expm(tau*H_i), simulation)  \n",
    "origin = scipy.linalg.expm(1j*t*H.construct_h_matrix(qubit_num=qubits))  # Calculate the original circuit of the target Hamiltonian\n",
    "print(f'The simulation circuit matrix of qDRIFT is: \\n {simulation} \\nThe original circuit matrix is: \\n {origin}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can calculate the simulation error $\\Vert e^{iHt}-U_{circuit} \\Vert $ between the unitary circuit sampled from qDRIFT and the original circuit, note that the norm here is the spectral norm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulation error: 0.0309\n"
     ]
    }
   ],
   "source": [
    "distance = 0.5*np.linalg.norm(origin-simulation, ord=2)\n",
    "print(f'Simulation error: {distance:.4f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, we can take a specific quantum state to test. Without losing generality, we assume that the initial quantum state is a zero state, that is, $\\rho(0)  = | 0 \\rangle \\langle 0 | $. We can let the quantum state evolve through the original circuit and the simulation circuit respectively. At the time of $t $, the quantum state is $\\rho(t)_{origin}$ and $\\rho(t)_{qDRIFT}$, we can use the fidelity between these two quantum states to measure the effect of simulation circuits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial State: \n",
      " [[1.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j]]\n",
      "The fidelity between two states is 0.9989\n"
     ]
    }
   ],
   "source": [
    "rho_0 = pq.state.zero_state(qubits).numpy() # Generate the zero state density matrix\n",
    "print(f'Initial State: \\n {rho_0}')\n",
    "\n",
    "rho_t_origin = pq.state.to_state(origin @ rho_0 @ origin.T.conjugate())  # Evolve through the original circuit\n",
    "rho_t_qdrift = pq.state.to_state(simulation @ rho_0 @ simulation.T.conjugate())  # Evolve through the simulation circuit\n",
    "fidelity = pq.qinfo.state_fidelity(rho_t_origin, rho_t_qdrift)\n",
    "print(f'The fidelity between two states is {float(fidelity):.4f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can be found that the above tests meet our accuracy requirements. However, different from a sampled unitary circuit, we regard the qDRIFT sampling method as a quantum channel, that is, a mapping of the quantum state $\\rho $. The above experiment is only a specific instance of this channel. Next, we will analyze the performance of this channel. First, we can define a function to describe the qDRIFT channel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define qDRIFT channel\n",
    "def qdrift_channel(iter_num, sample_num, hamiltonian_list, coefficient_list, simulation_time, qubits, input_state):\n",
    "    '''\n",
    "    Input :\n",
    "        iter_num : the current number of iterations, as a label of recursion\n",
    "        sample_num : number of samples, N\n",
    "        hamiltonian_list : the list of Pauli string of the target Hamiltonian, H_j\n",
    "        coefficient_list : the coefficient list of sub-Hamiltonian, h_j\n",
    "        simulation_time : simulation time, t\n",
    "        qubits : the number of qubits \n",
    "        input_state : the input quantum state, which should be a density matrix\n",
    "    \n",
    "    Return :\n",
    "        The quantum state aftre the evolution of this channel, represented in density matrix\n",
    "    '''\n",
    "    lamda = coefficient_list.sum() \n",
    "    tau = lamda*simulation_time/sample_num\n",
    "    output = 0\n",
    "\n",
    "    if iter_num != 1:   # Enable recursion when iteration flag is not 1\n",
    "        input_state = qdrift_channel(iter_num-1, sample_num, hamiltonian_list,\n",
    "                                     coefficient_list, simulation_time, qubits, input_state)\n",
    "\n",
    "    # Calculate e^{iH\\tau} \\rho e^{-iH\\tau}                                 \n",
    "    for sub_H, sub_h in zip(hamiltonian_list, coefficient_list):\n",
    "        sub_H = pq.hamiltonian.Hamiltonian([sub_H]).construct_h_matrix(qubit_num=qubits)\n",
    "        unitary = scipy.linalg.expm(1j*tau*sub_H)  # Calculate e^{iH\\tau}\n",
    "        output += sub_h/lamda*unitary @ input_state @ unitary.conjugate().T\n",
    "    return output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can calculate the distance between the two channels through the diamond norm, but the solution of the diamond norm is a positive semidefinite programming problem, that is\n",
    "\n",
    "$$\n",
    "d_\\Diamond(\\mathcal{U}- \\mathcal{E})=\\sup_{\\Omega \\geq 0 \\atop \\rho \\geq 0}\\{\\text{Tr}[\\Omega (\\Gamma_\\mathcal{U}-\\Gamma_\\mathcal{E})]: \\Omega \\leq \\rho \\otimes \\mathbb{I},\\text{Tr} (\\rho)=1\\},\n",
    "\\tag{17}\n",
    "$$\n",
    "where $\\Gamma_ \\mathcal{U} $ and $\\Gamma_ \\mathcal{E}$ are Choi representations of the original channel and the simulation channel. There are many forms of positive semidefinite programming and Choi representation of diamond norm, and interested readers can read [6-8] for more details. The Choi representation we use here is:\n",
    "$$\n",
    "\\Gamma_\\mathcal{P}=\\sum_{i,j=0}^{d-1} |i\\rangle \\langle j| \\otimes \\mathcal{P}(|i\\rangle \\langle j|),\n",
    "\\tag{18}\n",
    "$$\n",
    "where $\\mathcal {P} $ is the quantum channel and $d$ is the dimension of the input quantum state of the quantum channel. Here we first calculate the Choi representation of the two channels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# There is how to calculate the Choi representation of the original channel and the qDRIFT channel, \n",
    "# and the diamond norm can be calculated under this representation\n",
    "choi_qdrift = 0\n",
    "choi_origin = 0\n",
    "channel = scipy.linalg.expm(1j*t*H.construct_h_matrix(qubit_num=qubits))\n",
    "for i in range(2 ** qubits):\n",
    "    for k in range(2 ** qubits):\n",
    "        choi_temp = np.zeros((2 ** qubits, 2 ** qubits))\n",
    "        choi_temp[i][k] = 1  # Generate |i\\rangle \\langle k|\n",
    "\n",
    "        # Calculate the Choi matrix of channel E\n",
    "        # Calculate \\mathcal{E}(|i\\rangle \\langle k|)\n",
    "        choi_temp_qdrift = qdrift_channel(gate_counts, gate_counts, H_j, h_j, t, qubits, choi_temp)  \n",
    "        # Calculate |i\\rangle \\langle k| \\otimes \\mathcal{E}(|i\\rangle \\langle k|)\n",
    "        choi_qdrift += np.kron(choi_temp, choi_temp_qdrift)\n",
    "\n",
    "        # Calculate the Choi matrix of channel U\n",
    "        # Calculate \\mathcal{U}(|i\\rangle \\langle k|)\n",
    "        choi_temp_origin = channel @ choi_temp @ channel.T.conjugate()\n",
    "        # Calculate |i\\rangle \\langle k| \\otimes \\mathcal{U}(|i\\rangle \\langle k|)\n",
    "        choi_origin += np.kron(choi_temp, choi_temp_origin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can calculate the diamond norm according to formula (17) and find the diamond distance of the two channels. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The distance between the two channels is: 0.0764\n"
     ]
    }
   ],
   "source": [
    "diamond_distance = 0.5 * pq.qinfo.diamond_norm(paddle.to_tensor(choi_origin-choi_qdrift))\n",
    "print(f'The distance between the two channels is: {diamond_distance:.4f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The calculation results are in line with our expectations. Note that this value represents the expectance of the worst performance of the channel sampling for a specific simulation circuit instance, so it can not guarantee that each sampled circuit can achieve this accuracy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "Quantum simulation itself is a relatively broad topic, and its application is also very extensive. This tutorial introduces the basic theory of product formula and the qDRIFT method, but qDRIFT is not the only method for random product formulas. As a branch of the method of quantum simulation using product formula, random product formula has many methods worth exploring."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Reference\n",
    " \n",
    "[1] Lloyd, Seth. \"Universal quantum simulators.\" [Science (1996): 1073-1078](https://www.jstor.org/stable/2899535).\n",
    "\n",
    "[2] Childs, Andrew M., et al. \"Toward the first quantum simulation with quantum speedup.\" [Proceedings of the National Academy of Sciences 115.38 (2018): 9456-9461](https://www.pnas.org/content/115/38/9456.short).\n",
    "\n",
    "[3] Nielsen, Michael A., and Isaac Chuang. \"Quantum computation and quantum information.\" (2002): 558-559.\n",
    "\n",
    "[4] Campbell, E. . \"Random Compiler for Fast Hamiltonian Simulation.\" [Physical Review Letters 123.7(2019):070503.1-070503.5](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.123.070503).\n",
    "\n",
    "[5] Khatri, Sumeet, and Mark M. Wilde. \"Principles of quantum communication theory: A modern approach.\" [arXiv preprint arXiv:2011.04672 (2020).](https://arxiv.org/abs/2011.04672)\n",
    "\n",
    "[6] Watrous, J. . [The Theory of Quantum Information](https://cs.uwaterloo.ca/~watrous/TQI/).  2018.\n",
    "\n",
    "[7] Watrous, J. . \"Simpler semidefinite programs for completely bounded norms.\" [Chicago Journal of Theoretical Computer Science (2012).](https://arxiv.org/abs/1207.5726)\n",
    "\n",
    "[8] Watrous, J. . \"Semidefinite Programs for Completely Bounded Norms.\" [Theory of Computing 5.1(2009):217-238.](https://arxiv.org/abs/0901.4709)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.7.13 ('py3.7_pq2.2.1')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.13"
  },
  "vscode": {
   "interpreter": {
    "hash": "4e4e2eb86ad73936e915e7c7629a18a8ca06348106cf3e66676b9578cb1a47dd"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}