tutorial/quantum_simulation/SSVQE_EN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Subspace-search Variational Quantum Eigensolver\n",
    "\n",
    "<em> Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved. </em>\n",
    "\n",
    "## Overview\n",
    "\n",
    "- In this tutorial, we will show how to train a quantum neural network (QNN) through Paddle Quantum to find the entire energy spectrum of a quantum system.\n",
    "\n",
    "- First, import the following packages."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "source": [
    "import numpy\n",
    "from numpy import pi as PI\n",
    "import paddle \n",
    "from paddle import matmul\n",
    "from paddle_quantum.circuit import UAnsatz\n",
    "from paddle_quantum.utils import random_pauli_str_generator, pauli_str_to_matrix, dagger"
   ],
   "outputs": [],
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:12:27.747028Z",
     "start_time": "2021-04-30T09:12:25.171248Z"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Background\n",
    "\n",
    "- Variational Quantum Eigensolver (VQE) [1-3] is one of the most promising applications for near-term quantum computing. One of the its powerful versions is SSVQE [4], which can be used to find the ground state and the **excited state** of a physical system's Hamiltonian. Mathematically, one can interpret it as solving the eigenvalues and eigenvectors of a Hermitian matrix. The set of eigenvalues of the Hamiltonian is called the energy spectrum.\n",
    "- Next, we will use a brief example to demonstrate how to solve this problem by training a QNN, that is, to solve the energy spectrum of a given Hamiltonian $H$.\n",
    "\n",
    "## Hamiltonian \n",
    "\n",
    "- For a specific molecule that needs to be analyzed, we need its geometry, charge, and spin multiplicity to obtain the Hamiltonian (in Pauli products form) describing the system. Specifically, through our built-in quantum chemistry toolkit, fermionic-to-qubit mapping technology can be used to output the qubit Hamiltonian.\n",
    "- As a simple demonstration of SSVQE, we provide a random 2-qubit Hamiltonian."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "source": [
    "N = 2  # Number of qubits\n",
    "SEED = 14  # Fixed random seed"
   ],
   "outputs": [],
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:12:27.773417Z",
     "start_time": "2021-04-30T09:12:27.752568Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "source": [
    "# Generate random Hamiltonian represented by Pauli string\n",
    "numpy.random.seed(SEED)\n",
    "hamiltonian = random_pauli_str_generator(N, terms=10)\n",
    "print(\"Random Hamiltonian in Pauli string format = \\n\", hamiltonian)\n",
    "\n",
    "# Generate matrix representation of Hamiltonian\n",
    "H = pauli_str_to_matrix(hamiltonian, N)"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "Random Hamiltonian in Pauli string format = \n",
      " [[0.9152074787317819, 'x1,y0'], [-0.2717604556798945, 'z0'], [0.3628495008719168, 'x0'], [-0.5050129214094752, 'x1'], [-0.6971554357833791, 'y0,x1'], [0.8651151857574237, 'x0,y1'], [0.7409989105435002, 'y0'], [-0.39981603921243236, 'y0'], [0.06862640764702, 'z0'], [-0.7647553733438246, 'y1']]\n"
     ]
    }
   ],
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:12:27.804177Z",
     "start_time": "2021-04-30T09:12:27.779339Z"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Building a quantum neural network\n",
    "\n",
    "- To implement SSVQE, we first need to design a QNN $U(\\theta)$ (parameterized quantum circuit). In this tutorial, we provide a predefined universal quantum circuit template suitable for 2 qubits. Theoretically, this template has enough expressibility to simulate arbitrary 2-qubit unitary operation [5]. The specific implementation requires 3 $CNOT$ gates plus 15 single-qubit rotation gates $\\in \\{R_y, R_z\\}$.\n",
    "\n",
    "- One can randomly initialize the QNN parameters ${\\bf{\\vec{\\theta }}}$ containing 15 parameters."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "source": [
    "THETA_SIZE = 15 # The number of parameters in the quantum neural network\n",
    "\n",
    "def U_theta(theta, N):\n",
    "    \"\"\"\n",
    "    U_theta\n",
    "    \"\"\"\n",
    "    # Initialize the quantum neural network according to the number of qubits/network width\n",
    "    cir = UAnsatz(N)\n",
    "    \n",
    "    # Call the built-in quantum neural network template\n",
    "    cir.universal_2_qubit_gate(theta, [0, 1])\n",
    "\n",
    "    # Return the circuit of the quantum neural network\n",
    "    return cir"
   ],
   "outputs": [],
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:12:27.822250Z",
     "start_time": "2021-04-30T09:12:27.809696Z"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Training model and loss function\n",
    "\n",
    "- After setting up the Hamiltonian and the quantum neural network architecture, we will further define the parameters to be trained, the loss function and optimization methods. For a detailed inspection of the theory of SSVQE, please refer to the original paper [4].\n",
    "\n",
    "- By acting the quantum neural network $U(\\theta)$ on a set of orthogonal initial states (one can take the computational basis $\\{|00\\rangle, |01\\rangle, |10\\rangle, |11 \\rangle \\}$), we will get the output states $\\{\\left| {\\psi_1 \\left( {\\bf{\\theta }} \\right)} \\right\\rangle, \\left| {\\psi_2 \\left( {\\bf{\\theta }} \\right)} \\right\\rangle, \\left| {\\psi_3 \\left( {\\bf{\\theta }} \\right)} \\right\\rangle, \\left| {\\psi_4 \\left( {\\bf{\\theta }} \\right)} \\right\\rangle \\}$.\n",
    "\n",
    "- Further, the loss function in the SSVQE model generally consists of expectation value of each output quantum state $\\left| {\\psi_k \\left( {\\bf{\\theta }} \\right)} \\right\\rangle$ given the Hamiltonian $H$. More specifically, it's the weighted summation of the energy expectation value. In this example, the default weight vector is $\\vec{w} = [4, 3, 2, 1]$.\n",
    "\n",
    "- The loss function is defined as:\n",
    "\n",
    "$$\n",
    "\\mathcal{L}(\\boldsymbol{\\theta}) = \\sum_{k=1}^{2^n}w_k*\\left\\langle {\\psi_k \\left( {\\bf{\\theta }} \\right)} \\right|H\\left| {\\psi_k \\left( {\\bf{\\theta }} \\right)} \\right\\rangle. \\tag{1}\n",
    "$$"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "source": [
    "class Net(paddle.nn.Layer):\n",
    "    def __init__(self, shape, dtype='float64'):\n",
    "        super(Net, self).__init__()\n",
    "        \n",
    "        # Initialize the theta parameter list and fill the initial value with the uniform distribution of [0, 2*pi]\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, H, N):\n",
    "        \n",
    "        # Build quantum neural network\n",
    "        cir = U_theta(self.theta, N)\n",
    "        U = cir.U\n",
    "        \n",
    "        # Calculate the loss function\n",
    "        loss_struct = paddle.real(matmul(matmul(dagger(U), H), U))\n",
    "\n",
    "        # Enter the computational basis to calculate the expected value \n",
    "        # which is equivalent to taking the diagonal element of U^dagger*H*U\n",
    "        loss_components = []\n",
    "        for i in range(len(loss_struct)):\n",
    "            loss_components.append(loss_struct[i][i])\n",
    "        \n",
    "        # Weighted summation of loss function\n",
    "        loss = 0\n",
    "        for i in range(len(loss_components)):\n",
    "            weight = 4 - i\n",
    "            loss += weight * loss_components[i]\n",
    "        \n",
    "        return loss, loss_components, cir"
   ],
   "outputs": [],
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:12:28.432737Z",
     "start_time": "2021-04-30T09:12:28.423798Z"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Hyper-parameters\n",
    "\n",
    "Before training the quantum neural network, we also need to set up several hyper-parameters, mainly the learning rate LR, the number of iterations ITR. Here we set the learning rate to be LR = 0.3 and the number of iterations ITR = 50. One can adjust these hyper-parameters accordingly and check how they influence the training performance."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "source": [
    "ITR = 100  # Set the total number of iterations of training\n",
    "LR = 0.3  # Set the learning rate"
   ],
   "outputs": [],
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:12:29.579180Z",
     "start_time": "2021-04-30T09:12:29.575632Z"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Training process\n",
    "\n",
    "- After setting all the parameters of SSVQE model, we need to convert all the data into Tensor in the PaddlePaddle, and then train the quantum neural network.\n",
    "- We use Adam Optimizer in training, and one can also call other optimizers provided in PaddlePaddle."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "source": [
    "paddle.seed(SEED)\n",
    "\n",
    "# We need to convert numpy.ndarray to Tensor supported in Paddle\n",
    "hamiltonian = paddle.to_tensor(H)\n",
    "\n",
    "# Determine the parameter dimension of the network\n",
    "net = Net(shape=[THETA_SIZE])\n",
    "\n",
    "# We use Adam optimizer for better performance\n",
    "# One can change it to SGD or RMSprop.\n",
    "opt = paddle.optimizer.Adam(learning_rate=LR, parameters=net.parameters())\n",
    "\n",
    "# Optimization loop\n",
    "for itr in range(1, ITR + 1):\n",
    "\n",
    "    # Forward propagation calculates the loss function and returns the estimated energy spectrum\n",
    "    loss, loss_components, cir = net(hamiltonian, N)\n",
    "\n",
    "    # Under the dynamic graph mechanism, use back propagation to minimize the loss function\n",
    "    loss.backward()\n",
    "    opt.minimize(loss)\n",
    "    opt.clear_grad()\n",
    "\n",
    "    # Print training results\n",
    "    if itr% 10 == 0:\n",
    "        print('iter:', itr,'loss:','%.4f'% loss.numpy()[0])\n",
    "    if itr == ITR:\n",
    "        print(\"\\nThe trained circuit:\")\n",
    "        print(cir)"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "iter: 10 loss: -4.5668\n",
      "iter: 20 loss: -5.3998\n",
      "iter: 30 loss: -5.6210\n",
      "iter: 40 loss: -5.8872\n",
      "iter: 50 loss: -5.9246\n",
      "iter: 60 loss: -5.9471\n",
      "iter: 70 loss: -5.9739\n",
      "iter: 80 loss: -5.9833\n",
      "iter: 90 loss: -5.9846\n",
      "iter: 100 loss: -5.9848\n",
      "\n",
      "The trained circuit:\n",
      "--U----x----Rz(-1.18)----*-----------------x----U--\n",
      "       |                 |                 |       \n",
      "--U----*----Ry(-0.03)----x----Ry(2.362)----*----U--\n",
      "                                                   \n"
     ]
    }
   ],
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:12:44.010556Z",
     "start_time": "2021-04-30T09:12:41.952650Z"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Benchmarking\n",
    "\n",
    "We have now completed the training of the quantum neural network, and we will verify the results by comparing them with theoretical values.\n",
    "- The theoretical Hamiltonian eigenvalues are solved by the linear algebra package in NumPy;\n",
    "- We compare the energy of each energy level obtained by training QNN with the theoretical value.\n",
    "- It can be seen that the training output is very close to the exact value."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "source": [
    "def output_ordinalvalue(num):\n",
    "    r\"\"\"\n",
    "    Convert to ordinal value\n",
    "\n",
    "    Args:\n",
    "        num (int): input number\n",
    "\n",
    "    Return:\n",
    "        (str): output ordinal value\n",
    "    \"\"\"\n",
    "    if num == 1:\n",
    "        return str(num) + \"st\"\n",
    "    elif num == 2:\n",
    "        return str(num) + \"nd\"\n",
    "    elif num == 3:\n",
    "        return str(num) + \"rd\"\n",
    "    else:\n",
    "        return str(num) + 'th'\n",
    "\n",
    "for i in range(len(loss_components)):\n",
    "    if i == 0:\n",
    "        print('The estimated ground state energy is: ', loss_components[i].numpy())\n",
    "        print('The theoretical ground state energy is: ', numpy.linalg.eigh(H)[0][i])\n",
    "    else:\n",
    "        print('The estimated {} excited state energy is: {}'.format(\n",
    "            output_ordinalvalue(i), loss_components[i].numpy())\n",
    "        )\n",
    "        print('The theoretical {} excited state energy is: {}'.format(\n",
    "            output_ordinalvalue(i), numpy.linalg.eigh(H)[0][i])\n",
    "        )"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "The estimated ground state energy is:  [-2.18762366]\n",
      "The theoretical ground state energy is:  -2.18790201165885\n",
      "The estimated 1st excited state energy is: [-0.13721024]\n",
      "The theoretical 1st excited state energy is: -0.13704127143749587\n",
      "The estimated 2nd excited state energy is: [0.85251457]\n",
      "The theoretical 2nd excited state energy is: 0.8523274042087416\n",
      "The estimated 3rd excited state energy is: [1.47231932]\n",
      "The theoretical 3rd excited state energy is: 1.4726158788876045\n"
     ]
    }
   ],
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:12:45.991342Z",
     "start_time": "2021-04-30T09:12:45.976287Z"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "_______\n",
    "\n",
    "## References\n",
    "\n",
    "[1] Peruzzo, A. et al. A variational eigenvalue solver on a photonic quantum processor. [Nat. Commun. 5, 4213 (2014).](https://www.nature.com/articles/ncomms5213)\n",
    "\n",
    "[2] McArdle, S., Endo, S., Aspuru-Guzik, A., Benjamin, S. C. & Yuan, X. Quantum computational chemistry. [Rev. Mod. Phys. 92, 015003 (2020).](https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.92.015003)\n",
    "\n",
    "[3] Cao, Y. et al. Quantum chemistry in the age of quantum computing. [Chem. Rev. 119, 10856–10915 (2019).](https://pubs.acs.org/doi/abs/10.1021/acs.chemrev.8b00803)\n",
    "\n",
    "[4] Nakanishi, K. M., Mitarai, K. & Fujii, K. Subspace-search variational quantum eigensolver for excited states. [Phys. Rev. Res. 1, 033062 (2019).](https://journals.aps.org/prresearch/pdf/10.1103/PhysRevResearch.1.033062)\n",
    "\n",
    "[5] Vatan, F. & Williams, C. Optimal quantum circuits for general two-qubit gates. [Phys. Rev. A 69, 032315 (2004).](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.69.032315)"
   ],
   "metadata": {}
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3.8.10 64-bit ('paddle_quantum_test': conda)"
  },
  "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.8.10"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "interpreter": {
   "hash": "1ba3360425d54dc61cc146cb8ddc529b6d51be6719655a3ca16cefddffc9595a"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}