tutorials/quantum_simulation/QDRIFT_CN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 利用 qDRIFT 模拟时间演化\n",
    "<em> Copyright (c) 2022 Institute for Quantum Computing, Baidu Inc. All Rights Reserved. </em>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 概述\n",
    "量子力学中系统的能量由哈密顿量算符 $H$ 描述，它决定了系统演化的性质，模拟哈密顿量的时间演化，在建模复杂的化学和物理系统方面具有巨大的实用价值。然而，由于系统的自由度随系统（如量子比特数）增大呈指数级增加，导致一般情况下无法利用经典计算机有效模拟量子系统。当前使用量子计算机模拟哈密顿量的主要技术是利用乘积式方法（product formula）模拟时间演化，本教程将介绍关于 product formula 的一些基础理论和方法，和基于 product formula 改进的 quantum stochastic drift protocol (qDRIFT) —— 一种随机的 product formula 方法，并在文末进行了代码演示。\n",
    "\n",
    "## 利用 Product Formula 模拟时间演化\n",
    "\n",
    "根据量子力学的基本公理，在确定了一个系统的哈密顿量 $H$ 之后，该系统随时间演化的过程可以由如下方程描述\n",
    "\n",
    "$$\n",
    "i \\hbar \\frac{d}{d t} | \\psi \\rangle = H | \\psi \\rangle,\n",
    "\\tag{1}\n",
    "$$\n",
    "\n",
    "其中 $\\hbar$ 为约化普朗克常数。因此，对于一个不含时的哈密顿量，系统的时间演化方程可以写为\n",
    "\n",
    "$$\n",
    "|\\psi(t) \\rangle = U(t) | \\psi (0) \\rangle, ~ U(t) = e^{- i H t}.\n",
    "\\tag{2}\n",
    "$$\n",
    "\n",
    "这里我们取自然单位 $\\hbar=1$，$U(t)$ 为时间演化算符。利用量子电路来模拟时间演化过程的核心思想是利用量子电路构建出的酉变换模拟和近似该时间演化算符。Seth Lloyd 在其 1996 年的文章中指出，可以将一整段的演化时间 $t$ 拆分为 $r$ 份较短的“时间片段”来减小模拟时间演化的误差 [1]。考虑一个一般的哈密顿量形式 $H = \\sum_{k=1}^{L} H_k$，其中 $H_k$ 是作用在部分系统上的子哈密顿量。我们考虑每个子哈密顿量 $H_k$ 的演化算符为$e^{-i H_k t}$，我们依次模拟每个子哈密顿量可以得到 $\\prod_{k=1}^{L} e^{-i H_k t}$。通过泰勒展开，可以发现\n",
    "\n",
    "$$\n",
    "e^{-iHt} = \\prod_{k=1}^{L} e^{-i H_k t} + O(t^2).\n",
    "\\tag{3}\n",
    "$$\n",
    "\n",
    "那么，我们令 $\\tau = t/r$，并考虑演化算符 $\\left(e^{-iH \\tau}\\right)^r$，就可以推导出\n",
    "\n",
    "$$\n",
    "e^{-iHt} = \\left(e^{-iH \\tau}\\right)^r = \\left(\\prod_{k=1}^{L} e^{-i H_k \\tau} + O(\\tau^2) \\right)^r = \\left(\\prod_{k=1}^{L} e^{-i H_k \\tau} \\right)^r + O\\left(\\frac{t^2}{r}\\right).\n",
    "\\tag{4}\n",
    "$$\n",
    "\n",
    "上式告诉我们，只要将一整段演化时间拆为足够多的“片段”，就可以达到任意高的模拟精度，这就是 product formula 的基本思想。不过，(4) 中给出的只是一个粗略的估计。如果我们想要估计达到某一模拟精度所需要的量子电路深度，就需要推导其更严格的误差上界。具体地，我们令 $U_{circuit}$ 代表我们构造的电路，$\\Vert \\cdot \\Vert$ 为 Schatten-$\\infty$ 范数，即[谱范数](https://en.wikipedia.org/wiki/Schatten_norm)。那么其模拟误差 $\\epsilon$ 可以写为\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\epsilon\\left(e^{-iHt}, U_{circuit}\\right)  & = \\Vert e^{-iHt} - U_{circuit}\\Vert .\n",
    "\\end{aligned}\n",
    "\\tag{5}\n",
    "$$\n",
    "下面，我们展示一个比较简略的误差上界计算过程，我们不加证明地列出 (6)、(7) 两个结论，会在证明 (8) 时用到，感兴趣的读者可以参考 [2] 中的 F.1 节获取证明细节。\n",
    "$$\n",
    "\\left\\Vert \\mathcal{R}_k \\left( \\prod_{k=1}^{L} e^{-i H_k \\tau} \\right) \\right\\Vert\n",
    "\\leq\n",
    "\\mathcal{R}_k \\left( e^{\\vert \\tau \\vert   \\sum_{k=1}^{L} \\Vert H_k \\Vert } \\right),\n",
    "\\tag{6}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\vert \\mathcal{R}_k(e^\\alpha) \\vert \\leq \\frac{\\vert \\alpha \\vert^{k+1}}{(k+1)!}  e^{ \\vert \\alpha \\vert }, ~\n",
    "\\forall \\alpha \\in \\mathbb{C},\n",
    "\\tag{7}\n",
    "$$\n",
    "\n",
    "其中 $\\mathcal{R}_k(f)$为函数 $f$ 泰勒展开至 $k$ 阶之后的余项，例如 $\\mathcal{R}_1 (e^x)=\\mathcal{R}_1 (\\sum_{j=0}^\\infty \\frac{x^n}{n!})=\\sum_{j=2}^\\infty \\frac{x^n}{n!}$。\n",
    "令 $\\Lambda = \\max_k \\Vert H_k \\Vert$，考虑完整的演化时间 $t = r \\cdot \\tau$，那么模拟长度为 $t$ 的时间演化算符时的误差为：\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\left \\Vert \\left ( e^{-i\\tau \\sum_{k=1}^L H_k  }\\right)^r - \\left (\\prod_{k=1}^{L} e^{-i H_k \\tau} \\right)^r \\right \\Vert \\leq &\n",
    "r \\left \\Vert e^{-i\\tau \\sum_{k=1}^L H_k } - \\prod_{k=1}^{L} e^{-i H_k \\tau } \\right \\Vert  \\\\\n",
    "=& r \\left \\Vert \\mathcal{R}_1 \\left(  e^{-i\\tau \\sum_{k=1}^L H_k} \\right)- \\mathcal{R}_1 \\left( \\prod_{k=1}^{L} e^{-i H_k \\tau } \\right) \\right \\Vert \\\\\n",
    "\\leq& r \\left \\Vert \\mathcal{R}_1 \\left(  e^{-i\\tau \\sum_{k=1}^L H_k} \\right) \\right \\Vert+ r\\left \\Vert \\mathcal{R}_1 \\left( \\prod_{k=1}^{L} e^{-i H_k \\tau } \\right) \\right \\Vert \\\\\n",
    "\\leq& 2r \\left \\Vert \\mathcal{R}_1 \\left(  e^{-i |\\tau | \\sum_{k=1}^L \\Vert H_k \\Vert} \\right) \\right \\Vert \\\\\n",
    "\\leq& 2r \\left \\Vert \\mathcal{R}_1 \\left(  e^{-i |\\tau | L \\Lambda} \\right) \\right \\Vert \\\\\n",
    "\\leq& r (  \\tau L \\Lambda )^2 e^{\\vert \\tau \\vert L \\Lambda } \\\\\n",
    "=&\\frac{(  t L \\Lambda )^2}{r} e^{\\frac{\\vert t \\vert L \\Lambda}{r} }.\n",
    "\\end{aligned}\n",
    "\\tag{8}\n",
    "$$\n",
    "\n",
    "其中这里用到了量子电路中误差线性累积的结论，即 $\\Vert U^r - V^r \\Vert \\leq r\\Vert U - V \\Vert$，不熟悉这一结论的读者可以参考 [3] 中的 4.5.3 节；也用到了 (7) 式中 $k=1$ 时的结论。至此，我们就计算出了 product formula 对于一段完整的演化时间 $t$ 的模拟误差上界，即 (4) 式中的二阶项 $O(t^2/r)$。 \n",
    "\n",
    "\n",
    "在得到了模拟误差上界的基础上，便可以进一步计算达到一定精度 $\\epsilon$ 时所需要的电路深度的下界。从 (8) 中我们不难发现，式子里含有 $L$ 项，这就意味着，随着哈密顿量项数的增加，若需控制误差上界，则时间片段的划分必须越来越细，这就使得电路深度增加。本文所要介绍的 qDRIFT 在一定程度上解决了该问题。qDRIFT 着眼于哈密顿量本身的系数，将其建模为一个概率分布，每次从该概率分布中采样酉门并重复一定的次数，从而构成量子电路，最终在给定的精度下，其量子电路的深度将不显含哈密顿量项数 $L$。下面我们将详细介绍它。\n",
    "\n",
    "\n",
    "## 利用 qDRIFT 模拟时间演化\n",
    "首先，我们给出目标哈密顿量的形式\n",
    "$$\n",
    "H=\\sum_{j=1}^L h_j H_j,\n",
    "\\tag{9}\n",
    "$$\n",
    "它含有 $L$ 项子哈密顿量 $H_j$，值得注意的是，这里的 $H_j$ 是已经被归一化了的，也就是 $\\Vert H_j \\Vert = 1$，其中 $\\Vert\\cdot\\Vert$ 为 Schatten-$\\infty$ 范数 ，$h_j$ 是每个子哈密顿量的系数，它是一个正实数。通过该系数我们便可以构造一个离散的概率分布，以单个系数在整个哈密顿量系数总和的占比作为每个酉门被采样的概率，也就是 $p_j =h_j / \\lambda$，其中 $\\lambda =\\sum _j h_j$ 是系数和，如此采样重复 $N$ 次（为了与 product formula 对照，我们取 $ N=Lr$），我们就得到一个由 $j$ 排列的有序列表，并可以根据该排列构造酉门 $U_j = e^{i\\tau H_j}$ 。假设 $L=3$ ，$r=2$，我们可以根据上述概率分布采样一个有序列表形如\n",
    "$$\n",
    "[ 3, 1, 2 ,3 ,3 ,1 ],\n",
    "$$\n",
    "那么就可以据此构造量子电路\n",
    "$$\n",
    "U_{circuit} = e^{i\\tau H_1}e^{i\\tau H_3}e^{i\\tau H_3}e^{i\\tau H_2}e^{i\\tau H_1}e^{i\\tau H_3},\n",
    "$$\n",
    "$\\tau=t\\lambda /N$，这就是 qDRIFT 模拟哈密顿量的一个实现。\n",
    "\n",
    "qDRIFT 的实现流程非常简单，而它的优势在于，在给定目标精度 $\\epsilon$ 时其酉门数的复杂度为 $O((\\lambda t)^2 /\\epsilon) $，可以看到这是一个不含 $L$ 的结果，也就是说，其酉门数量与哈密顿量的项数不显式相关，这在哈密顿量项数很大时可以有效地缩减模拟电路的长度。接下来我们将给出证明。\n",
    "\n",
    "我们将根据概率分布进行采样的过程建模为一个量子信道，我们用花体字母 $\\mathcal{E}$ 和 $\\mathcal{U}$ 分别来表示通过 qDRIFT 建立的信道和所要模拟的信道，并且用 $\\mathcal{E}_N$ 和 $\\mathcal{U}_N$ 代表其各自信道对量子态 $\\rho$ 的 $N$ 次作用中的一次作用，即\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "&\\mathcal{U}_N (\\rho) = e^{\\frac{it}{N}H} \\rho e^{\\frac{-it}{N}H}= e^{\\frac{t}{N}\\mathcal{L}} (\\rho),\n",
    "\\\\\n",
    "&\\mathcal{E}_N (\\rho)=\\sum_{j}p_j e^{i\\tau H_j} \\rho e^{-i\\tau H_j}=\\sum_{j} p_j e^{\\tau \\mathcal{L}_j}(\\rho).\n",
    "\\end{aligned}\n",
    "\\tag{10}\n",
    "$$\n",
    "\n",
    "这里我们引入 Liouvillian 表示 ，即对量子信道 $\\mathcal{P}(\\rho)=e^{iHt}\\rho e^{-iHt}$ 有\n",
    "$$\n",
    "\\mathcal{P}(\\rho)=e^{iHt}\\rho e^{-iHt}=e^{t\\mathcal{L}}(\\rho)=\\sum_{k=0}^\\infty \\frac{t^k \\mathcal{L}^k (\\rho)}{k!},\n",
    "\\tag{11}\n",
    "$$\n",
    "\n",
    "\n",
    "其中 $\\mathcal{L}(\\rho)=i(H\\rho - \\rho H)$ ，同理有 $\\mathcal{L}_j(\\rho)=i(H_j\\rho - \\rho H_j)$ 。需要注意的是，$\\mathcal{L}$ 的迭代规则遵循 $\\mathcal{L}^{n+1}(\\rho)=i(H\\mathcal{L}^n(\\rho)-\\mathcal{L}^n(\\rho)H)$。具体来说，$\\mathcal{U}_N = \\sum_{n=0}^\\infty \\frac{t^n\\mathcal{L}^n}{n!N^n}$，$\\mathcal{E}_N =\\sum_{j}p_j \\sum_{n=0}^\\infty \\frac{\\lambda^n t^n \\mathcal{L}_j^n}{n!N^n}$。接下来我们该如何度量两个信道的距离呢？这里引入[菱形范数](https://en.wikipedia.org/wiki/Diamond_norm) (diamond norm) 的定义式\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\Vert \\mathcal{P} \\Vert_\\Diamond :=\\sup_{\\rho ; \\Vert \\rho \\Vert _1 =1}\\Vert (\\mathcal{P} \\otimes \\mathbb{I})(\\rho )\\Vert _1 .\n",
    "\\end{aligned}\n",
    "\\tag{12}\n",
    "$$\n",
    "其中 $\\mathbb{I}$ 为与 $\\mathcal{P}$ 空间相同大小的单位信道，$\\Vert \\cdot \\Vert_1$ 为 Schatten-$1$ 范数，即[迹范数](https://en.wikipedia.org/wiki/Schatten_norm)。我们使用菱形范数定义两个量子信道的距离\n",
    "$$\n",
    "\\begin{aligned}\n",
    "d_\\Diamond (\\mathcal{E},\\mathcal{U}) &=\\frac{1}{2} \\Vert \\mathcal{E} -\\mathcal{U} \\Vert_\\Diamond\n",
    "\\\\\n",
    "&=\\sup_{\\rho ; \\Vert \\rho \\Vert _1 =1} \\frac{1}{2} \\Vert ((\\mathcal{E}-\\mathcal{U}) \\otimes \\mathbb{I})(\\rho )\\Vert _1 .\n",
    "\\end{aligned}\n",
    "\\tag{13}\n",
    "$$\n",
    "菱形范数代表了在所有量子态中能够分辨两个信道的最大可能性，它的值越大，两个信道被区分的可能性就越大，也就代表了两个信道距离远，模拟效果差；反之它的值小，就代表模拟效果好。接着，我们可以去计算单次作用的信道的距离上界\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    \\Vert \\mathcal{U}_N-\\mathcal{E}_N \\Vert_\\Diamond &= \\left\\Vert \\sum_{n=2}^\\infty \\frac{t^n\\mathcal{L}^n}{n!N^n}-\\sum_{j}\\frac{h_j}{\\lambda} \\sum_{n=2}^\\infty \\frac{\\lambda^n t^n \\mathcal{L}_j^n}{n!N^n} \\right\\Vert_\\Diamond\\\\\n",
    "    &\\leq \\sum_{n=2}^\\infty \\frac{t^n\\Vert \\mathcal{L}^n \\Vert_\\Diamond }{n!N^n} + \\sum_{j}\\frac{h_j}{\\lambda} \\sum_{n=2}^\\infty \\frac{\\lambda^n t^n \\Vert\\mathcal{L}_j^n \\Vert_\\Diamond }{n!N^n}\\\\\n",
    "    &\\leq \\sum_{n=2}^\\infty \\frac{1}{n!}\\left( \\frac{2\\lambda t}{N}\\right)^n+\\sum_{j}\\frac{h_j}{\\lambda} \\sum_{n=2}^\\infty \\frac{1}{n!}\\left( \\frac{2\\lambda t}{N}\\right)^n\\\\\n",
    "    &=2\\sum_{n=2}^\\infty \\frac{1}{n!}\\left( \\frac{2\\lambda t}{N}\\right)^n .\n",
    "\\end{aligned}\n",
    "\\tag{14}\n",
    "$$\n",
    "其中这里用到了结论 $\\Vert \\mathcal{L} \\Vert_\\Diamond \\leq 2\\Vert H\\Vert \\leq 2\\lambda$ ，同理有 $\\Vert \\mathcal{L}_j \\Vert_\\Diamond \\leq 2\\Vert H_j\\Vert \\leq 2$ [4]。接着，我们可以利用上文中提到的 (7) 的结论，令 $k=1$，$\\alpha=2\\lambda t /N$ 便可得到\n",
    "\n",
    "$$\n",
    "d_\\Diamond (\\mathcal{U}_N,\\mathcal{E}_N) \\leq \\frac{2\\lambda^2 t^2}{N^2} e^{2\\lambda t/N} ,\n",
    "\\tag{15}\n",
    "$$\n",
    "然后再次利用 $\\Vert U^r - V^r \\Vert \\leq r\\Vert U - V \\Vert$ 这一结论（需要注意，式子中的 $U$ 与 $V$ 本是线性算子，但对于量子信道 $\\mathcal{U}$ 和 $\\mathcal{E}$ 依然适用，感兴趣的读者可以参考 [6] 中的第 3.3.2 节获取证明细节），且通常情况下 $2\\lambda t \\ll N$，便可推出\n",
    "$$\n",
    "\\begin{aligned}\n",
    "d_\\Diamond (\\mathcal{U},\\mathcal{E}) &\\leq N d_\\Diamond (\\mathcal{U}_N, \\mathcal{E}_N)\\\\\n",
    "    &=\\frac{2\\lambda^2 t^2}{N} e^{2\\lambda t/N} \\approx \\frac{2\\lambda^2 t^2}{N}.\n",
    "\\end{aligned}\n",
    "\\tag{16}\n",
    "$$\n",
    "因此 $ N \\sim O((\\lambda t)^2 /\\epsilon)$。 由上式可以看出，在满足 $\\lambda \\ll \\Lambda L$ 的条件下（回忆一下，$\\Lambda = \\max_k \\Vert H_k \\Vert$，qDRIFT 将哈密顿量写为 $H=\\sum_{j=1}^L h_j H_j$，那么对应的 $\\Lambda = \\max_k h_k $），其距离将不与 $L$ 显式相关，这也就可以在 $L$ 较大即情况较为复杂时，不会带来量子电路深度的显著增加，可以有效控制酉门的数量。很多物理系统的哈密顿量都满足 $\\lambda \\ll \\Lambda L$，如乙烷、二氧化碳的电子结构，但并非所有情况都满足，若 $\\lambda = \\Lambda L$ 或 $\\lambda = \\Lambda \\sqrt{L}$ 时，它们的上界分别为 $O(L^2(\\Lambda t)^2 /\\epsilon)$ 和 $O(L(\\Lambda t)^2 /\\epsilon)$ ，可以看到，它们仍然随着哈密顿量项数增大而增大。感兴趣的读者可以参考 [4] 获取更多细节。\n",
    "\n",
    "\n",
    "## 代码实现\n",
    "我们将结合实际代码实现 qDRIFT。我们将首先演示其采样结果的性能，再计算其信道的模拟误差。首先我们需要导入需要的包。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "\n",
    "import math\n",
    "import numpy as np                     \n",
    "import scipy                                           \n",
    "import paddle_quantum as pq         \n",
    "import paddle\n",
    "\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")   # 隐藏 warnings\n",
    "np.set_printoptions(suppress=True, linewidth=np.nan)        # 启用完整显示，便于在终端 print 观察矩阵时不引入换行符\n",
    "pq.set_backend('density_matrix')    # 使用密度矩阵表示"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们假设系统由 2 个 qubits 组成，我们可以利用量桨的 `hamiltonian` 模块构造一个哈密顿量项数为 $L=4$ 的哈密顿量，为了演示 qDRIFT 的效果，我们选择一组满足 $\\lambda \\ll \\Lambda L$ 的参数，这便是我们的目标哈密顿量，具体如下\n",
    "$$\n",
    "\\begin{aligned}\n",
    "H&=I \\otimes X + 0.05 * X \\otimes Z + 0.05 * I \\otimes Y + 0.05 * X \\otimes X   .\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "目标哈密顿量为: \n",
      " [[ 0.  +0.j    1.  -0.05j  0.05+0.j    0.05+0.j  ]\n",
      " [ 1.  +0.05j  0.  +0.j    0.05+0.j   -0.05+0.j  ]\n",
      " [ 0.05+0.j    0.05+0.j    0.  +0.j    1.  -0.05j]\n",
      " [ 0.05+0.j   -0.05+0.j    1.  +0.05j  0.  +0.j  ]]\n"
     ]
    }
   ],
   "source": [
    "qubits = 2  # 设置量子比特数\n",
    "H_j = [(1.0, 'I0,X1'),  # 构造哈密顿量的泡利串\n",
    "       (0.05, 'X0,Z1'),\n",
    "       (0.05, 'I0,Y1'),\n",
    "       (0.05, 'X0,X1'), ]\n",
    "\n",
    "H = pq.hamiltonian.Hamiltonian(H_j)    \n",
    "print(f'目标哈密顿量为: \\n {H.construct_h_matrix(qubit_num=qubits)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来，我们根据 $\\lambda = \\sum_j h_j$，$ p_j=h_j/\\lambda $ 计算概率。在本次实验中，假设我们的目标精度 $\\epsilon=0.1$，模拟时间 $t=1$，也就是说，我们需要采样 $N=\\lceil \\frac{2\\lambda^2 t^2}{\\epsilon}\\rceil = 27$ 次。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "达到 0.1 的精度需要 27 个酉门\n"
     ]
    }
   ],
   "source": [
    "h_j = np.array(H.coefficients)  # 获取系数\n",
    "lamda = h_j.sum()\n",
    "p_j = h_j/lamda  # 计算离散概率分布\n",
    "accuracy = 0.1\n",
    "t = 1\n",
    "gate_counts = math.ceil(2 * lamda**2 * t**2 / accuracy)\n",
    "\n",
    "print(f'达到 {accuracy} 的精度需要 {gate_counts} 个酉门')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接着，我们将根据概率分布 $p_j$ 独立采样 27 次，并根据该采样结果构造酉电路。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "采样结果为:\n",
      " [1 1 1 1 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 4 1]\n",
      "qDRIFT 的模拟电路矩阵为: \n",
      " [[ 0.51998969-0.02531459j  0.03408183+0.85099285j -0.03524884+0.02376227j -0.0353899 +0.02366261j]\n",
      " [-0.03408183+0.85099285j  0.51998969+0.02531459j  0.0353899 +0.02366261j -0.03524884-0.02376227j]\n",
      " [-0.03524884+0.02376227j -0.0353899 +0.02366261j  0.51998969-0.02531459j  0.03408183+0.85099285j]\n",
      " [ 0.0353899 +0.02366261j -0.03524884-0.02376227j -0.03408183+0.85099285j  0.51998969+0.02531459j]] \n",
      "原始电路矩阵为: \n",
      " [[ 0.53752508-0.00075235j  0.04202098+0.83966719j -0.04201839+0.04202098j -0.00075235+0.02697398j]\n",
      " [-0.04202098+0.83966719j  0.53752508+0.00075235j  0.00075235+0.02697398j -0.04201839-0.04202098j]\n",
      " [-0.04201839+0.04202098j -0.00075235+0.02697398j  0.53752508-0.00075235j  0.04202098+0.83966719j]\n",
      " [ 0.00075235+0.02697398j -0.04201839-0.04202098j -0.04202098+0.83966719j  0.53752508+0.00075235j]]\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(666)  # 固定随机数初始位置，便于演示说明\n",
    "sample_list = np.random.choice(a=range(1, 5), size=gate_counts, replace=True, p=p_j)\n",
    "print(f'采样结果为:\\n {sample_list}')\n",
    "\n",
    "# 根据采样结果计算采样出来的酉电路\n",
    "simulation = np.identity(2 ** qubits)  # 生成单位矩阵\n",
    "tau = 1j*lamda*t/gate_counts\n",
    "for i in sample_list:\n",
    "    pauli_str_j = (1.0, H_j[i-1][1])   # 获取H_j，注意，应抛弃其原有系数\n",
    "    H_i = pq.hamiltonian.Hamiltonian([pauli_str_j]).construct_h_matrix(qubit_num=qubits)\n",
    "    simulation = np.matmul(scipy.linalg.expm(tau*H_i), simulation)  \n",
    "origin = scipy.linalg.expm(1j*t*H.construct_h_matrix(qubit_num=qubits))  # 计算目标哈密顿量的原始电路\n",
    "print(f'qDRIFT 的模拟电路矩阵为: \\n {simulation} \\n原始电路矩阵为: \\n {origin}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "然后我们便可计算出从 qDRIFT 采样出来的酉电路和原始电路之间的模拟误差 $\\Vert e^{iHt}-U_{circuit}\\Vert$，注意区分，这里的范数为谱范数。 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "模拟误差为: 0.0309\n"
     ]
    }
   ],
   "source": [
    "distance = 0.5 * np.linalg.norm(origin-simulation, ord=2)\n",
    "print(f'模拟误差为: {distance:.4f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "当然我们可以带入一个具体的量子态试验一下，不失一般性，我们假设初始量子态为零态，即 $\\rho(0)  = | 0 \\rangle \\langle 0 | $，本教程的实验我们均使用密度矩阵描述量子态。我们可以让量子态分别通过原始方法和 qDRIFT 模拟方法演化，到 $t$ 时刻量子态分别为 $\\rho(t)_{origin}$ 和 $\\rho(t)_{qDRIFT}$，最后可以比较这两个量子态的保真度来衡量模拟电路的效果。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "初始量子态为 \n",
      " [[1.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j]]\n",
      "两个量子态之间的保真度为0.9989\n"
     ]
    }
   ],
   "source": [
    "rho_0 = pq.state.zero_state(qubits).numpy() # 构造零态密度矩阵\n",
    "print(f'初始量子态为 \\n {rho_0}')\n",
    "\n",
    "rho_t_origin = pq.state.to_state(origin @ rho_0 @ origin.T.conjugate())  # 经过原始电路演化\n",
    "rho_t_qdrift = pq.state.to_state(simulation @ rho_0 @ simulation.T.conjugate())  # 经过模拟电路演化\n",
    "fidelity = pq.qinfo.state_fidelity(rho_t_origin, rho_t_qdrift)\n",
    "print(f'两个量子态之间的保真度为{float(fidelity):.4f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以发现，上面的测试均符合我们的精度要求。但区别于根据 qDRIFT 方法采样得到的某个具体的酉电路，我们将 qDRIFT 的采样方法看作是一个量子信道，也即对量子态 $\\rho$ 的一个映射。上面的实验只是这个信道的一次具体表达，我们接下来将分析这个信道的性能。我们可以定义一个函数，用于描绘 qDRIFT 信道。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 定义 qDRIFT 信道\n",
    "def qdrift_channel(iter_num, sample_num, hamiltonian_list, coefficient_list, simulation_time, qubits, input_state):\n",
    "    '''\n",
    "    输入 :\n",
    "        iter_num : 当前迭代次数，作为递归的标记\n",
    "        sample_num : 采样次数，即 N\n",
    "        hamiltonian_list : 目标哈密顿量的泡利串形式的列表,即 H_j\n",
    "        coefficient_list : 子哈密顿量的系数列表，即 h_j\n",
    "        simulation_time : 模拟时间，即 t\n",
    "        qubits : 系统的量子比特数\n",
    "        input_state : 输入的量子态，应为密度算子\n",
    "    \n",
    "    输出 :\n",
    "        经过该 qDRIFT 信道的量子态（密度算子表示）\n",
    "    '''\n",
    "    lamda = coefficient_list.sum() \n",
    "    tau = lamda*simulation_time/sample_num\n",
    "    output = 0\n",
    "\n",
    "    if iter_num != 1:   # 在迭代标志不为 1 的时候启用递归\n",
    "        input_state = qdrift_channel(iter_num-1, sample_num, hamiltonian_list,\n",
    "                                     coefficient_list, simulation_time, qubits, input_state)\n",
    "\n",
    "    # 计算 e^{iH\\tau} \\rho e^{-iH\\tau}                                 \n",
    "    for sub_H, sub_h in zip(hamiltonian_list, coefficient_list):\n",
    "        sub_H = pq.hamiltonian.Hamiltonian([sub_H]).construct_h_matrix(qubit_num=qubits)\n",
    "        unitary = scipy.linalg.expm(1j*tau*sub_H)  # 计算 e^{iH\\tau}\n",
    "        output += sub_h/lamda*unitary @ input_state @ unitary.conjugate().T\n",
    "    return output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接着我们便可以通过菱形范数计算两个信道的距离，不过菱形范数的求解可以转换为半正定规划问题，即\n",
    "\n",
    "$$\n",
    "d_\\Diamond(\\mathcal{U}- \\mathcal{E})=\\sup_{\\Omega \\geq 0 \\atop \\rho \\geq 0}\\{\\text{Tr}[\\Omega (\\Gamma_\\mathcal{U}-\\Gamma_\\mathcal{E})]: \\Omega \\leq \\rho \\otimes \\mathbb{I},\\text{Tr} (\\rho)=1\\},\n",
    "\\tag{17}\n",
    "$$\n",
    "其中 $\\Gamma_\\mathcal{U}$ 与 $\\Gamma_\\mathcal{E}$ 为原始信道和模拟信道的 Choi 表示。菱形范数的半正定规划和 Choi 表示有多种形式，感兴趣的读者可以阅读 [6-8] 获取更多细节。我们这里使用的 Choi 表示具体为\n",
    "$$\n",
    "\\Gamma_\\mathcal{P}=\\sum_{i,j=0}^{d-1} |i\\rangle \\langle j| \\otimes \\mathcal{P}(|i\\rangle \\langle j|),\n",
    "\\tag{18}\n",
    "$$\n",
    "其中 $\\mathcal{P}$ 为量子信道，$d$ 为该量子信道输入量子态的维度。这里我们首先计算两个信道的 Choi 表示。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 计算原始信道和 qDRIFT 信道的 Choi 表示，在该表示下可以进而计算菱形范数\n",
    "choi_qdrift = 0\n",
    "choi_origin = 0\n",
    "channel = scipy.linalg.expm(1j*t*H.construct_h_matrix(qubit_num=qubits))\n",
    "for i in range(2 ** qubits):\n",
    "    for k in range(2 ** qubits):\n",
    "        choi_temp = np.zeros((2 ** qubits, 2 ** qubits))\n",
    "        choi_temp[i][k] = 1  # 生成 |i\\rangle \\langle k|\n",
    "\n",
    "        # 分两步计算信道 E 的 Choi 表示\n",
    "        # 先计算 \\mathcal{E}(|i\\rangle \\langle k|）\n",
    "        choi_temp_qdrift = qdrift_channel(gate_counts, gate_counts, H_j, h_j, t, qubits, choi_temp)  \n",
    "        # 再计算 |i\\rangle \\langle k| \\otimes \\mathcal{E}(|i\\rangle \\langle k|）\n",
    "        choi_qdrift += np.kron(choi_temp, choi_temp_qdrift)\n",
    "\n",
    "        # 分两步计算信道 U 的 Choi 表示\n",
    "        # 先计算 \\mathcal{U}(|i\\rangle \\langle k|）\n",
    "        choi_temp_origin = channel @ choi_temp @ channel.T.conjugate()\n",
    "        # 再计算 |i\\rangle \\langle k| \\otimes \\mathcal{U}(|i\\rangle \\langle k|）\n",
    "        choi_origin += np.kron(choi_temp, choi_temp_origin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接着我们可以按照 (17) 式计算菱形范数，并求取两个信道的菱形距离。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "两个信道之间的距离为: 0.0764\n"
     ]
    }
   ],
   "source": [
    "print(f'两个信道之间的距离为: {0.5*pq.qinfo.diamond_norm(paddle.to_tensor(choi_origin-choi_qdrift)):.4f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以看到，计算结果是符合预期的。值得注意的是，该值代表了该信道采样为具体模拟电路的最差表现的期望值，它并不能保证每个采样出来的电路都能够达到该精度。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 小结\n",
    "\n",
    "量子模拟本身是一个比较宽泛的话题，其应用也十分广泛。本教程介绍了 product formula 的理论基础和 qDRIFT 方法，并给出了 qDRIFT 的实现例子。但 qDRIFT 并非随机的 product formula 的唯一方法。作为使用 product formula 进行量子模拟的方法的一个分支，随机的 product formula 还有诸多方法值得我们去探究。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## 参考资料\n",
    " \n",
    "[1] Lloyd, Seth. \"Universal quantum simulators.\" [Science (1996): 1073-1078](https://www.jstor.org/stable/2899535).\n",
    "\n",
    "[2] Childs, Andrew M., et al. \"Toward the first quantum simulation with quantum speedup.\" [Proceedings of the National Academy of Sciences 115.38 (2018): 9456-9461](https://www.pnas.org/content/115/38/9456.short).\n",
    "\n",
    "[3] Nielsen, Michael A., and Isaac Chuang. \"Quantum computation and quantum information.\" (2002): 558-559.\n",
    "\n",
    "[4] Campbell, E. . \"Random Compiler for Fast Hamiltonian Simulation.\" [Physical Review Letters 123.7(2019):070503.1-070503.5](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.123.070503).\n",
    "\n",
    "[5] Khatri, Sumeet, and Mark M. Wilde. \"Principles of quantum communication theory: A modern approach.\" [arXiv preprint arXiv:2011.04672 (2020).](https://arxiv.org/abs/2011.04672)\n",
    "\n",
    "[6] Watrous, J. . [The Theory of Quantum Information](https://cs.uwaterloo.ca/~watrous/TQI/).  2018.\n",
    "\n",
    "[7] Watrous, J. . \"Simpler semidefinite programs for completely bounded norms.\" [Chicago Journal of Theoretical Computer Science (2012).](https://arxiv.org/abs/1207.5726)\n",
    "\n",
    "[8] Watrous, J. . \"Semidefinite Programs for Completely Bounded Norms.\" [Theory of Computing 5.1(2009):217-238.](https://arxiv.org/abs/0901.4709)\n",
    "\n"
   ]
  }
 ],
 "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
}