{ "cells": [ { "cell_type": "markdown", "id": "7d09a637", "metadata": {}, "source": [ "# Hamiltonian Simulation with Product Formula\n", " Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved. " ] }, { "cell_type": "markdown", "id": "280b8be1", "metadata": {}, "source": [ "## Overview\n", "\n", "In quantum mechanics, a system's energy is described with a Hamiltonian operator $H$. Solving all or partial properties of the Hamiltonian of a given system constitutes a central problem in many disciplines, including condensed matter physics, computational chemistry, high-energy physics, etc. However, the degrees of freedom of a quantum system grow exponentially with its system size, which leads to the inability to effectively simulate quantum systems using classical computers - the quantum state of several hundred qubits cannot be directly stored even with all the storage resources in the world. Unlike classical computers, quantum computers perform all operations directly on the exponentially large Hilbert space, thus having a natural advantage over classical computer when simulating a quantum system. Matter of fact, designing a controlled quantum system to efficiently simulate quantum systems in nature was Feynman's original idea when he first introduced the concept of quantum computing in the 1980s:\n", " \n", "> _\"Nature isn't classical, dammit, and if you want to make a simulation of nature, you'd better make it quantum mechanical, and by golly it's a wonderful problem, because it doesn't look so easy.\"_\n", ">\n", "> --- \"Simulating physics with computers\", 1982, Richard P. Feynman [1]\n", "\n", "\n", "The development of universal quantum computers and a series of quantum simulators has made it possible to realize Feynman's vision. Digital quantum simulation on a universal quantum computer (i.e. quantum simulation by constructing quantum circuits through quantum gates) is considered to be to have the largest potential due to its scalability and generality.\n", "\n", "In this tutorial, we introduce Hamiltonian simulation in Paddle Quantum. It will be divided into three parts:\n", "1. How to construct a system's Hamiltonian using `Hamiltonian` class.\n", "2. How to create the time-evolving circuit with `construct_trotter_circuit()` function.\n", "3. The Suzuki product formula algorithm and how to create its corresponding circuit up to arbitrary order.\n", "\n", "## Define the system's Hamiltonian \n", "Before demoing how to construct a time-evolving circuit, we will first introduce to readers how to construct a `Hamiltonian` object in Paddle Quantum. Users could create a `Hamiltonian` object by specifying a list of Pauli string containing the coefficients and Pauli operators of each term. First let's consider a simple Hamiltonian:\n", "\n", "$$\n", "H = Z \\otimes Z\n", "\\tag{1}\n", "$$\n", "\n", "This Hamiltonian describes a simple interaction between two qubits: when both qubits are in $|0\\rangle$ or $|1\\rangle$ , the energy of the system is $+1$; conversely when the two qubits are in different states, the energy of the system is $-1$.\n", "\n", "The user could construct this Hamiltonian by:" ] }, { "cell_type": "code", "execution_count": 14, "id": "b1546d21", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.0 Z0, Z1\n" ] } ], "source": [ "import warnings\n", "warnings.filterwarnings(\"ignore\")\n", "from paddle_quantum.hamiltonian import Hamiltonian\n", "\n", "h = Hamiltonian([[1, 'Z0, Z1']])\n", "print(h)" ] }, { "cell_type": "markdown", "id": "99ff7830", "metadata": {}, "source": [ "The `Hamiltonian` class in Paddle Quantum supports automatic merging of equal terms, addition, subtraction, indexing, and splitting:" ] }, { "cell_type": "code", "execution_count": 2, "id": "82fae2f8", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.0 Z0, Z1\n" ] } ], "source": [ "h = Hamiltonian([[0.5, 'Z0, Z1'], [0.5, 'Z1, Z0']], compress=True)\n", "print(h)" ] }, { "cell_type": "code", "execution_count": 3, "id": "37ea39a0", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "h + h: 2.0 Z0, Z1\n", "h * 2: 2.0 Z0, Z1\n", "h: 1.0 Z0, Z1\n" ] } ], "source": [ "print('h + h:', h + h)\n", "print('h * 2:', h * 2)\n", "print('h:', h[:])" ] }, { "cell_type": "markdown", "id": "34438c9e", "metadata": {}, "source": [ "The `decompose_pauli_words()` and `decompose_with_sites()` methods can decompose the Hamiltonian into more manageable forms:" ] }, { "cell_type": "code", "execution_count": 4, "id": "8f292d61", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Pauli words decomposition: ([1.0], ['ZZ'])\n", "Pauli with sites decomposition: ([1.0], ['ZZ'], [[0, 1]])\n" ] } ], "source": [ "print('Pauli words decomposition:', h.decompose_pauli_words())\n", "print('Pauli with sites decomposition:', h.decompose_with_sites())" ] }, { "cell_type": "markdown", "id": "10254485", "metadata": {}, "source": [ "In addition, `construct_h_matrix()` will construct its matrix in the $Z$ Pauli basis:" ] }, { "cell_type": "code", "execution_count": 5, "id": "73dac520", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "matrix([[ 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],\n", " [ 0.+0.j, -1.+0.j, 0.+0.j, 0.+0.j],\n", " [ 0.+0.j, 0.+0.j, -1.+0.j, 0.+0.j],\n", " [ 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j]], dtype=complex64)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "h.construct_h_matrix()" ] }, { "cell_type": "markdown", "id": "3a320cb9", "metadata": {}, "source": [ "## Simulate the time evolution\n", "\n", "According to one of the fundamental axioms of quantum mechanics, the evolution of a system over time can be described by\n", "\n", "$$\n", "i \\hbar \\frac{\\partial}{\\partial t} | \\psi \\rangle = H | \\psi \\rangle,\n", "\\tag{2}\n", "$$\n", "\n", "$\\hbar$ is the reduced Planck constant. This equation is the well-known Schrödinger equation. Thus, for a time independent Hamiltonian, the time evolution equation of the system can be written as\n", "\n", "$$\n", "|\\psi(t) \\rangle = U(t) | \\psi (0) \\rangle, ~ U(t) = e^{- i H t}.\n", "\\tag{3}\n", "$$\n", "\n", "Here we take the natural unit $\\hbar=1$ and $U(t)$ is the time evolution operator. The idea of simulating the time evolution process with quantum circuits is to approximate this time evolution operator using the unitary transformation constructed by quantum circuits. In Paddle Quantum, we provide the `construct_trotter_circuit(circuit, Hamiltonian)` function to construct a time-evolving circuit corresponding to a Hamiltonian. Now, let us test it with the Hamiltonian we just constructed:" ] }, { "cell_type": "code", "execution_count": 6, "id": "f983834e", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "--*-----------------*--\n", " | | \n", "--x----Rz(2.000)----x--\n", " \n" ] } ], "source": [ "from paddle_quantum.trotter import construct_trotter_circuit\n", "from paddle_quantum.ansatz import Circuit\n", "\n", "cir = Circuit()\n", "construct_trotter_circuit(cir, h, tau=1, steps=1) \n", "print(cir)" ] }, { "cell_type": "markdown", "id": "988ffa8b", "metadata": {}, "source": [ "We can see that a quantum circuit has been constructed for `h`, which can simulate the time evolution of arbitrary time length based on the input `tau`.\n", "\n", "By calculating the matrix form of this circuit, one can see that it is identical to the time evolution operator $e^{-iHt}$. Here, we use `gate_fidelity` to calculate the fidelity between the unitary matrix of the quantum circuit and the unitary matrix of the time evolution operator. These two processes are identical when the fidelity is equal to 1. We note that a more formal definition of simulation error will be introduced at the end of this section. For now, let's consider fidelity as the criteria of similarity between two evolution processes (unitary operators)." ] }, { "cell_type": "code", "execution_count": 15, "id": "52f3dff9", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The fidelity is: 0.471\n" ] } ], "source": [ "from scipy import linalg\n", "from paddle_quantum.qinfo import gate_fidelity\n", "\n", "# calculate the fidelity between e^{-iHt} and the unitary matrix of circuit\n", "print('The fidelity is: %.3f' \n", " % gate_fidelity(cir.unitary_matrix(), linalg.expm(-1 * 1j * h.construct_h_matrix())))" ] }, { "cell_type": "markdown", "id": "075faf9a", "metadata": {}, "source": [ "Actually, this is because any rotation associated with a tensor product of the pauli operators can be efficiently decomposed into a circuit. In this example, we could change the angle of the Rz gate to simulate any $e^{-i Z\\otimes Z t}$ evolutionary operator. Does this mean that the time-evolving operator of any Pauli Hamiltonian could be perfectly simulated with a quantum circuit? Unfortunately, the answer is negative. Let us consider a slightly more complicated Hamiltonian:\n", "\n", "$$\n", "H = Z \\otimes Z + X \\otimes I + I \\otimes X.\n", "\\tag{4}\n", "$$\n", "\n", "Similarly, let's use `construct_trotter_circuit` to construct its corresponding time-evolving circuit:" ] }, { "cell_type": "code", "execution_count": 8, "id": "8ec1213d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "--*-----------------*----Rx(2.000)--\n", " | | \n", "--x----Rz(2.000)----x----Rx(2.000)--\n", " \n", "The fidelity is: 0.681\n" ] } ], "source": [ "h_2 = Hamiltonian([[1, 'Z0, Z1'], [1, 'X0'], [1, 'X1']]) # no need to write out unit operator\n", "cir = Circuit()\n", "construct_trotter_circuit(cir, h_2, tau=1, steps=1)\n", "print(cir)\n", "print('The fidelity is: %.3f' \n", " % gate_fidelity(cir.unitary_matrix().numpy(), linalg.expm(-1 * 1j * h_2.construct_h_matrix())))" ] }, { "cell_type": "markdown", "id": "f75bd0c7", "metadata": {}, "source": [ "This time the fidelity is less than $1$, which means the circuit cannot correctly simulate the time-evolution process of the system.\n", "\n", "The reason is that, the unitary transformation of the circuit is $e^{-iZ\\otimes Z t} e^{-i (X\\otimes I + I\\otimes X)t}$, while the time evolution operator is $e^{-iZ\\otimes Z t - i(X\\otimes I + I\\otimes X)t}$. And for a quantum system, $e^{A+B} \\neq e^A e^B$ when $[A, B] \\neq 0$. Here, since $[X, Z] \\neq 0$, the corresponding unitary transformation of the circuit is not equal to the correct time-evolution operator.\n", "\n", "In addition to using the fidelity to describe the similarity between the quantum circuit and the time-evolving operator that one wishes to simulate, one can define the error $\\epsilon$ as follows\n", "\n", "$$\n", "\\epsilon(U) = \\Vert e^{-iHt} - U \\Vert,\n", "\\tag{5}\n", "$$\n", "\n", "where $\\Vert \\cdot \\Vert$ denotes the mode of the largest eigen (singular) value. Such a definition better describes the deviation of the quantum state under different evolution operators and it is a more rigorous way to define the simulation time evolution error. We note that the simulation error of this form will be used repeatedly in the next section." ] }, { "cell_type": "markdown", "id": "3109b3f5", "metadata": {}, "source": [ "### Product formula and Suzuki decomposition\n", "\n", "In 1996, Seth Lloyd showed that the error in simulating time evolution can be reduced by splitting a whole evolution time $t$ into $r$ shorter \"time blocks\" [2]. Consider a more general Hamiltonian of the form $H = \\sum_{k=1}^{L} h_k$, where $h_k$ acts on a part of the system. By Taylor expansion, it is not difficult to find that the simulation error is a second-order term, i.e.\n", "\n", "$$\n", "e^{-iHt} = \\prod_{k=1}^{L} e^{-i h_k t} + O(t^2).\n", "\\tag{6}\n", "$$\n", "\n", "Let $\\tau = t/r$ and consider the evolution operator $\\left(e^{-iH \\tau}\\right)^r$, then its error can be derived from \n", "\n", "$$\n", "e^{-iHt} = \\left(e^{-iH \\tau}\\right)^r = \\left(\\prod_{k=1}^{L} e^{-i h_k \\tau} + O(\\tau^2) \\right)^r = \\left(\\prod_{k=1}^{L} e^{-i h_k \\tau} \\right)^r + O\\left(\\frac{t^2}{r}\\right).\n", "\\tag{7}\n", "$$\n", "\n", "The above equation states that an arbitrarily high simulation accuracy can be achieved by splitting the whole evolution time into enough \"pieces\". This is the basic idea of the product formula. However, the error given in (7) is only a rough estimate. In practice, in order to estimate the depth of the quantum circuit required to achieve a certain simulation accuracy, an exact upper bound on the error needs to be further computed. In the following, we show a relatively abbreviated procedure for calculating the error upper bound and readers who are not interested in details can skip to the conclusion on the error bound in (11).\n", "\n", "Let us note the remainder Taylor expansion of the function $f$ up to order $k$ to be $\\mathcal{R}_k(f)$. And the two following statements are needed for the calculating of the error bound:\n", "\n", "$$\n", "\\left\\Vert \\mathcal{R}_k \\left( \\prod_{k=1}^{L} \\exp (-i h_k \\tau) \\right) \\right\\Vert\n", "\\leq\n", "\\mathcal{R}_k \\left( \\exp \\left( \\sum_{k=1}^{L} \\vert \\tau \\vert \\cdot \\Vert h_k \\Vert \\right) \\right),\n", "\\tag{8}\n", "$$\n", "\n", "$$\n", "\\vert \\mathcal{R}_k(\\exp (\\alpha)) \\vert \\leq \\frac{\\vert \\alpha \\vert^{k+1}}{(k+1)!} \\exp ( \\vert \\alpha \\vert ), ~\n", "\\forall \\alpha \\in \\mathbb{C}.\n", "\\tag{9}\n", "$$\n", "\n", "We omit the proofs of these two statements due to length limitations and interested reader could refer to Section F.1 in [3]. As defined in (5), the simulation error can be written as\n", "\n", "$$\n", "\\epsilon\\left(e^{-iH\\tau}, U_{\\rm circuit}\\right) = \\left \\Vert \\exp\\left(-i\\sum_{k=1}^L h_k \\tau\\right) - \\prod_{k=1}^{L} \\exp\\left(-i h_k \\tau \\right) \\right \\Vert.\n", "\\tag{10}\n", "$$\n", "\n", "We already know that simulation error is the reminder of the time-evolving operators' Taylor expansion to the first order. Then using (8), (9) and the triangular inequality, the upper bound on the error in (10) can be calculated as follows\n", "\n", "$$\n", "\\begin{aligned}\n", "\\left \\Vert \\exp\\left(-i\\sum_{k=1}^L h_k \\tau\\right) - \\prod_{k=1}^{L} \\exp\\left(-i h_k \\tau \\right) \\right \\Vert\n", "=~&\n", "\\left \\Vert \\mathcal{R}_1 \\left( \\exp \\left( -i \\sum_{k=1}^{L} h_k \\tau \\right) - \\prod_{k=1}^{L} \\exp (-i h_k \\tau) \\right) \\right \\Vert\n", "\\\\\n", "\\leq~&\n", "\\left \\Vert \\mathcal{R}_1 \\left( \\exp \\left( -i \\sum_{k=1}^{L} h_k \\tau \\right) \\right) \\right \\Vert\n", "+\n", "\\left \\Vert \\mathcal{R}_1 \\left( \\prod_{k=1}^{L} \\exp (-i h_k \\tau) \\right) \\right \\Vert\n", "\\\\\n", "\\leq~ &\n", "\\left \\Vert \\mathcal{R}_1 \\left( \\exp \\left( \\vert \\tau \\vert \\cdot \\left \\Vert \\sum_{k=1}^{L} h_k \\right \\Vert \\right) \\right) \\right \\Vert\n", "+ \n", "\\left \\Vert \\mathcal{R}_1 \\left( \\exp \\sum_{k=1}^{L} \\left( \\vert \\tau \\vert \\cdot \\Vert h_k \\Vert \\right) \\right) \\right \\Vert\n", "\\\\\n", "\\leq~&\n", "2 \\mathcal{R}_1 \\left( \\exp ( \\vert \\tau \\vert L \\Lambda ) \\right )\n", "\\\\\n", "\\leq~&\n", " ( \\vert \\tau \\vert L \\Lambda )^2 \\exp ( \\vert \\tau \\vert L \\Lambda ),\n", "\\end{aligned}\n", "\\tag{11}\n", "$$\n", "\n", "with $\\Lambda = \\max_k \\Vert h_k \\Vert$. Considering the complete evolution time $t = r \\cdot \\tau$, the error when simulating a time evolution operator of length $t$ is\n", "\n", "$$\n", "\\begin{aligned}\n", "\\left \\Vert \\left ( \\exp\\left(-i\\sum_{k=1}^L h_k \\tau \\right)\\right)^r - \\left (\\prod_{k=1}^{L} \\exp\\left(-i h_k \\tau \\right) \\right)^r \\right \\Vert\n", "\\leq ~&\n", "r \\left \\Vert \\exp\\left(-i\\sum_{k=1}^L h_k \\tau\\right) - \\prod_{k=1}^{L} \\exp\\left(-i h_k \\tau \\right) \\right \\Vert\n", "\\\\\n", "=~& r ( \\tau L \\Lambda )^2 \\exp ( \\vert \\tau \\vert L \\Lambda )\n", "\\\\\n", "=~& \\frac{( t L \\Lambda )^2}{r} \\exp \\left( \\frac{\\vert t \\vert L \\Lambda}{r} \\right).\n", "\\end{aligned}\n", "\\tag{12}\n", "$$\n", "\n", "Here we use the conclusion of linear accumulation of errors in quantum circuits, i.e. $\\Vert U^r - V^r \\Vert \\leq r\\Vert U - V \\Vert$, and the reader who is not familiar with this conclusion can refer to Section 4.5.3 in [4]. At this point, we have calculated an upper bound on the simulation error of the product formula for a complete period of evolution time $t$, i.e., the second-order term $O(t^2/r)$ in Eq. (7). \n", "\n", "In fact, we can further optimize the simulation accuracy for the time-evolution operator $e^{-iH\\tau}$ within each \"time block\" by the Suzuki decomposition. For the Hamiltonian $H = \\sum_{k=1}^{L} h_k$, the Suzuki decomposition of the time evolution operator can be written as\n", "\n", "$$\n", "\\begin{aligned}\n", "S_1(\\tau) &= \\prod_{k=0}^L \\exp ( -i h_k \\tau),\n", "\\\\\n", "S_2(\\tau) &= \\prod_{k=0}^L \\exp ( -i h_k \\frac{\\tau}{2})\\prod_{k=L}^0 \\exp ( -i h_k \\frac{\\tau}{2}),\n", "\\\\\n", "S_{2k}(\\tau) &= [S_{2k - 2}(p_k\\tau)]^2S_{2k - 2}\\left( (1-4p_k)\\tau\\right)[S_{2k - 2}(p_k\\tau)]^2,\n", "\\end{aligned}\n", "\\tag{13}\n", "$$\n", "\n", "For $k > 1, k\\in\\mathbb{Z}$, where $p_k = 1/(4 - 4^{1/(2k - 1)})$. The previously derived product formula actually uses only the first-order Suzuki decomposition $S_1(\\tau)$ to simulate each \"time block\". Therefore it's also known as the first-order Suzuki product formula, or simply the first-order product formula. In some scenarios, the Suzuki product formula is also referred to as the Trotter-Suzuki decomposition. For higher-order product formulas, using similar calculations as in (10-12), it can be shown that the error bound on the error for the $2k$-th order product formula is:\n", "\n", "$$\n", "\\epsilon\\left(e^{-iHt}, \\left(S_{2k}(\\tau)\\right)^r\\right)\n", "\\leq\n", "\\frac{(2L5^{k-1}\\Lambda\\vert t \\vert)^{2k+1}}{3r^{2k}} \\exp \\left( \\frac{2L5^{k-1}\\Lambda\\vert t \\vert}{r} \\right),\n", "~ k > 1.\n", "\\tag{14}\n", "$$\n", "\n", "With the upper bound on the simulation error obtained, it is possible to further calculate the lower bound on the circuit depth required to reach a certain minimum accuracy $\\epsilon$. It should be noted that the error bounds given in (12) and (14) are calculated very loosely. In recent years, many works have gone further to give tighter upper bounds [3, 5-6]. In addition, product formulas that are not based on the Suzuki decomposition have also been proposed [7]." ] }, { "cell_type": "markdown", "id": "f3bf6fb2", "metadata": {}, "source": [ "![image.png](./figures/trotter_suzuki_circuit.png)\n", "