tutorial/quantum_simulation/ClassicalShadow_Application_EN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4f560a55",
   "metadata": {},
   "source": [
    "# Estimation of Quantum State Properties Based on the Classical Shadow"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a894cbba",
   "metadata": {},
   "source": [
    "<em> Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved. </em>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4cdaa03",
   "metadata": {},
   "source": [
    "## Overview"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b3ea48a",
   "metadata": {},
   "source": [
    "In [The Classical Shadow of Unknown Quantum States](./ClassicalShadow_Intro_EN.ipynb), we introduced some theoretical knowledge of the classical shadow and showed how to construct the classical shadow of a quantum state $\\rho$. According to the theoretical derivation of the classical shadow, it is very suitable for estimating the linear properties of quantum states. At present, its basic applications are: quantum fidelity estimation, entanglement verification, the estimation of local observables' expectation values, and the estimation of global observables' expectation values [1]. Among them, the estimation of observables' expectation values widely appears in current quantum algorithms, such as [Variational Quantum Eigensolver](./VQE_EN.ipynb)  (VQE) for estimating the ground state energy of complex molecular Hamiltonian. Next, we will focus on the estimation of observables' expectation values algorithms based on the classical shadow and show how to use the shadow function in Paddle Quantum to do so."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db8b954f",
   "metadata": {},
   "source": [
    "## Estimation of Observables' Expectation Values"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2869ea5",
   "metadata": {},
   "source": [
    "### Problem description"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6adea50",
   "metadata": {},
   "source": [
    "In the field of quantum chemistry, one of the core tasks is to solve for ground state energy of the Hamiltonian $\\hat{H}$ on a closed physical system on a quantum scale and its corresponding ground state. The main method is to prepare a parameterized trial wave function $|\\Psi(\\theta)\\rangle$. Then, the parameter $\\theta$ is being continuously adjusted and optimized to minimize the expected value $\\langle\\Psi(\\theta)|\\hat{H}| \\Psi(\\theta)\\rangle$ using classical optimization algorithms (such as gradient descent method). The principle of this scheme is based on Rayleigh-Ritz variational principle,\n",
    "\n",
    "$$\n",
    "E_{0}=\\min _{\\theta}\\langle\\Psi(\\theta)|\\hat{H}| \\Psi(\\theta)\\rangle \\tag{1}\n",
    "$$\n",
    "\n",
    "where $E_{0}$ represents the ground state energy of the system. Numerically, the problem can be understood as solving for the minimum eigenvalue $\\lambda_{\\min }$ of a discretized Hamiltonian $\\hat{H}$ (Hermitian matrix) and its corresponding eigenvector $\\left|\\Psi_{0}\\right\\rangle$. Where the classical shadow comes into play is to calculate the $\\langle\\Psi(\\theta)|\\hat{H}| \\Psi(\\theta)\\rangle = \\operatorname{tr}(\\hat{H}\\rho )$ part in every optimization iteration where $ \\rho = | \\Psi(\\theta)\\rangle\\langle\\Psi(\\theta)|$.\n",
    "\n",
    "The problem is then transformed into: for a quantum state $\\rho$ of $n$ qubits and an observable (Hamiltonian) $\\hat{H}$ that can be written as a linear combination of a set of Pauli operators $\\{I, X, Y, Z\\}^{\\otimes n}$,\n",
    "\n",
    "$$\n",
    "\\hat{H}=\\sum_{Q \\in\\{I, X, Y, Z\\} ^{\\otimes n}} \\alpha_{Q} Q \\quad \\text{where} \\quad \\alpha_{Q} \\in \\mathbb{R}, \\tag{2}\n",
    "$$\n",
    "\n",
    "how to estimate the observable's expectation value $\\operatorname{tr}(\\hat{H}\\rho )$ using the classical shadow?\n",
    "\n",
    "\n",
    "The most intuitive method is to use each term of the Hamiltonian as a measurement base, and the corresponding Pauli measurements are made for the quantum state $\\rho$. A certain number of repetitions are made for each of the measurements, and then the measurement results are being processed to obtain the estimation value. Here we refer to this method as the item-by-item measurement method.\n",
    "\n",
    "Readers can see that when both $n$ and the number of terms of the Hamiltonian $\\hat{H}$ are small, we can get $\\operatorname{tr}(\\hat{H}\\rho )$ through the item-by-item measurement method. However, when $n$ is large and the number of terms of $\\hat{H}$ increases, the cost of the this method will increase significantly. The classical-shadow-based method that will be introduced can obtain the same precision estimation of $\\operatorname{tr}(\\hat{H}\\rho )$ with less cost."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff8af318",
   "metadata": {},
   "source": [
    "### The improved algorithm based on the classical shadow"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3db3328",
   "metadata": {},
   "source": [
    "When constructing the classical shadow, a critical step is to uniformly and randomly sample the unitary transformation from the fixed set. In [The Classical Shadow of Unknown Quantum States](./ClassicalShadow_Intro_EN.ipynb), we showed the case when the selected set is the Clifford group. When the selected set is a Clifford group on single qubit, the process of sampling and measurement is equivalent to making Pauli measurements on the quantum states. The classical shadow algorithm using random Pauli measurements (CS) is provided in Paddle Quantum. Briefly, in the CS algorithm, we repeatedly choose a Pauli basis for each qubit uniformly at random to measure the quantum state $\\rho$ and estimate the observable's expectation value based on the measurement results, and the reader may refer to [1-2] to learn the details of the principle. Further, when the Pauli measurement base is no longer selected uniformly and randomly, improved algorithms are proposed [2-3]. Relevant algorithm functions are also provided in Paddle Quantum: Locally-biased classical shadows (LBCS) [2], Adaptive Pauli shadows (APS) [3]. Readers can refer to [1-3] to learn these algorithms in detail."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b772de7",
   "metadata": {},
   "source": [
    "## Paddle Quantum Implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e094e752",
   "metadata": {},
   "source": [
    "In Paddle Quantum, we provide the shadow function, which mainly includes two parts to use the above three algorithms to estimate the expectation value of the observable and obtain the classical shadow data of unknown quantum states. Next, we will show how to implement finding the ground state energy estimation of $H_{2}$ and $LiH$ based on the shadow function in Paddle Quantum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "efc769b4",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Import all the dependencies\n",
    "import numpy as np\n",
    "from numpy import pi as PI\n",
    "import paddle\n",
    "from paddle_quantum.circuit import UAnsatz\n",
    "from paddle_quantum.VQE.chemistrysub import H2_generator\n",
    "from paddle_quantum.utils import Hamiltonian"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a58e8871",
   "metadata": {},
   "source": [
    "###  Estimate the ground state energy of hydrogen molecule  ($H_{ 2}$)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebda5aeb",
   "metadata": {},
   "source": [
    "Import the Hamiltonian of $H_2$ with 4 qubits (for details, please refer to [Variational Quantum Eigensolver](./VQE_CN.ipynb) tutorial to obtain $H_2$ molecular Hamiltonian)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "70b8fdef",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H2 hamiltonian =  -0.04207897647782277 I0\n",
      "0.17771287465139946 Z0\n",
      "0.1777128746513994 Z1\n",
      "-0.2427428051314046 Z2\n",
      "-0.24274280513140462 Z3\n",
      "0.17059738328801055 Z0, Z1\n",
      "0.04475014401535163 Y0, X1, X2, Y3\n",
      "-0.04475014401535163 Y0, Y1, X2, X3\n",
      "-0.04475014401535163 X0, X1, Y2, Y3\n",
      "0.04475014401535163 X0, Y1, Y2, X3\n",
      "0.12293305056183797 Z0, Z2\n",
      "0.1676831945771896 Z0, Z3\n",
      "0.1676831945771896 Z1, Z2\n",
      "0.12293305056183797 Z1, Z3\n",
      "0.1762764080431959 Z2, Z3\n"
     ]
    }
   ],
   "source": [
    "# Set up the Hamiltonian of hydrogen molecule\n",
    "H2_pauli_str, H2_qubit = H2_generator()\n",
    "# Construct a Hamiltonian class instance using the H2_pauli_str\n",
    "H2_hamiltonian = Hamiltonian(H2_pauli_str)\n",
    "print('H2 hamiltonian = ', H2_hamiltonian)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d6a6532",
   "metadata": {},
   "source": [
    "To show how to estimate the ground state energy using the classical-shadow-based algorithm, we first get the estimated ground state of $H_{2}$ using the VQE algorithm in Paddle Quantum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "1bc042c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def U_theta(theta, hamiltonian, N, D):\n",
    "    \"\"\"\n",
    "    Quantum Neural Network\n",
    "    \"\"\"\n",
    "    \n",
    "    # Initialize the quantum neural network according to the number of qubits N\n",
    "    cir = UAnsatz(N)\n",
    "    \n",
    "    # Built-in {R_y + CNOT} circuit template\n",
    "    cir.real_entangled_layer(theta[:D], D)\n",
    "    \n",
    "    # Lay R_y gates in the last row\n",
    "    for i in range(N):\n",
    "        cir.ry(theta=theta[D][i][0], which_qubit=i)\n",
    "        \n",
    "    # The quantum neural network acts on the default initial state |0000>\n",
    "    cir.run_state_vector()\n",
    "    \n",
    "    # Calculate the expected value of a given Hamiltonian\n",
    "    expectation_val = cir.expecval(hamiltonian)\n",
    "\n",
    "    return expectation_val, cir\n",
    "\n",
    "class StateNet(paddle.nn.Layer):\n",
    "    \"\"\"\n",
    "    Construct the model net\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, shape, dtype=\"float64\"):\n",
    "        super(StateNet, self).__init__()\n",
    "        \n",
    "        # Assign the theta parameter list to be the trainable parameter list of the circuit\n",
    "        self.theta = self.create_parameter(shape=shape, \n",
    "                                           default_initializer=paddle.nn.initializer.Uniform(low=0.0, high=2*PI),\n",
    "                                           dtype=dtype, is_bias=False)\n",
    "        \n",
    "    # Define loss function and forward propagation mechanism\n",
    "    def forward(self, hamiltonian, N, D):\n",
    "        \n",
    "        # Calculate the loss function/expected value\n",
    "        loss, cir = U_theta(self.theta, hamiltonian, N, D)\n",
    "\n",
    "        return loss, cir"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c53b204",
   "metadata": {},
   "source": [
    "After building the quantum neural network of the VQE algorithm and defining the loss function, we can estimate the ground state energy of $H_{2}$ and the quantum circuit corresponding to the ground state by training the quantum neural network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "40036ae4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter: 20 loss: -1.0865\n",
      "iter: 20 Ground state energy: -1.0865 Ha\n",
      "iter: 40 loss: -1.1284\n",
      "iter: 40 Ground state energy: -1.1284 Ha\n",
      "iter: 60 loss: -1.1355\n",
      "iter: 60 Ground state energy: -1.1355 Ha\n",
      "iter: 80 loss: -1.1360\n",
      "iter: 80 Ground state energy: -1.1360 Ha\n",
      "\n",
      "The trained circuit:\n",
      "--Ry(-1.57)----*--------------x----Ry(4.713)----*--------------x----Ry(3.129)--\n",
      "               |              |                 |              |               \n",
      "--Ry(5.091)----x----*---------|----Ry(1.768)----x----*---------|----Ry(3.711)--\n",
      "                    |         |                      |         |               \n",
      "--Ry(0.997)---------x----*----|----Ry(4.696)---------x----*----|----Ry(3.180)--\n",
      "                         |    |                           |    |               \n",
      "--Ry(4.753)--------------x----*----Ry(4.701)--------------x----*----Ry(-0.37)--\n",
      "                                                                               \n"
     ]
    }
   ],
   "source": [
    "ITR = 80  # Set the number of optimization iterations\n",
    "LR = 0.4   # Set the learning rate\n",
    "D = 2      # Set the depth of the repetitive calculation module in QNN\n",
    "N = H2_hamiltonian.n_qubits \n",
    "\n",
    "# Determine the parameter dimension of the network\n",
    "net = StateNet(shape=[D + 1, N, 1])\n",
    "\n",
    "# Generally speaking, we use Adam optimizer to obtain relatively good convergence,\n",
    "# You can change it to SGD or RMS prop.\n",
    "opt = paddle.optimizer.Adam(learning_rate=LR, parameters=net.parameters())\n",
    "\n",
    "# Record optimization results\n",
    "summary_iter, summary_loss = [], []\n",
    "\n",
    "# Optimization loop\n",
    "for itr in range(1, ITR + 1):\n",
    "\n",
    "    # Forward propagation to calculate loss function\n",
    "    loss, cir = net(H2_hamiltonian, N, D)\n",
    "\n",
    "    # Use back propagation to minimize the loss function\n",
    "    loss.backward()\n",
    "    opt.minimize(loss)\n",
    "    opt.clear_grad()\n",
    "\n",
    "    # Record optimization results\n",
    "    summary_loss.append(loss.numpy())\n",
    "    summary_iter.append(itr)\n",
    "\n",
    "    # Print result\n",
    "    if itr % 20 == 0:\n",
    "        print(\"iter:\", itr, \"loss:\", \"%.4f\" % loss.numpy())\n",
    "        print(\"iter:\", itr, \"Ground state energy:\", \"%.4f Ha\" \n",
    "                                            % loss.numpy())\n",
    "    if itr == ITR:\n",
    "        print(\"\\nThe trained circuit:\") \n",
    "        print(cir)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6c9d23f",
   "metadata": {},
   "source": [
    "#### Introduction to shadow function"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9891733",
   "metadata": {},
   "source": [
    "At this point, we've obtained the quantum circuit for generating the ground state of $H_{2}$. We can run the `shadow_trace` function on this circuit to obtain the estimated ground state energy. In the `shadow_trace` function, our inputs are the Hamiltonian to be estimated, the number of samples, and the method selected. You can choose the sampling mode you want by specifying the parameter `method`.  Among them, CS has a broader application range and the fastest speed, but its estimation accuracy may be slightly poor; LBCS has a higher accuracy, but it runs slower as the number of terms of Hamiltonian grows; APS also has higher accuracy, but it runs slower when the number of qubits is large."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3b31e57d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H2 ground state energy =  -1.1360331063822373\n",
      "H2 ground state energy CS=  -1.1652093330967717\n",
      "H2 ground state energy LBCS=  -1.1213982275348622\n",
      "H2 ground state energy APS=  -1.137720214050516\n"
     ]
    }
   ],
   "source": [
    "# The actual energy value corresponding to the estimated ground state\n",
    "H2_energy = cir.expecval(H2_hamiltonian).numpy()[0]\n",
    "\n",
    "# Sampling times\n",
    "sample = 1500 \n",
    "# Three algorithms are used to estimate the expectation value of observable.\n",
    "H2_energy_CS = cir.shadow_trace(H2_hamiltonian, sample, method=\"CS\")\n",
    "H2_energy_LBCS = cir.shadow_trace(H2_hamiltonian, sample, method=\"LBCS\")\n",
    "H2_energy_APS = cir.shadow_trace(H2_hamiltonian, sample, method=\"APS\")\n",
    "\n",
    "print('H2 ground state energy = ', H2_energy)\n",
    "print('H2 ground state energy CS= ', H2_energy_CS)\n",
    "print('H2 ground state energy LBCS= ', H2_energy_LBCS)\n",
    "print('H2 ground state energy APS= ', H2_energy_APS)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cb856da",
   "metadata": {},
   "source": [
    "Now let's use the item-by-item measurement method to estimate the ground state energy of $H_{2}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3ae0242c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H2 ground state energy traditional =  -1.1367622727687419\n"
     ]
    }
   ],
   "source": [
    "# Use the item-by-item measurement method to estimate the ground state energy\n",
    "H2_energy_traditional = cir.expecval(H2_hamiltonian, shots=100).numpy()[0]\n",
    "print('H2 ground state energy traditional = ', H2_energy_traditional)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2e1ad02",
   "metadata": {},
   "source": [
    "We can see that under 1500 samples, the estimated ground state energy by these three algorithms are very close to the energy of the estimated ground state by VQE. The item-by-item measurement method is to make 100 measurements for each item of the Hamiltonian, and there are 15 items of the Hamiltonian of $H_{2}$, which is equivalent to 1500 measurements in total. The difference between the obtained results and the energy of the estimated ground state by VQE is also tiny. In this small-scale situation, the classical-shadow-based algorithms do not show significant advantages over the item-by-item measurement method. But in large-scale qubits scenarios, this algorithm requires only constant-level growth of the number of the Hamiltonian terms. In contrast, the item-by-item measurement method or some other methods require polynomial-level or even exponential-level growth in the number of samples to get the same accuracy[1]. In fact, it is pointed out in [2] that for CS algorithm and LBCS algorithm, the average error of estimation $\\epsilon$, variance $\\operatorname {var}(\\nu)$ and number of samples $S$ are related as follows:\n",
    "\n",
    "$$\n",
    "S = O(\\epsilon^{-2} \\operatorname{var}(\\nu) ),  \\tag{3}\n",
    "$$\n",
    "\n",
    "where the variance $\\operatorname{var} (\\nu) $ is independent of the number of samples, and the variance is related to the number of terms of the Hamiltonian. Therefore, knowing our expected precision, we can calculate the required number of samples. Similarly, our average error can also be defined according to the number of samples: \n",
    "\n",
    "$$\n",
    "\\epsilon = \\sqrt{\\frac{\\operatorname{var}}{S}}.  \\tag{4}\n",
    "$$\n",
    "\n",
    "We can be seen that in the above experiments, the errors obtained by the CS algorithm and LBCS algorithm are within the accuracy of theoretical estimation in reference [3]."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e0e174d",
   "metadata": {},
   "source": [
    "At the same time, Paddle Quantum provides a sampling function ` shadow_sample`, which supports pre-sampling of unknown quantum states, which is convenient for readers to explore other applications of the classical shadow. The specific usage is as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "fc6b465c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('yxxz', '0010'), ('xyzx', '0101'), ('yxyx', '1010'), ('yxyy', '1101'), ('yyyx', '0101'), ('xyxz', '0010'), ('zyxz', '1010'), ('xzyz', '0110'), ('zzzx', '1101'), ('xyxx', '0101')]\n"
     ]
    }
   ],
   "source": [
    "from paddle_quantum.shadow import shadow_sample\n",
    "\n",
    "# Run the circuit in the vector form\n",
    "H2_rho = np.array(cir.run_state_vector())\n",
    "# Get the data of classical shadow and output it in the form of a list\n",
    "H2_sample_data_CS = shadow_sample(H2_rho, H2_qubit, sample_shots=10, mode='state_vector', \n",
    "                                  hamiltonian=H2_hamiltonian, method='CS')\n",
    "print(H2_sample_data_CS)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e74ffd18",
   "metadata": {},
   "source": [
    "###  Estimate the ground state energy of lithium hydride  ($LiH$)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4480cb6",
   "metadata": {},
   "source": [
    "Next, we consider the ground state energy of $LiH$. First, we load from a pre-computed file to generate the molecular Pauli Hamiltonian of $LiH$ with 12 qubits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6916c8a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('./LiH_hamiltonian.txt', 'r') as lih_file:\n",
    "    unprocessed_pauli_str = lih_file.read()\n",
    "    LiH_pauli_str = [term.split(maxsplit=1) for term in unprocessed_pauli_str.split('\\n')]\n",
    "    LiH_pauli_str = [[float(term[0]), term[1]] for term in LiH_pauli_str]\n",
    "    LiH_hamiltonian = Hamiltonian(LiH_pauli_str)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61d8fd59",
   "metadata": {},
   "source": [
    "Then, we also use the VQE algorithm to obtain the ground state energy and the ground state of $LiH $ and then use the classical shadow algorithm to estimate the ground state energy. Due to the large size of the $LiH$ molecular Hamiltonian, it will take a long time to train the VQE circuit. Hence we provide a set of pre-trained parameters, using which the users could test the classical-shadow-based algorithms on $LiH$ Hamiltonian directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "fe326409",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The pre-trained VQE gets a ground state energy of: -7.7720 \n"
     ]
    }
   ],
   "source": [
    "# Load the pre-trained parameters\n",
    "pretrained_parameters = paddle.load('LiH_VQE_parameters.pdtensor')\n",
    "N = LiH_hamiltonian.n_qubits\n",
    "# Run the VQE circuit with pre-trained parameters\n",
    "energy, cir = U_theta(pretrained_parameters, LiH_hamiltonian, N, D)\n",
    "print('The pre-trained VQE gets a ground state energy of: %.4f ' % energy.numpy())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f8ff501",
   "metadata": {},
   "source": [
    "Once we have the circuit corresponding to the $LiH$ ground state, we can directly use the `shadow_trace` function to  perform random measurements. Also, since this molecular Hamiltonian has 631 terms, we specify `sample = 1262` for the function `shadow_trace` and `shots = 2` for the function `expecval` in order to ensure that the number of measurements is the same for both types of methods.\n",
    "\n",
    "Since ground state of $LiH$ has 12 qubits, there are $3^{12}$ possible combinations of measurements when doing different Pauli measurements on the ground state of $LiH$. So it is too random to just perform 1262 samples to get the valuation. Thus, we run each of the above four methods 20 times. Then take the mean of these 20 samples of data as the estimation for each algorithm, and calculate the sample variance for a simple comparison of the algorithms.（It may take at least 1 hour to run the following code blocks.）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "8f3684ce",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LiH ground state energy =  -7.771980394176657\n",
      "ave LiH ground state energy CS =  -7.835791570579005\n",
      "ave LiH ground state energy LBCS =  -7.7622296662623445\n",
      "ave LiH ground state energy APS =  -7.762836542787509\n",
      "ave LiH ground state energy traditional =  -7.8964746269601465\n",
      "time =  4206.216086864471\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "\n",
    "begin = time.time()\n",
    "estimator_CS = []\n",
    "estimator_LBCS = []\n",
    "estimator_APS = []\n",
    "estimator_traditional = []\n",
    "\n",
    "# The actual energy value corresponding to the estimated ground state\n",
    "LiH_energy = cir.expecval(LiH_hamiltonian).numpy()[0]\n",
    "\n",
    "# Number of repetition times\n",
    "n = 20 \n",
    "\n",
    "for i in range(n):\n",
    "    LiH_energy_CS = cir.shadow_trace(LiH_hamiltonian, 1262, method=\"CS\")\n",
    "    LiH_energy_LBCS = cir.shadow_trace(LiH_hamiltonian, 1262, method=\"LBCS\")\n",
    "    LiH_energy_APS = cir.shadow_trace(LiH_hamiltonian, 1262, method=\"APS\")\n",
    "    LiH_energy_traditional = cir.expecval(LiH_hamiltonian, shots=2).numpy()[0]\n",
    "\n",
    "    estimator_CS.append(LiH_energy_CS) \n",
    "    estimator_LBCS.append(LiH_energy_LBCS) \n",
    "    estimator_APS.append(LiH_energy_APS) \n",
    "    estimator_traditional.append(LiH_energy_traditional) \n",
    "\n",
    "ave_LiH_energy_CS = np.mean(estimator_CS)\n",
    "ave_LiH_energy_LBCS = np.mean(estimator_LBCS)\n",
    "ave_LiH_energy_APS = np.mean(estimator_APS)\n",
    "ave_LiH_energy_traditional = np.mean(estimator_traditional)\n",
    "end = time.time() \n",
    "\n",
    "print(\"LiH ground state energy = \", LiH_energy)\n",
    "print(\"ave LiH ground state energy CS = \", ave_LiH_energy_CS)\n",
    "print(\"ave LiH ground state energy LBCS = \", ave_LiH_energy_LBCS)\n",
    "print(\"ave LiH ground state energy APS = \", ave_LiH_energy_APS)\n",
    "print('ave LiH ground state energy traditional = ', ave_LiH_energy_traditional)\n",
    "print('time = ', end-begin)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89469273",
   "metadata": {},
   "source": [
    "From the results, the mean values obtained by the classical-shadow-based algorithms are closer to the energy of the estimated ground state of $LiH$ by VQE than the item-by-item measurement method, and the errors of the algorithms are all within the accuracy of theoretical estimation in reference [3]. So, what about the sample variance of each algorithm?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "742f4f15",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LiH variance CS =  0.034596840359755784\n",
      "LiH variance LBCS =  0.016602696085670984\n",
      "LiH variance APS =  0.0016603026356630662\n",
      "LiH variance traditional =  0.13200055652163223\n"
     ]
    }
   ],
   "source": [
    "# Calculate the sample variance\n",
    "variance_CS = []\n",
    "variance_LBCS = []\n",
    "variance_APS = []\n",
    "variance_traditional = []\n",
    "\n",
    "for i in range(n):\n",
    "    variance_CS.append((estimator_CS[i] - ave_LiH_energy_CS) ** 2)\n",
    "    variance_LBCS.append((estimator_LBCS[i] - ave_LiH_energy_LBCS) ** 2)\n",
    "    variance_APS.append((estimator_APS[i] - ave_LiH_energy_APS) ** 2)\n",
    "    variance_traditional.append((estimator_traditional[i] - ave_LiH_energy_traditional) ** 2)\n",
    "\n",
    "var_CS = sum(variance_CS)/(n-1)\n",
    "var_LBCS = sum(variance_LBCS)/(n-1)\n",
    "var_APS = sum(variance_APS)/(n-1)\n",
    "var_traditional = sum(variance_traditional)/(n-1)\n",
    "\n",
    "print('LiH variance CS = ', var_CS)\n",
    "print('LiH variance LBCS = ', var_LBCS)\n",
    "print('LiH variance APS = ', var_APS)\n",
    "print('LiH variance traditional = ', var_traditional)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d40ef806",
   "metadata": {},
   "source": [
    "It can be seen that the APS algorithm has the lowest sample variance, followed by the LBCS algorithm, then the CS algorithm, and finally the item-by-item measurement method. Accordingly, we can find that after the increase in the number of terms of the Hamiltonian, the classical-shadow-based algorithm has higher accuracy and more stability at the same cost than the item-by-item measurement method. Among them, the APS algorithm is the most stable one.\n",
    "\n",
    "It is worth mentioning that for the classical shadow algorithm, the scene of 12 qubits still can not show its significant advantages compared with some existing algorithms. In large-scale systems with more qubits, its advantages can be better demonstrated [1]."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa9e6a26",
   "metadata": {},
   "source": [
    "## Conclusion"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91a60770",
   "metadata": {},
   "source": [
    "This tutorial discusses how to use the classical-shadow-based algorithms to estimate the observable's expectation value and shows how to use the `shadow` function in Paddle Quantum. It can be seen that the improved algorithm based on the classical shadow can get a good estimate of the observable's expectation value. Compared with the item-by-item measurement method, when the number of samples are the same, its estimated value is more accurate and the algorithm is more stable. In the large-scale qubit scenario, the number of samples required by the classical shadow method is only a constant increase in some tasks. So its role in the NISQ (noisy intermediate-scale quantum) era will continue to be explored."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5fcdcad3",
   "metadata": {},
   "source": [
    "_______\n",
    "\n",
    "## References"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "947b52db",
   "metadata": {},
   "source": [
    "[1] Huang, Hsin-yuan, R. Kueng and J. Preskill. “Predicting many properties of a quantum system from very few measurements.” [Nature Physics (2020): 1-8.](https://www.nature.com/articles/s41567-020-0932-7?proof=t)\n",
    "\n",
    "[2] Hadfield, Charles, et al. \"Measurements of quantum hamiltonians with locally-biased classical shadows.\" [arXiv preprint arXiv:2006.15788 (2020).](https://arxiv.org/abs/2006.15788)\n",
    "\n",
    "[3] Hadfield, Charles. \"Adaptive Pauli Shadows for Energy Estimation.\" [arXiv preprint arXiv:2105.12207 (2021).](https://arxiv.org/abs/2105.12207)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}