GIBBS_Tutorial_CN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 吉布斯态的制备 （Gibbs State Preparation）\n",
    "\n",
    "<em> Copyright (c) 2020 Institute for Quantum Computing, Baidu Inc. All Rights Reserved. </em>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 概览\n",
    "- 在这个案例中，我们将展示如何通过Paddle Quantum训练量子神经网络来制备量子吉布斯态。\n",
    "\n",
    "- 让我们通过下面几行代码引入必要的library和package。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import scipy\n",
    "\n",
    "from numpy import array, concatenate, zeros\n",
    "from numpy import pi as PI\n",
    "from numpy import trace as np_trace\n",
    "\n",
    "from paddle import fluid\n",
    "from paddle.complex import matmul, trace\n",
    "from paddle_quantum.circuit import UAnsatz\n",
    "from paddle_quantum.state import density_op\n",
    "from paddle_quantum.utils import state_fidelity, partial_trace, pauli_str_to_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 背景\n",
    "\n",
    "量子计算中的前沿方向包括量子机器学习和量子优化，在这两个方向中，特定量子态的制备是非常重要的问题。特别的，吉布斯态（Gibbs state）的制备是实现诸多量子算法所必须的一个步骤并且广泛应用于：\n",
    "- 量子机器学习中受限波尔兹曼机的学习 [1]\n",
    "- 解决凸优化和半正定规划等优化问题 [2]\n",
    "- 组合优化问题 [3]\n",
    "\n",
    "具体的吉布斯态定义如下：给定一个 $n$ 量子位的哈密顿量 $H$（一般来说这是一个$2^n\\times2^n$的厄米矩阵），其在温度 $T$ 下的吉布斯态为 \n",
    "$$\n",
    "\\rho_G = \\frac{{{e^{ - \\beta H}}}}{{tr({e^{ - \\beta H}})}}\n",
    "$$\n",
    "\n",
    "其中 ${e^{ - \\beta H}}$ 是矩阵 $ - \\beta H$ 的矩阵指数，$\\beta  = \\frac{1}{{kT}}$ 是系统的逆温度参数，其中 $T$ 是温度参数，$k$ 是玻尔兹曼常数 (这个例子中我们取 $k = 1$)。作为一个上手的例子，这里我们首先考虑一个3量子比特的哈密顿量及其吉布斯态。\n",
    "\n",
    "$$\n",
    "H = -Z \\otimes Z \\otimes I - I \\otimes Z \\otimes Z - Z \\otimes I \\otimes Z,\\,\n",
    "I=\\left [\n",
    "\\begin{matrix}\n",
    "1 & 0  \\\\\n",
    "0 & 1  \\\\\n",
    "\\end{matrix} \n",
    "\\right ],\n",
    "Z=\\left [\n",
    "\\begin{matrix}\n",
    "1 & 0  \\\\\n",
    "0 & -1  \\\\\n",
    "\\end{matrix} \n",
    "\\right ].$$\n",
    "\n",
    "\n",
    "这个例子中，我们将逆温度参数设置为 $\\beta = 1.5$。此外，为了方便测试结果，我们按照定义提前生成好了理想情况的吉布斯态。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "name": "#%% \n"
    }
   },
   "outputs": [],
   "source": [
    "N = 4        # 量子神经网络的宽度\n",
    "N_SYS_B = 3  # 用于生成吉布斯态的子系统B的量子比特数   \n",
    "SEED = 14    # 固定随机种子"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "beta = 1.5   # 设置逆温度参数 beta\n",
    "\n",
    "# 生成用泡利字符串表示的特定的哈密顿量\n",
    "H = [[-1.0, 'z0,z1'], [-1.0, 'z1,z2'], [-1.0, 'z0,z2']]\n",
    "\n",
    "# 生成哈密顿量的矩阵信息\n",
    "hamiltonian = pauli_str_to_matrix(H, N_SYS_B)\n",
    "\n",
    "# 生成理想情况下的目标吉布斯态 rho\n",
    "rho_G = scipy.linalg.expm(-1 * beta * hamiltonian) / np_trace(scipy.linalg.expm(-1 * beta * hamiltonian))\n",
    "\n",
    "# 设置成 Paddle quantum 所支持的数据类型\n",
    "hamiltonian = hamiltonian.astype(\"complex128\")\n",
    "rho_G = rho_G.astype(\"complex128\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 搭建量子神经网络\n",
    "\n",
    "- 在这个案例中，我们将通过训练量子神经网络QNN（也可以理解为参数化量子电路）来训练吉布斯态。这里，我们提供一个简单的4量子位的量子电路如下：\n",
    "![Ugibbs.jpg](https://release-data.cdn.bcebos.com/PIC%2FUgibbs.jpg)\n",
    "\n",
    "- 我们需要预设一些电路的参数，比如电路有4量子比特，其中第1个量子位是辅助系统，第2-4个量子位是用以产生吉布斯态的子系统。\n",
    "\n",
    "- 初始化其中的变量参数，${\\bf{\\theta }}$ 代表我们量子神经网络中的参数组成的向量。\n",
    "         "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "接下来我们根据上图中的电路设计，通过 Paddle Quantum 的 `UAnsatz` 函数和内置的 `real_entangled_layer(theta, D)` 电路模板来高效搭建量子神经网络。      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def U_theta(initial_state, theta, N, D):\n",
    "    \"\"\"\n",
    "    Quantum Neural Network\n",
    "    \"\"\"\n",
    "    \n",
    "    # 按照量子比特数量/网络宽度初始化量子神经网络\n",
    "    cir = UAnsatz(N)\n",
    "    \n",
    "    # 内置的 {R_y + CNOT} 电路模板\n",
    "    cir.real_entangled_layer(theta[:D], D)\n",
    "    \n",
    "    # 铺上最后一列 R_y 旋转门\n",
    "    for i in range(N):\n",
    "        cir.ry(theta=theta[D][i][0], which_qubit=i)\n",
    "        \n",
    "    # 量子神经网络作用在给定的初始态上\n",
    "    final_state = cir.run_density_matrix(initial_state)\n",
    "\n",
    "    return final_state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 配置训练模型 - 损失函数\n",
    "- 现在我们已经有了数据和量子神经网络的架构，我们将进一步定义合适的训练参数、模型和损失函数来达到我们的目标。\n",
    "- 具体的我们参考的是论文[4]中的方法，核心思想是利用吉布斯态达到了最小自由能的性质。\n",
    "- 通过作用量子神经网络 $U(\\theta)$ 在初始态上，我们可以得到输出态 $\\left| {\\psi \\left( {\\bf{\\theta }} \\right)} \\right\\rangle $, 其在第2-4个量子位的态记为 $\\rho_B(\\theta)$。\n",
    "- 设置训练模型中的的损失函数。在吉布斯态学习中，我们利用冯诺依曼熵函数的截断来进行自由能的估计，相应的损失函数参考[4]可以设为 \n",
    "$loss= {L_1} + {L_2} + {L_3}$，其中 ${L_1}= tr(H\\rho_B)$, ${L_2} = 2{\\beta^{-1}}{tr}(\\rho_B^2)$ , $L_3 = - {\\beta ^{ - 1}}\\big(tr(\\rho_B^3) + 3\\big)/2$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "class Net(fluid.dygraph.Layer):\n",
    "    \"\"\"\n",
    "    Construct the model net\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, shape, param_attr=fluid.initializer.Uniform(low=0.0, high=2*PI, seed=SEED),\n",
    "                 dtype='float64'):\n",
    "        super(Net, self).__init__()\n",
    "        \n",
    "        # 初始化 theta 参数列表，并用 [0, 2*pi] 的均匀分布来填充初始值\n",
    "        self.theta = self.create_parameter(shape=shape, attr=param_attr, dtype=dtype, is_bias=False)\n",
    "        \n",
    "        # 初始化 rho = |0..0><0..0| 的密度矩阵\n",
    "        self.initial_state = fluid.dygraph.to_variable(density_op(N))\n",
    "\n",
    "    # 定义损失函数和前向传播机制\n",
    "    def forward(self, H, N, N_SYS_B, D):\n",
    "\n",
    "        # 施加量子神经网络\n",
    "        rho_AB = U_theta(self.initial_state, self.theta, N, D)\n",
    "\n",
    "        # 计算偏迹 partial trace 来获得子系统B所处的量子态 rho_B\n",
    "        rho_B = partial_trace(rho_AB, 2 ** (N - N_SYS_B), 2 ** (N_SYS_B), 1)\n",
    "        \n",
    "        # 计算三个子损失函数\n",
    "        rho_B_squre = matmul(rho_B, rho_B)\n",
    "        loss1 = (trace(matmul(rho_B, H))).real\n",
    "        loss2 = (trace(rho_B_squre)).real * 2 / beta\n",
    "        loss3 = - ((trace(matmul(rho_B_squre, rho_B))).real + 3) / (2 * beta)\n",
    "        \n",
    "        # 最终的损失函数\n",
    "        loss = loss1 + loss2 + loss3  \n",
    "\n",
    "        return loss, rho_B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 配置训练模型 - 模型参数\n",
    "\n",
    "在进行量子神经网络的训练之前，我们还需要进行一些训练的超参数设置，主要是学习速率 (LR, learning rate)、迭代次数(ITR, iteration)和量子神经网络计算模块的深度 (D, Depth)。这里我们设定学习速率为0.5, 迭代次数为50次。读者不妨自行调整来直观感受下超参数调整对训练效果的影响。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "ITR = 50 # 设置训练的总迭代次数\n",
    "LR = 0.5 # 设置学习速率\n",
    "D = 1    # 设置量子神经网络中重复计算模块的深度 Depth"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 进行训练\n",
    "\n",
    "- 当训练模型的各项参数都设置完成后，我们将数据转化为Paddle动态图中的变量，进而进行量子神经网络的训练。\n",
    "- 训练过程中我们用的是[Adam Optimizer](https://www.paddlepaddle.org.cn/documentation/docs/zh/api_cn/optimizer_cn/AdagradOptimizer_cn.html)，也可以调用Paddle中提供的其他优化器。\n",
    "- 我们将训练过程中的结果依次输出。\n",
    "- 特别的我们依次输出了我们学习到的量子态$\\rho_B(\\theta)$与吉布斯态$\\rho_G$的保真度，保真度越高说明QNN输出的态越接近于吉布斯态。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "name": "#%% \n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter: 5 loss: -2.5615 fid: 0.8631\n",
      "iter: 10 loss: -3.1189 fid: 0.9504\n",
      "iter: 15 loss: -3.3290 fid: 0.9732\n",
      "iter: 20 loss: -3.3502 fid: 0.9846\n",
      "iter: 25 loss: -3.3446 fid: 0.9868\n",
      "iter: 30 loss: -3.3630 fid: 0.9873\n",
      "iter: 35 loss: -3.3937 fid: 0.9916\n",
      "iter: 40 loss: -3.4087 fid: 0.9948\n",
      "iter: 45 loss: -3.4108 fid: 0.9954\n",
      "iter: 50 loss: -3.4110 fid: 0.9953\n"
     ]
    }
   ],
   "source": [
    "# 初始化paddle动态图机制\n",
    "with fluid.dygraph.guard():\n",
    "    \n",
    "    # 我们需要将 Numpy array 转换成 Paddle 动态图模式中支持的 variable\n",
    "    H = fluid.dygraph.to_variable(hamiltonian)\n",
    "\n",
    "    # 确定网络的参数维度\n",
    "    net = Net(shape=[D + 1, N, 1])\n",
    "\n",
    "    # 一般来说，我们利用Adam优化器来获得相对好的收敛，当然你可以改成SGD或者是RMS prop.\n",
    "    opt = fluid.optimizer.AdamOptimizer(learning_rate=LR, parameter_list=net.parameters())\n",
    "\n",
    "     # 优化循环\n",
    "    for itr in range(1, ITR + 1):\n",
    "        \n",
    "        # 前向传播计算损失函数并返回生成的量子态 rho_B\n",
    "        loss, rho_B = net(H, N, N_SYS_B, D)\n",
    "        \n",
    "        # 在动态图机制下，反向传播极小化损失函数\n",
    "        loss.backward()\n",
    "        opt.minimize(loss)\n",
    "        net.clear_gradients()\n",
    "\n",
    "        # 转换成 Numpy array 用以计算量子态的保真度 F(rho_B, rho_G)\n",
    "        rho_B = rho_B.numpy()\n",
    "        fid = state_fidelity(rho_B, rho_G)\n",
    "        \n",
    "        # 打印训练结果\n",
    "        if itr % 5 == 0:\n",
    "            print('iter:', itr, 'loss:', '%.4f' % loss.numpy(), 'fid:', '%.4f' % fid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 总结\n",
    "根据上面训练得到的结果，通过大概50次迭代，我们就能达到高于99.5%保真度的高精度吉布斯态，高效并精确地完成了吉布斯态的制备。我们可以通过print函数来输出学习到的量子神经网络的参数和它的输出态。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 参考文献\n",
    "\n",
    "[1] [Kieferová, M. & Wiebe, N. Tomography and generative training with quantum Boltzmann machines. Phys. Rev. A 96, 062327 (2017).](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.96.062327)\n",
    "\n",
    "[2] [Brandao, F. G. S. L. & Svore, K. M. Quantum Speed-Ups for Solving Semidefinite Programs. in 2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS) 415–426 (IEEE, 2017). ](https://ieeexplore.ieee.org/abstract/document/8104077)\n",
    "\n",
    "[3] [Somma, R. D., Boixo, S., Barnum, H. & Knill, E. Quantum Simulations of Classical Annealing Processes. Phys. Rev. Lett. 101, 130504 (2008).](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.101.130504)\n",
    "\n",
    "[4] [Wang, Y., Li, G. & Wang, X. Variational quantum Gibbs state preparation with a truncated Taylor series. arXiv:2005.08797 (2020).](https://arxiv.org/pdf/2005.08797.pdf)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}