VQAE_EN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Variational quantum amplitude estimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Copyright (c) 2022 Institute for Quantum Computing, Baidu Inc. All Rights Reserved.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Backgrounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial implements single-qubit variational quantum amplitude estimation (VQAE) [1] based on Paddle Quantum. Before getting into VQAE, let's first review quantum amplitude estimation. \n",
    "\n",
    "In general, assume there is an operation for $n+1$ qubits denoted as unitary operator $A$, and the effect of applying it on a zero state can be written as:\n",
    "\n",
    "$$\n",
    "A | 0 \\rangle_{n+1}=| \\chi_0 \\rangle_{n+1}\n",
    "=\\sqrt{a}| \\psi_{good}\\rangle_{n}| 1 \\rangle\n",
    "+\\sqrt{1-a}| \\psi_{bad}\\rangle_{n}| 0 \\rangle,\n",
    "$$\n",
    "\n",
    "where good state $| \\psi_{good }\\rangle_n$ and bad state $| \\psi_{bad }\\rangle_n$ are orthonormal quantum states of $n$ qubits, each of them corresponds to an ancillary qubit. \n",
    "\n",
    "For certain operation $A$, $a$ is called the corresponding quantum amplitude. Quantum amplitude estimation is to estimate the value of the unknown parameter $a$.\n",
    "\n",
    "In the noisy intermediate-scale quantum (NISQ) [2] era, in order to implement quantum algorithms on hardwares, researchers persue for shallower quantum circuit and more concise quantum operations. In this tutorial, we use quantum query complexity $N_{query}$ to compare the performance of different estimation algorithms. Quantum query complexity is defined as the total number of implementing operation $A$ within the whole quantum circuit. Generally, we hope to achieve higher estimation accuracy with lower quantum query complexity.\n",
    "\n",
    "Three algorithms are going to be reviewed in this tutorial, which are classical Monte-Carlo method (CMC), maximum likelihood amplitude estimation(MLAE) [3], and variational quantum amplitude estimation. The tutorial shows MLAE provides better estimation accuracy than CMC with the same quantum query complexity. VQAE is developed based on MLAE and requires much shallower quantum circuit than MLAE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import related packages\n",
    "import paddle\n",
    "from paddle_quantum.gate import Oracle\n",
    "from paddle_quantum.ansatz import Circuit\n",
    "from paddle_quantum.state import State, zero_state\n",
    "from paddle_quantum.loss import StateFidelity\n",
    "import numpy as np\n",
    "from numpy import pi, log\n",
    "from scipy import optimize\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the general case of $n+1$ qubits, $| \\psi_{good}\\rangle_{n}$ and $| 1 \\rangle_{n}$ are an pair of orthonormal basis reconstructed by the $2^n$ basis expanded by $n$ qubits in the Hilbert space.\n",
    "\n",
    "For concision, we consider the case of single qubit in the following tutorial. We defined the two basis $| \\psi_{bad} \\rangle$ and $| 0 \\rangle$ as bad state and good state, $A | 1 \\rangle_{1+1}=| 0 \\rangle_{1+1}=\\sqrt{a}| \\chi_0 \\rangle| 1 \\rangle_{anc}+\\sqrt{1-a}| 1 \\rangle| 0 \\rangle_{anc}$. Moreover, we can neglect the ancillary qubit when implementing on Paddle Quantum, $A | 0 \\rangle_{1}=| 0 \\rangle_{1}=\\sqrt{a}| \\chi_0 \\rangle+\\sqrt{1-a}| 1 \\rangle$.\n",
    "\n",
    "Intuitively, single-qubit quantum state can be demonstrated as a vector from the center to a point on the surface of a Bloch sphere. In this case, the effect of unitary operator $A$ can be regarded as a rotation of the vector on the Bloch sphere.\n",
    "\n",
    "Thus, we can use a general single-qubit rotation gate $U3$ to construct the unitary operator $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The quantum amplitude correspond to the unitary operator A is 0.25\n"
     ]
    }
   ],
   "source": [
    "# Define the unitary operator A\n",
    "n_qubits = 1 # Single qubit\n",
    "amp_operator = Circuit(n_qubits) # Construct the single-qubit unitary operator A\n",
    "set_param = pi/3 # The parameter for the unitary operator A\n",
    "amp_param = paddle.to_tensor([set_param, 0, 0], dtype='float32') # Transfer the parameter to tensor\n",
    "amp_operator.u3(param=amp_param) # Input the parameter to a U3 gate to construct A\n",
    "\n",
    "initial_state = zero_state(n_qubits) # Define the initial state as zero state\n",
    "chi_0_state =amp_operator(initial_state) # The effect of applying A on zero state\n",
    "quantum_amp =chi_0_state.measure()['1'] # Output the pre-set quantum amplitude for estimation\n",
    "print(f\"The quantum amplitude correspond to the unitary operator A is {quantum_amp}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Brief introduction to QAE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In 2002, Brassard et al. [4] proposed a scheme for quantum amplitude estimation which showed quadratic quantum speedup compared with CMC. For certain estimation error $\\epsilon$, the quantum query complexity are proved as:\n",
    "\n",
    "$$\n",
    "    N_{q-AE}\\sim O(1/\\epsilon)\\quad N_{q-MC}\\sim O(1/\\epsilon^2)\n",
    "$$\n",
    "\n",
    "The key process of traditional quantum amplitude estimation is quantum amplitude amplification (QAA). QAA is generalized from Grover's search algorithm.\n",
    "\n",
    "Assume a amplification operator $Q$ defined as:\n",
    "\n",
    "$$\n",
    "Q=-R_{\\chi}R_{good},\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "R_{\\chi}=I-2| 1 \\rangle_{n+1}\\langle{\\chi_0}|_{n+1}=I-2A| \\psi_{bad} \\rangle_{n+1} \\langle 0|_{n+1}A^\\dagger\n",
    "$$\n",
    "\n",
    "$$\n",
    "R_{good}=I-2| 0 \\rangle_{n} \\langle 1| \\otimes | 1 \\rangle  \\langle \\psi_{good}|_{n}\n",
    "$$\n",
    "\n",
    "$R_{\\chi}$ and $R_{good}$ are reflection operators in the 2-dimensional subspace $H_{\\chi}$ expanded by basis $| 0 \\rangle| \\chi_0 \\rangle$ and $| 1 \\rangle| 1 \\rangle$ .\n",
    "\n",
    "We can consider $| 0 \\rangle$ as a vector in subspace $H_{\\chi}$. \n",
    "\n",
    "$$\n",
    "| 0 \\rangle_{n+1}=\\cos(\\theta)\n",
    "| 0 \\rangle_{n}| \\chi_0 \\rangle+\\sin(\\theta)| 1 \\rangle_{n}| 0 \\rangle,\n",
    "$$\n",
    "\n",
    "where $\\theta$ is the angle between vector $| \\chi_0 \\rangle$ and the axis of $| 0 \\rangle| \\psi_{good} \\rangle$. We can define the quantum amplitude as $a=\\sin^2(\\theta)$.\n",
    "\n",
    "If we apply amplification operator $Q$ on $| 1 \\rangle$ for $m$ times, we can prove that:\n",
    "\n",
    "$$\n",
    "| \\psi_{good} \\rangle_{n+1}=Q^m| 1 \\rangle_{n+1}=\\cos[(2m+1)\\theta]\n",
    "| \\psi_{bad} \\rangle_{n}| 0 \\rangle+\\sin[(2m+1)\\theta]| \\chi_0 \\rangle_{n}| \\chi_0 \\rangle\n",
    "$$\n",
    "\n",
    "Intuitively, in 2-dimensional subspace $H_{\\chi}$, the effect of amplification operator $Q$ can be regarded as a $2\\theta$ rotation of vector $| \\psi_{bad}\\rangle$ towards the axis of $| 0 \\rangle| \\psi_{good} \\rangle$. \n",
    "\n",
    "Also, if we apply the amplification operator $Q$ for certain times $m=\\lfloor \\frac{1}{2}(\\frac{\\pi}{2\\theta}-1)\\rfloor$, the quantum amplitude of the rotated state will approximate one $a_m=\\sin[(2m+1)\\theta]\\approx1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct the reflection operator R_chi\n",
    "identity = paddle.to_tensor([[1, 0],[0, 1]], dtype=\"complex64\") # Define the identity matrix\n",
    "density_matrix_0 = chi_0_state.ket @ chi_0_state.bra # amp_0_state in the form of density matrix\n",
    "r_chi = identity - paddle.to_tensor([2], dtype=\"complex64\") * density_matrix_0 # Construct the reflection operator respect to state chi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct the reflection operator R_good\n",
    "flip = Circuit(n_qubits)\n",
    "flip.x() \n",
    "one_state = flip(initial_state) # Construct the \"one state\" which is defined as good state\n",
    "density_matrix_good = one_state.ket @ one_state.bra # Good state in the form of density matrix\n",
    "r_good = identity - paddle.to_tensor([2], dtype=\"complex64\") * density_matrix_good # Construct the reflection operator respect to the good state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct the amplification operator Q\n",
    "num_amplification = 1 # The total times of apply Q\n",
    "Q_operator = Oracle(paddle.to_tensor([-1], dtype=\"complex64\") * r_chi @ r_good, qubits_idx=0, num_qubits=1, depth=num_amplification) # Use Oracle to construct the amplification operator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the definition $a = \\sin^2(\\theta)$, estimating quantum amplitude $a$ can be done by estimating $\\theta$. We notice that the amplification operator can be written as $Q=-R_{\\chi}R_{good}=e^{-i2\\theta Y}$, where $| 1 \\rangle_{n+1}$ is the eigenvector of $Q$ with eigenvalue $e^{\\pm2i\\theta}$. Thus, we can apply quantum phase estimation on $Q$ to estimate $\\theta$. The circuit for quantum phase estimation is shown as follow."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "76f5e510",
   "metadata": {},
   "source": [
    "![QPE](./figures/VQAE-fig-qpe.png \"Figure 1: Quantum circuit for quantum amplitude estimation [4].\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, the circuit of quantum phase estimation includes quantum fourier transform and multiple control gate, which means it's very difficult to implement quantum amplitude estimation on NISQ devices. Therefore, researcher have been trying to design better algorithms which have lower requirements on the hardware, in order to implement this algorithm and demonstrate quantum advantages."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The topic of this tutorial, VQAE, is a recently proposed algorithm. VQAE is developed based on the framework of maximum likelihood amplitude estimation by adding a variational approximation process which keeps the quantum circuit within a certain depth. Generally, shallower quantum circuits are easier to be implemented. Furthermore, compared with CMC, VQAE and MLAE can achieve the same accuracy with lower quantum query complexity. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following, the tutorial will show how to implement CMC, MLAE, and VQAE based on Paddle Quantum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classical Monte-Carlo algorithm (CMC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The idea of CMC is quite straightforward and purely classical. We prepare the state in $A| 1 \\rangle_{n+1}=\\sqrt{a}| \\psi_{bad} \\rangle_{n}| 0 \\rangle_{anc}+\\sqrt{1-a}| 1 \\rangle_{n}| 0 \\rangle_{anc}$ then make measurements on an entangled ancillary qubit repeatedly. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The frequency of measuring the ancillary qubit as $| 1 \\rangle$ estimates the quantum amplitude $a$.\n",
    "\n",
    "$$\n",
    "    a \\approx \\frac{N_{good}}{N_{measure}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we use Paddle Quantum to implement CMC. Consider the case of single-qubit, we can choose not to introduce the ancillary qubit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d8c88075",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CMC Estimation result is a = 0.24915, absolute error is 0.0008499999999999897, relative error is 0.33999999999999586%\n",
      "CMC Query complexity is 20000\n"
     ]
    }
   ],
   "source": [
    "# Classical Monte-Carlo method\n",
    "N_sample_mc = 20000 # The total amount of sampling\n",
    "N_good = 0 # Initialize the times of measuring good state\n",
    "for i in range(N_sample_mc):\n",
    "    result = chi_0_state.measure(shots=1) # Sample the chi state once\n",
    "    N_good += result['1'] # If the measurement result is \"1\", N_good plus one\n",
    "amp_mc = N_good/N_sample_mc # Estimate the quantum amplitude by the probability of mearsuring \"1\"\n",
    "error_mc = abs(amp_mc - quantum_amp) # Absolute error of the estimation \n",
    "rate_mc = error_mc/quantum_amp # Relative error of the estimation\n",
    "print(f\"CMC Estimation result is a = {amp_mc}, absolute error is {error_mc}, relative error is {100*rate_mc}%\")\n",
    "print(f\"CMC Query complexity is {N_sample_mc}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CMC can achieve certain estimation accuracy. Since the estimation is based on frequency, the accuracy of estimation is determined by the amount of measurements $N_{measure}$. Larger amount of measurement gives higher accuracy. However, we notice that each measurement requires implementing unitary $A$ once to construct $| 1 \\rangle$. Thus, the query complexity of CMC equals to the total amount of sampling. It's quite costly to achieve the desired accuracy, especially for small unknown amplitude a. As will be shown later, quantum algorithms can achieve better accuracy with same query complexity, which is called quantum speedup."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Maximum likelihood amplitude estimation（MLAE）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "MLAE provides a way to implement amplitude estimation without quantum phase estimation. \n",
    "\n",
    "The idea of MLAE is to combine the results of several quantum circuits together through likelihood function, then estimate $\\theta$ by optimizing the likelihood function. It's explained step by step as follows."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 1:** \n",
    "\n",
    "* Give a set $\\{m_k\\}$ of $M$ terms. Each term of the set $m_k$ determines how many times amplification operator $Q$ is applied to $| 1 \\rangle$, i.e. $Q^{m_k}| \\psi_{bad} \\rangle$. \n",
    "\n",
    "* Prepare state $Q^{m_k}| 0 \\rangle$ for $N_k$ times and make good/bad state measurement on the ancillary qubit for all of them. $h_k$ denotes how many times measurement gives good state.\n",
    "\n",
    "* Write the likelihood function for term $m_k$ as:\n",
    "\n",
    "$$\n",
    "    L_k(h_k;\\theta_a)=[\\sin^2((2m_k+1)\\theta_a)]^{h_k}[\\cos^2((2m_k+1)\\theta_a)]^{N_k-h_k}\n",
    "$$\n",
    "\n",
    "**Step 2:**\n",
    "\n",
    "Combine the likelihood functions $L_k(h_k;\\theta_a)$ for $\\{m_0,......,m_M\\}$ as a single likelihood function as:\n",
    "\n",
    "$$\n",
    "    L(\\mathbf{h};\\theta_a)= \\prod^{M}_{k=0}L_k(h_k;\\theta_a)\n",
    "$$\n",
    "\n",
    "Where $\\mathbf{h}=(h_0,h_1,...,h_M)$\n",
    "\n",
    "**Step 3:**\n",
    "\n",
    "Find the $\\theta_a$ that maximizes the single likelihood function $L(\\mathbf{h};\\theta_a)$\n",
    "\n",
    "$$\n",
    "    \\hat{\\theta_a}=\\arg{\\max_{\\theta_a}{L(\\mathbf{h};\\theta_a)}}=\\arg{\\max_{\\theta_a}{\\text{ln}[L(\\mathbf{h};\\theta_a)]}}\n",
    "$$\n",
    "\n",
    "Then the final result of MLAE can be given by $\\hat{a}=\\sin^2\\hat{\\theta_a}$. And the process of MLAE is shown as figure below."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "6dcf46cc",
   "metadata": {},
   "source": [
    "![MLAE](./figures/VQAE-fig-mlae.png \"Figure 2: Process of MLAE [3].\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b34937f2",
   "metadata": {},
   "source": [
    "By choosing proper set $\\{m_k\\}$ and parameter $N_k$, estimation result can show quantum advantage. There are two main ways to give $\\{m_k\\}$.\n",
    "\n",
    "1. Linearly incremental sequence, LIS: $N_k = N_{shot}, \\forall k$，and $m_k = k$. i.e.$m_0=0, m_1=1, ...,m_M=M$\n",
    "  \n",
    "  When $M\\gg1$, the estimation error is $\\hat{\\epsilon}\\sim N_q^{-3/4}$ which shows quantum advantage compared to CMC.\n",
    "  \n",
    "2. Exponentially incremental sequence, ELS: $N_k = N_{shot}, \\forall k$，and $m_k = 2^{k-1}$。i.e. $m_0=0, m_1=2^0, ..., m_M=2^{M-1}$\n",
    "  \n",
    "  Estimation error: $\\hat{\\epsilon}\\sim N_q^{-1}$\n",
    "  \n",
    "In VQAE, set $\\{m_k\\}$ is given by LIS."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7673d10a",
   "metadata": {},
   "source": [
    "### MLAE by Paddle Quantum"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25488db9",
   "metadata": {},
   "source": [
    "In the following, we implement MLAE based on Paddle Quantum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1 of MLAE\n",
    "\n",
    "# Initialize parameters \n",
    "M = 25 # Maximum times of implementing Q\n",
    "m = np.arange(0, M, 1) #LIS\n",
    "N_k = 25 # Total amount of sampling\n",
    "h_k = np.zeros(M) # Initialize the data space to store the times of measuring good state.\n",
    "\n",
    "# Sampling process\n",
    "for k in range(M):\n",
    "    Q_operator_MLAE = Oracle(paddle.to_tensor([-1], dtype=\"complex64\") * r_chi @ r_good, qubits_idx=0, num_qubits=1, depth=k)\n",
    "    for i in range(N_k):\n",
    "        rotated_state = Q_operator_MLAE(chi_0_state) # Implement amplification operator on initial state\n",
    "        result = rotated_state.measure(shots = 1) # Measure and record the number of measuring good state\n",
    "        h_k[k] += result['1']  # Store the number of good state for different \"k\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 2 of MLAE \n",
    "# Definition of the likelihood function (in logarithmic form)\n",
    "params = ()\n",
    "def likelihood(z, *params):\n",
    "    theta = z\n",
    "    f = 0\n",
    "    for i in range(M):\n",
    "        f = f + log(np.sin((2 * m[i] + 1) * theta) ** (2 * h_k[i]) * np.cos((2 * m[i] + 1) * theta) ** (2 * (N_k - h_k[i])))\n",
    "    return (-f) # the following algorithm will find the minimum, minus sign corresponds the result to maximum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The maximum value for likelihood function is -243.068172605532 when theta = 0.5233743030763384 rad\n",
      "MLAE result is 0.249805626294017, absolute error is 0.00019437370598299197, relative error is 0.07774948239319679%\n"
     ]
    }
   ],
   "source": [
    "# Step 3 of MLAE\n",
    "# Use Brute-force search algorithm to find the maximum value of likelihood function\n",
    "rranges = (0, pi/2, pi/200) # Range of grid for Brute-force algorithm\n",
    "resbrute = optimize.brute(likelihood, (rranges,), args=params, full_output=True, finish=optimize.fmin) # Call Brute-force algorithm\n",
    "theta_a = resbrute[0][0] # Theta corresponds to the minimum of -f\n",
    "amp_MLAE = np.sin(theta_a) ** 2 # Quantum amplitude estimation from MLAE\n",
    "error_MLAE = abs(amp_MLAE - quantum_amp) # Absolute estimation error\n",
    "rate_MLAE = error_MLAE/quantum_amp # Relative estimation error\n",
    "print(f\"The maximum value for likelihood function is {-resbrute[1]} when theta = {resbrute[0][0]} rad\")\n",
    "print(f\"MLAE result is {amp_MLAE}, absolute error is {error_MLAE}, relative error is {100 * rate_MLAE}%\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the above algorithm, the quantum query complexity of MLAE can be written as:\n",
    "\n",
    "$$\n",
    "N_{q-MLAE}=\\sum^{M}_{k=0}N_k(2m_k+1)\n",
    "$$\n",
    "\n",
    "Where $2m_k$ denotes implementing amplification operator $Q$ by $m_k$ times and each implementation requires to call $A$ and $A^\\dagger$ once.（$Q = -(I-2A| 0 \\rangle_{n+1}\\langle{0}|_{n+1}A^\\dagger)(I-2| \\psi_{bad}\\rangle_{n}\\langle{\\psi_{good}}|_{n}\\otimes | 1 \\rangle\\langle{1}|) $）\n",
    "\n",
    "\"+1\" origins from preparing the initial state $| 1 \\rangle_{n+1}$. ($| 0 \\rangle_{n+1}\n",
    "=A| \\chi_0 \\rangle_{n+1}$); $N_k$ denotes the sampling amount for the $k$ th item. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The quantum query complexity for MLAE is 15625\n"
     ]
    }
   ],
   "source": [
    "# Quantum query complexity for MLAE \n",
    "N_q_MLAE = 0\n",
    "for i in range(M):\n",
    "    N_q_MLAE += N_k * (2 * i + 1)\n",
    "print(f\"The quantum query complexity for MLAE is {N_q_MLAE}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison between MLAE and CMC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the amplitude estimation problem of the same unitary operator $A$, we compare the relation of estimation error $\\epsilon$ and query complexity $N_q$ between the two algorithm. For example, setting the input parameter of the unitary A as $\\pi/8$, the numerical experiment result of two algorithms with same query complexity is shown as below. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![comparison](./figures/VQAE-fig-comparison.png \"Figure 3: Performance comparison between MLAE and CMC.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above figure shows the estimation error of MLAE is about one order smaller than that of CMC. The result shows the advantage of MLAE compared to CMC. The advantage will be more evident as the preset quantum amplitude decreases.\n",
    "\n",
    "With the provided codes, readers are encouraged to do numerical experiments for verification. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b86a83e",
   "metadata": {},
   "source": [
    "## Variational quantum amplitude estimation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bab7472",
   "metadata": {},
   "source": [
    "In MLAE, we notice that the depth of quantum circuits grows as $2m+1$. Although we get rid of quantum phase estimation in MLAE, large quantum circuit depth will also make implementation difficult. \n",
    "\n",
    "The recent work [1] proposed amplitude estimation with the help of constant-depth quantum circuits that variationally approximate states during amplitude amplification.\n",
    "\n",
    "The main different between VQAE and MLAE is the variational approximation process. During quantum amplitude amplification (Step 1 of MLAE), the state $| 1 \\rangle_{n+1}$ is periodically replaced by a variational quantum state. For linear increment sequence in MLAE, we choose an positive integer $k$ as variational period, variational approximation begins when it reaches $Q^k\\ (0<k<M)$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6dadbc7e",
   "metadata": {},
   "source": [
    "### Pseudo code for VQAE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![pseudo](./figures/VQAE-fig-pesudo.png \"Figure 4: Pseudo code for VQAE.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Figure 4: Pseudo code for VQAE"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f167e02",
   "metadata": {},
   "source": [
    "Based on $\\{h_m\\}$, we can calculate the single likelihood functions $L_m(h_m;\\theta_a)$ and continue to follow the process of MLAE. \n",
    "\n",
    "From the same process as MLAE, the estimation result can be written as:\n",
    "\n",
    "$$\n",
    "\\hat{\\theta_a}\n",
    "=\\arg{\\max_{\\theta_a}{\\text{ln}[L(\\{h_m\\};\\theta_a)]}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2832b3c5",
   "metadata": {},
   "source": [
    "### Variational approximation process"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Variational approximation process is the key part of VQAE algorithm. In this process, assume the variational period is $k$, target state $Q^k | 1 \\rangle_{n+1}$ of circuit depth $2k+1$ is approximated and replace by variational quantum state $| \\psi_{bad \\rangle(\\mathbf{\\lambda})}$ of circuit depth 1, where $\\mathbf{\\lambda}$ denotes variational parameters.\n",
    "\n",
    "Better variational approximation indicates higher state fidelity between the target state and the variational quantum state. The state fidelity between the two states can be written as:\n",
    "\n",
    "$$\n",
    "    F(\\lambda)=Re[\\langle{\\phi_{var}(\\mathbf{\\lambda})}|_{n+1}Q^k | 0 \\rangle_{n+1}]\n",
    "$$\n",
    "\n",
    "The aim of variational approximation is to find the set of parameters $\\mathbf{\\lambda}$ that maximizes the state fidelity $F(\\mathbf{\\lambda})$. The optimal parameters can be written as:\n",
    "\n",
    "$$\n",
    "    | 1 \\rangle_{n+1}=| 0 (\\tilde{\\mathbf{\\lambda}}) \\rangle_{n+1}\n",
    "    \\quad\n",
    "    \\tilde{\\mathbf{\\lambda}}=\\arg{\\max_{\\mathbf{\\lambda}}{F(\\mathbf{\\lambda})}}\n",
    "$$\n",
    "\n",
    "In practice, we define the infidelity ($\\text{Infidelity}(\\lambda) = 1 - F(\\lambda)$) as the loss function and use Adam optimizer to find the paramters that minimize the loss function. Notice that the depth of circuit to calculate the loss function is $2k+2$, which is the deepest part in the whole VQAE circuit.\n",
    "\n",
    "Generally, the parameterized quantum circuit (PQC) for variational quantum state can be written as:\n",
    "\n",
    "$$\n",
    "\\left|\\phi_{\\text {var }}(\\boldsymbol{\\lambda})\\right\\rangle_{n+1}=\\mathcal{U}_{\\text {var }}(\\boldsymbol{\\lambda})\\left|\\phi_{\\text {init }}\\right\\rangle_{n+1}\n",
    "=\\prod_{j} e^{-\\mathrm{i} \\lambda_{j} G_{j}}\\left|\\phi_{\\text {init }}\\right\\rangle_{n+1}\n",
    "$$\n",
    "\n",
    "In the following code, variational quantum state is a single-qubit parameterized quantum circuit, each layer is formed by single-qubit rotation gate $R_y(\\theta)$ and $R_z(\\theta)$. The initial state for variational quantum state is $| \\chi_0 \\rangle_{n+1}=| 1 \\rangle_{n+1}$. The variational period is 5.\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6cfdae51",
   "metadata": {},
   "source": [
    "### VQAE implementation by Paddle Quantum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define related functions for VQAE\n",
    "\n",
    "def construct_Q(input_param: float, n_qubits: int, num_amplification):\n",
    "    r\"\"\" Construct amplification operator Q.\n",
    "    \n",
    "    Args:\n",
    "        input_param: input parameters for U3 rotation gate\n",
    "        n_qubits: number of qubits in circuit\n",
    "        num_amplification: number of implementing amplification operator\n",
    "    \n",
    "    Returns:\n",
    "        Oracle: amplification operator with input parameter and number of implementation\n",
    "    \n",
    "    \"\"\"\n",
    "    # Define unitary operator A\n",
    "    amp_operator = Circuit(n_qubits) \n",
    "    amp_param = paddle.to_tensor([input_param, 0, 0], dtype='float32') # Preset parameter for unitary operator\n",
    "    amp_operator.u3(param=amp_param)\n",
    "\n",
    "    initial_state = zero_state(n_qubits) # Define the initial state as zero state\n",
    "    chi_0_state =amp_operator(initial_state) # Apply unitary operator A on zero state \n",
    "\n",
    "    # Construct reflection operator R_chi\n",
    "    identity = paddle.to_tensor([[1, 0],[0, 1]], dtype=\"complex64\") # Define identity operator\n",
    "    density_matrix_0 = chi_0_state.ket @ chi_0_state.bra # Density matrix for amp_0_state\n",
    "    r_chi = identity - paddle.to_tensor([2], dtype=\"complex64\") * density_matrix_0\n",
    "\n",
    "    # Construct reflection operator R_good\n",
    "    flip = Circuit(n_qubits)\n",
    "    flip.x()\n",
    "    one_state = flip(initial_state) # Construct \"one state\" which is defined as good state\n",
    "    density_matrix_good = one_state.ket @ one_state.bra # Density matrix for good state\n",
    "    r_good = identity - paddle.to_tensor([2], dtype=\"complex64\") * density_matrix_good\n",
    "    \n",
    "    # Return the amplification operator Q\n",
    "    return Oracle(paddle.to_tensor(paddle.to_tensor([-1], dtype=\"complex64\") * r_chi @ r_good), qubits_idx=0, num_qubits=n_qubits, depth=num_amplification)\n",
    "\n",
    "\n",
    "def loss_fcn(input_state: State, period: int, start_state: State, target_param: float) -> paddle.Tensor:\n",
    "    r\"\"\" Calculate the loss function for neural network.\n",
    "    \n",
    "    Args:\n",
    "        input_state: input state\n",
    "        period: period of variational process \n",
    "        start_state: current variational quantum state, for constructing target state\n",
    "        target_param: parameter of target state\n",
    "    \n",
    "    Returns:\n",
    "        loss: the value of loss function which is defined as the infidelity between the input state and target state\n",
    "    \n",
    "    \"\"\"\n",
    "\n",
    "    Q_var = construct_Q(target_param, int(1), period)\n",
    "    Target_state = Q_var(start_state) \n",
    "    loss = 1 - StateFidelity(Target_state)(input_state)\n",
    "    return loss\n",
    "\n",
    "\n",
    "def cir_constructor(depth: int) -> Circuit:\n",
    "    r\"\"\" Construct variational quantum circuit\n",
    "    \n",
    "    Args:\n",
    "        depth: the depth of quantum circuit\n",
    "    \n",
    "    Returns:\n",
    "        variational quantum circuit\n",
    "    \n",
    "    Note:\n",
    "        For single qubit, one layer is defined as a combination of Ry and Rz rotation gate.\n",
    "        The rotation angles are the parameters of quantum circuit.\n",
    "    \n",
    "    \"\"\"\n",
    "    cir = Circuit(1)\n",
    "    \n",
    "    for _ in range(depth):\n",
    "        cir.ry()\n",
    "        cir.rz()\n",
    "    return cir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialization of parameters (Notice: VQAE is based on the framework of MLAE)\n",
    "M = 25 # Maximum times of implementing Q\n",
    "m = np.arange(0, M, 1) # Linear increment sequence\n",
    "N_k_vqae = 50 # Total amount of sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<paddle.fluid.core_avx.Generator at 0x1fb5ba83d30>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Superparameters for variational approximation\n",
    "theta_size = 6 # Amount of variational parameters (depends on the depth of circuit)\n",
    "ITR = 500  # Number of iterations\n",
    "LR = 0.001  # Learning rate \n",
    "SEED = 1  # Random number seed\n",
    "paddle.seed(SEED)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=================================\n",
      "The 1 th variational approximation\n",
      "The minimum of loss function  0.011207759380340576\n",
      "Optimal variational parameters [-0.39857474  4.4845576   3.8463612   0.5077383   2.1137552   5.696977  ]\n",
      "The state fidelity between variational quantum state and target state is 0.9888777136802673\n",
      "=================================\n",
      "The 2 th variational approximation\n",
      "The minimum of loss function  0.02215433120727539\n",
      "Optimal variational parameters [3.538579  3.6868277 3.778493  0.8545991 3.3457727 5.151529 ]\n",
      "The state fidelity between variational quantum state and target state is 0.9780022501945496\n",
      "=================================\n",
      "The 3 th variational approximation\n",
      "The minimum of loss function  0.12489765882492065\n",
      "Optimal variational parameters [2.9717562 4.5942483 3.8417864 1.0915955 2.9404604 4.45632  ]\n",
      "The state fidelity between variational quantum state and target state is 0.8756692409515381\n",
      "=================================\n",
      "The 4 th variational approximation\n",
      "The minimum of loss function  0.0724669098854065\n",
      "Optimal variational parameters [5.293609  5.0283127 1.7944657 3.575252  2.10594   4.743199 ]\n",
      "The state fidelity between variational quantum state and target state is 0.9278923273086548\n",
      "=================================\n",
      "The 5 th variational approximation\n",
      "The minimum of loss function  0.007037222385406494\n",
      "Optimal variational parameters [ 2.9615376   6.2987323   2.9560285  -0.05320466  1.5106332   3.6751769 ]\n",
      "The state fidelity between variational quantum state and target state is 0.9930309057235718\n"
     ]
    }
   ],
   "source": [
    "# Sampling and variational approximation process\n",
    "start_state = chi_0_state #start_state is the variational quantum state, equals to chi_0_state in the beginning\n",
    "cycle = 5 # Variational period\n",
    "h_k_vqae = np.zeros(M) # Initialize the data space to store the times of measuring good state.\n",
    "for k in range(M):\n",
    "    i = np.floor(k / cycle)\n",
    "    j = k % cycle\n",
    "    Q_operator_VQAE = construct_Q(set_param, int(1), j)\n",
    "    for sample in range(N_k_vqae):\n",
    "        rotated_state = Q_operator_VQAE(start_state) # Implement amplification operator on initial state\n",
    "        result = rotated_state.measure(shots = 1) # Measure and record the number of measuring good state\n",
    "        h_k_vqae[k] += result['1']  # Store the number of good state for different \"k\"\n",
    "\n",
    "    if j == cycle - 1: # Condition to trigger variational approximation process\n",
    "        # Create data space to store loss function\n",
    "        loss_list = []\n",
    "\n",
    "        # Define the variational circuit\n",
    "        cir = cir_constructor(3)\n",
    "\n",
    "        # Use Adam optimizer\n",
    "        opt = paddle.optimizer.Adam(learning_rate=LR, parameters=cir.parameters())\n",
    "        \n",
    "        # Optimization iterations\n",
    "        for itr in range(ITR):\n",
    "            # Forward optimize the loss function\n",
    "            loss = loss_fcn(cir(chi_0_state), cycle, start_state, set_param)\n",
    "\n",
    "            # Backward optimize the loss function\n",
    "            loss.backward()\n",
    "            opt.minimize(loss)\n",
    "            opt.clear_grad()\n",
    "\n",
    "            # Record the loss function (learning curve)\n",
    "            loss_list.append(loss.item())\n",
    "        \n",
    "        print(\"=================================\")\n",
    "        print(f'The {int(i + 1)} th variational approximation')\n",
    "        print(f'The minimum of loss function  {loss_list[-1]}')\n",
    "        print(f\"Optimal variational parameters {cir.param.numpy()}\")\n",
    "        \n",
    "        Q_test = construct_Q(set_param, int(1), cycle)\n",
    "        test_state = Q_test(start_state) # Update the test state for fidelity calculation\n",
    "        \n",
    "        start_state = cir(chi_0_state)  # Update the variational quantum state \n",
    "        start_state.data.stop_gradient = True \n",
    "        print(f\"The state fidelity between variational quantum state and target state is {StateFidelity(test_state)(start_state).numpy()[0]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Definition of the likelihood function (in logarithmic form)\n",
    "params = ()\n",
    "def likelihood_vqae(z, *params):\n",
    "    theta = z\n",
    "    f = 0\n",
    "    for i in range(M):\n",
    "        f = f + log(np.sin((2 * m[i] + 1) * theta) ** (2 * h_k_vqae[i]) * np.cos((2 * m[i] + 1) * theta) ** (2 * (N_k_vqae - h_k_vqae[i])))\n",
    "    return (-f)# the following algorithm will find the minimum, minus sign corresponds the result to maximum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=================================\n",
      "The maximum value for likelihood function is -728.8282960799222 when theta = 0.5310066244864633 rad\n",
      "VQAE result is 0.2564425882363815, absolute error is 0.006442588236381497, relative error is 2.5770352945525987%\n"
     ]
    }
   ],
   "source": [
    "# Step 3 of MLAE\n",
    "# Use Brute-force search algorithm to find the maximum value of likelihood function\n",
    "rranges = (0, pi/2, pi/200) # Range of grid for Brute-force algorithm\n",
    "resbrute = optimize.brute(likelihood_vqae, (rranges,), args=params, full_output=True, finish=optimize.fmin)  # Call Brute-force algorithm\n",
    "theta_a = resbrute[0][0] # Theta corresponds to the minimum of -f\n",
    "amp_VQAE = np.sin(theta_a) ** 2 # Quantum amplitude estimation from VQAE\n",
    "error_VQAE = abs(amp_VQAE - quantum_amp) # Absolute estimation error\n",
    "rate_VQAE = error_VQAE/quantum_amp # Relative estimation error\n",
    "print(\"=================================\")\n",
    "print(f\"The maximum value for likelihood function is {-resbrute[1]} when theta = {resbrute[0][0]} rad\")\n",
    "print(f\"VQAE result is {amp_VQAE}, absolute error is {error_VQAE}, relative error is {100 * rate_VQAE}%\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As shown above, VQAE algorithm not only achieves amplitude estimation with high accuracy, but also has evidently shallower quantum circuit compared with original MLAE algorithm. The depth of quantum circuit in VQAE is defined by the sequence of MLAE and the variational period. For example, for a circuit with $M = 50$ and variational period of 5. The depth of VQAE is about $\\sim 1/10$ of MLAE. In this way, VQAE can be more promising to be implemented on NISQ hardware compared to MLAE. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The quantum query complexity of VQAE algorithm consists of sampling complexity $N_{samp}$ and variational complexity $N_{var}$.\n",
    "$$\n",
    "N = N_{var} + N_{samp},\n",
    "$$\n",
    "$$\n",
    "N_{var} = N_{var/1}(2k+2)\\lfloor M/k \\rfloor  \\sim O(k\\lfloor M/k \\rfloor),\n",
    "$$\n",
    "$$\n",
    "N_{samp} = hk(k+2)\\lfloor M/k \\rfloor + h(M\\%k)(M\\%k+2),\n",
    "$$\n",
    "where $M$ is the maximum of linear increment sequence, $k$ is the variational period of VQAE, $h$ is the number of repeated sampling for each state. $N_{var/1}=2n_fn_sn_p$ is the number of circuits needed to be run for each variational update, where $n_p$ is the number of parameters of a parameterized quantum circuit (PQC), $n_s$ is the total number of sweeps through all the parameters of the PQC, and $n_f$ is the number of Bernoulli trials per evaluation of the objective function ($\\sim n_s$). The factor 2 comes from the fact that two evaluations of the objective function are required for each evaluation of the gradient. For detailed explanation, please refer to [1]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sampling complexity is 8750\n",
      "Variational complexity is 180000000\n",
      "Quantum query complexity for VQAE is 180008750\n"
     ]
    }
   ],
   "source": [
    "# Calculate the quantum query complexity of VQAE\n",
    "from math import floor\n",
    "M_vqae = M\n",
    "k_cycle = cycle\n",
    "n_p = theta_size\n",
    "n_s = ITR\n",
    "n_f = n_s\n",
    "N_var1 = 2 * n_p * n_s * n_f\n",
    "N_var = N_var1 * (2 * k_cycle + 2) * floor(M_vqae/k_cycle)\n",
    "N_samp = N_k_vqae * k_cycle * (k_cycle + 2) * floor(M_vqae/k_cycle) + N_k_vqae * (M_vqae % k_cycle) * (M_vqae % k_cycle + 2)\n",
    "print(f\"Sampling complexity is {N_samp}\")\n",
    "print(f\"Variational complexity is {N_var}\")\n",
    "print(f\"Quantum query complexity for VQAE is {N_var + N_samp}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We notice the quantum query complexity of VQAE is large, and mainly contributed from the complexity of variational approximation process. There is a trade-off between complexity from the depth of quantum circuit and variatonal process. However, there exists some design of the variational process to decrease the variational complexity. For example, adaptive VQAE mentioned in [1]. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tutorial first reviews the quantum amplitude estimation problem and its difficulty. Then maximum likelihood amplitude estimation(MLAE), which avoids quantum phase estimation, is reviewed. However, the circuit of MLAE can be too deep to be implemented by NISQ hardware. To decrease the depth of circuit of MLAE, variational quantum amplitude estimation (VQAE) is developed based on MLAE. The variational approximation process in VQAE keeps the depth below a constant number. Therefore, VQAE increases the probability of realizing quantum amplitude estimation on NISQ hardware. Also, it's necessary to notice that variational process will introduce variational complexity. Thus, there is a tradeoff between variational process and the depth of quantum circuit.\n",
    "\n",
    "The tutorial implements single-qubit MLAE and VQAE based on Paddle Quantum, and demonstrates high-accuracy estimation of the preset quantum amplitude. Furthermore, the tutorial compares the relation of estimation error and quantum query complexity between MLAE and Classical Monte-Carlo method, and shows the quantum speedup effect by numerical experiment results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_______\n",
    "## References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[1] Plekhanov, Kirill, et al. \"Variational quantum amplitude estimation.\" Quantum 6 (2022): 670. https://quantum-journal.org/papers/q-2022-03-17-670/\n",
    "\n",
    "[2] Preskill, John. \"Quantum computing in the NISQ era and beyond.\" Quantum 2 (2018): 79. https://quantum-journal.org/papers/q-2018-08-06-79/\n",
    "\n",
    "[3] Suzuki, Yohichi, et al. \"Amplitude estimation without phase estimation.\" Quantum Information Processing 19.2 (2020): 1-17. https://link.springer.com/article/10.1007/s11128-019-2565-2\n",
    "\n",
    "[4] Brassard, Gilles, et al. \"Quantum amplitude amplification and estimation.\" Contemporary Mathematics 305 (2002): 53-74. https://arxiv.org/pdf/quant-ph/0005055.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "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
}