{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# HHL 算法\n", "\n", "[![下载Notebook](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/website-images/master/resource/_static/logo_notebook.svg)](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/notebook/master/mindquantum/zh_cn/case_library/mindspore_hhl_algorithm.ipynb)\n", "[![下载样例代码](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/website-images/master/resource/_static/logo_download_code.svg)](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/notebook/master/mindquantum/zh_cn/case_library/mindspore_hhl_algorithm.py)\n", "[![查看源文件](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/website-images/master/resource/_static/logo_source.svg)](https://gitee.com/mindspore/docs/blob/master/docs/mindquantum/docs/source_zh_cn/case_library/hhl_algorithm.ipynb)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 概述\n", "\n", "**HHL算法的目的:** 给定一个 Hermitian $A$ 和一个单位向量 $\\vec{b}$,\n", "求解方程 $A \\vec{x} = \\vec{b}$ 。\n", "\n", "令 $\\left| b \\right\\rangle = \\sum_{j=1}^{N} b_{j } \\left| j \\right\\rangle$,\n", "算法的关键在于在 $\\left|b \\right\\rangle$ 上模拟 $e^{i A \\Delta t }$ 。\n", "\n", "设 $A$ 的谱分解如下\n", "\n", "$$A = \\sum_{j=1}^{N} \\lambda_{j } \\left| u_{j} \\right\\rangle \\left\\langle u_{j} \\right|$$\n", "\n", "那么\n", "\n", "$$e^{i A \\Delta t } = \\sum_{j=1}^{N} e^{i \\lambda_{j} \\Delta t } \\left| u_{j} \\right\\rangle \\left\\langle u_{j} \\right|$$\n", "\n", "将 $\\left| b \\right\\rangle$ 用基 $\\left\\lbrace \\left| u_{j} \\right\\rangle \\right\\rbrace$ 表示 : $\\left| b \\right\\rangle = \\sum_{j=1}^{N} \\beta_{j} \\left| u_{j} \\right\\rangle$ 。\n", "\n", "使用量子相位估计 QPE(Quantum Phase Estimation),其中 $U = e^{i A \\Delta t}$,作用结果如下:\n", "\n", "$$\\left| b \\right\\rangle \\xrightarrow{\\text{QPE}} \\sum_{j=1}^{N} \\beta_{j} \\left| u_{j} \\right\\rangle | \\widetilde{ \\lambda_j } \\rangle$$\n", "\n", "其中 $| \\widetilde{\\lambda_j} \\rangle$ 是特征值 $\\lambda_{j}$ 的估计(波浪线表示估计值,下文同理)。\n", "\n", "很容易看出方程的解就是\n", "\n", "$$\n", "\\left| x \\right\\rangle = A^{-1 } \\left| b \\right\\rangle = \\sum_{j}^{ } \\beta_{j } (\\lambda_{j})^{-1 } \\left| u_{j} \\right\\rangle\n", "$$\n", "\n", "所以我们只需要将信息从量子态 $| \\widetilde{\\lambda_j} \\rangle$ 中提取出来即可。" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## HHL 算法步骤\n", "\n", "下面详细介绍 HHL 算法的详细步骤,包含详细的数学推导。" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### 数据预处理\n", "\n", "前面提到 $A$ 需要是 Hermitian 的,也就是 $A^\\dagger = A$。其实这个条件可以放宽,如果 $A$ 不是 Hermitian 的,构造 $\\tilde{A}$ 如下:\n", "\n", "$$\n", "\\tilde{A} = \\begin{pmatrix}\n", "0 & A \\\\\n", "A^{\\dagger } & 0 \\\\\n", "\\end{pmatrix}\n", "$$\n", "\n", "然后求解方程\n", "\n", "$$\n", "\\tilde{A} \\vec{y} = \\begin{pmatrix}\n", "\\vec{b} \\\\\n", "0 \\\\\n", "\\end{pmatrix}\n", "$$\n", "\n", "很容易验证方程的解 $\\vec{y}$ 一定具有如下的形式\n", "\n", "$$\n", "\\vec{y} = \\begin{pmatrix}\n", "0 \\\\\n", "\\vec{x} \\\\\n", "\\end{pmatrix}\n", "$$\n", "\n", "其中 $A \\vec{x} = \\vec{b}$ 。\n", "\n", "对于 $\\vec{b}$,因为 HHL 是量子算法,所以需要输入的是量子态 $|b\\rangle = \\sum_{j=1}^N b_j |j\\rangle$ 。关于如何制备这样量子态,如何将经典信息编码到量子信息,不是本算法关心的重点,这里只假设制备好了这样的量子态。" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### 制备算子 $U = e^{i A \\Delta t}$\n", "\n", "如何高效的制备算子 $U = e^{i A \\Delta t}$ 是核心问题。后面将会看到,整个 HHL 算法的复杂度主要取决于算子 $U$ ,所以需要高效的模拟 $U$ 。\n", "\n", "制备这个时间演化算子 $U = e^{i A \\Delta t}$ 属于量子系统模拟(哈密顿量模拟)问题,本身是一个重要的量子问题,这里并不展开讨论。\n", "\n", "在原论文中,对矩阵 $A$ 作出了**稀疏性**的限制,其目的就是能够高效的模拟 $e^{i A \\Delta t}$。对于一般的稠密矩阵 $A$,模拟其演化复杂度可能很高。" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "这里不去展开证明 QPE 的正确性,读者只需要知道 QPE 的输入是一个酉算子 $U$ 和其特征向量 $|u\\rangle$,设 $U |u\\rangle = e^{2 \\pi i \\varphi} |u\\rangle$ ,输出是 $\\varphi$ 的估计 $\\tilde{\\varphi}$ 。\n", "\n", "具体来说,使用 $t$ 个辅助量子比特, QPE 的作用如下\n", "\n", "$$\\left| b \\right\\rangle \\left| 0 \\right\\rangle^{\\otimes t} \\left| 0 \\right\\rangle = \\sum_{j=1}^{N} \\beta_{j} \\left| u_{j} \\right\\rangle \\left| 0 \\right\\rangle^{\\otimes t} \\left| 0 \\right\\rangle\n", "\\xrightarrow{\\text{QPE}} \\sum_{j=1}^{N} \\beta_{j} \\left| u_{j} \\right\\rangle \\left| \\widetilde{\\varphi_j} \\right\\rangle \\left| 0 \\right\\rangle$$\n", "\n", "其中 $\\left| \\widetilde{\\varphi_j} \\right\\rangle$ 是相位 $\\varphi_{j}$ 的估计。\n", "这里在 $t$ 个辅助比特后面还增加了一个辅助比特,是用于后面的旋转步骤。\n", "\n", "为了更好的理解,这里举一个例子:\n", "取\n", "\n", "$$A = Z = \\left| 0 \\right\\rangle \\left\\langle 0 \\right| - \\left| 1 \\right\\rangle \\left\\langle 1 \\right|$$\n", "\n", "那么\n", "\n", "$$U = e^{i A \\Delta t } = e^{i \\Delta t} \\left| 0 \\right\\rangle \\left\\langle 0 \\right| + e^{-i \\Delta t} \\left| 1 \\right\\rangle \\left\\langle 1 \\right|$$\n", "\n", "两个相位分别是 $\\varphi_{1} = \\Delta t / 2 \\pi$ 和 $\\varphi_{2} = -\\Delta t / 2 \\pi$ 。\n", "\n", "如果我们使用 $t = 4$ 个辅助比特,并且取 $\\Delta t = \\pi / 4$,那么 $\\varphi_{1} = 1 / 8$ 和 $\\varphi_{2} = - 1 / 8$ 。\n", "\n", "两个相位的估计值分别是 $\\widetilde{\\varphi_1} = 0.0010 \\times 2^4 = 2$ 和 $\\widetilde{\\varphi_2} = 0.1110 \\times 2^4 = 14$ 。\n", "因为相位是 $[0, 1)$ 之间的小数 $0.q_{t_1}q_{t_2}\\ldots q_{t_t}$,\n", "乘上 $2^t$ 得到其对应的量子态 $|q_{t_1} \\ldots q_{t_t} \\rangle$ 。\n", "需要指出的是,因为相位是模 $1$ 的小数,$-0.0010$ 被映射到了 $0.1110$,我们想要还原相位(为了还原 $\\lambda$),首先通过选取足够小的 $\\Delta t$ 使得 $|\\varphi| < 1 / 2$,这样如果 $|q_{t_1} \\ldots q_{t_t}\\rangle < 2^{t-1}$ 对应正数,$|q_{t_1}\\ldots q_{t_t}\\rangle > 2^{t-1}$ 对应负数。\n", "\n", "通过简单的定量计算,我们得到如下的映射关系:\n", "\n", "$$\n", "\\frac{\\lambda \\Delta t}{2 \\pi} = \\varphi = \\begin{cases}\n", "\\tilde{\\varphi} / 2^{t} & , \\varphi > 0 & , \\tilde{\\varphi} < 2^{t-1} \\\\\n", "\\tilde{\\varphi} / 2^{t} - 1 & , \\varphi < 0 & , \\tilde{\\varphi} > 2^{t-1}\\\\\n", "\\end{cases}\n", "$$\n", "\n", "其中 $\\lambda$ 是 $A$ 的特征值,$\\varphi$ 是 $U = e^{i A \\Delta t}$ 的相位,$\\tilde{\\varphi}$ 是 $\\varphi$ 的估计。" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### 条件旋转\n", "\n", "经过 QPE 后,整个量子态如下\n", "\n", "$$\n", "\\sum_{j=1}^{N} \\beta_{j} \\left| u_{j} \\right\\rangle \\left| \\widetilde{\\varphi_j} \\right\\rangle \\left| 0 \\right\\rangle\n", "$$\n", "\n", "想要将 $\\lambda_j$ 的信息从量子态 $\\left| \\widetilde{\\varphi_j} \\right\\rangle$ 中提取出来,我们需要使用条件旋转门 $CR(k)$,其作用效果如下\n", "\n", "$$\n", "CR(k) \\left| \\tilde{\\varphi} \\right\\rangle \\left| b \\right\\rangle = \\begin{cases}\n", "\\left| \\tilde{\\varphi} \\right\\rangle \\left| b \\right\\rangle & , k \\neq \\tilde{\\varphi} \\\\\n", "\\left| \\tilde{\\varphi} \\right\\rangle R_y \\left( 2 \\arcsin \\frac{C}{\\lambda} \\right) \\left| b \\right\\rangle & , k = \\tilde{\\varphi} \\\\\n", "\\end{cases}\n", "$$\n", "\n", "简单来说,只有当 $k$ 选择了正确的 $\\tilde{\\varphi}$,才会对后面的量子比特作用旋转操作。\n", "\n", "因为我们不知道正确的 $\\tilde{\\varphi}$ 是什么,所以暴力的枚举所有的可能 $CR(k)$, $k = 1, \\ldots , 2^{t}-1$ 。\n", "\n", "旋转的效果是很简单的\n", "\n", "$$\n", "\\prod_{k=1}^{2^{t} - 1} I\\otimes CR(k) \\sum_{j=1}^{N} \\beta_{j} \\left| u_{j} \\right\\rangle \\left| \\widetilde{\\varphi_j} \\right\\rangle \\left| 0 \\right\\rangle\n", "= \\sum_{j=1}^{N} \\beta_{j} \\left| u_{j} \\right\\rangle \\left| \\widetilde{\\varphi_j} \\right\\rangle \\left( \\sqrt{1 - \\left( \\frac{C}{\\lambda_{j}} \\right)^{2}} \\left| 0 \\right\\rangle + \\frac{C}{\\lambda_{j}} \\left| 1 \\right\\rangle \\right)\n", "$$\n", "\n", "再作用一次逆 QPE,整体的量子态如下\n", "\n", "$$\n", "\\left| \\psi \\right\\rangle = \\sum_{j=1}^{N} \\beta_{j} \\left| u_{j} \\right\\rangle \\left| 0 \\right\\rangle^{\\otimes t} \\left( \\sqrt{1 - \\left( \\frac{C}{\\lambda_{j}} \\right)^{2}} \\left| 0 \\right\\rangle + \\frac{C}{\\lambda_{j}} \\left| 1 \\right\\rangle \\right)\n", "$$" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### 测量\n", "\n", "对最后一个量子比特进行测量,当测量结果是 $\\left| 1 \\right\\rangle$ 时,其概率是\n", "\n", "$$\n", "p_1 = C^{2} \\sum_{j=1}^{N} ( \\beta_{j} / \\lambda_{j} )^{2}\n", "$$\n", "\n", "测量后的量子态变为\n", "\n", "$$\n", "\\left| \\psi_{1} \\right\\rangle = \\frac{1}{\\left( \\sum_{j=1}^{N} \\left( \\beta_{j} / \\lambda_{j} \\right)^{2} \\right)^{-1}} \\sum_{j=1}^{N} \\frac{\\beta_{j}}{\\lambda_{j}} \\left| u_{j} \\right\\rangle \\left| 0 \\right\\rangle^{\\otimes t} \\left| 1 \\right\\rangle\n", "$$\n", "\n", "显然此时三个量子寄存器之间不存在纠缠,如果只看第一个量子寄存器,它已经是 $|x\\rangle$ 的状态了(忽略掉一个归一化系数)。\n", "\n", "需要指出的是,旋转操作里面有一个参数 $C$,我们看到 $C$ 不会影响最后结果的正确性,但是却实实在在的影响了得到结果的概率 $p_{1}$ 。\n", "我们肯定希望 $C$ 能够尽可能大,从而使 $p_1$ 尽可能大。但是 $C \\leq \\left\\vert \\lambda_{j} \\right\\vert$ ,它要比所有特征值的绝对值还要小,如果没有先验信息,不知道绝对值最小的 $\\lambda$ 有多大,那么就只能保守的取一个很小的 $C$,然后可以通过振幅放大技术来增大得到结果的概率。" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## MindQuantum 实现\n", "\n", "这里使用 MindQuantum 实现一个简单的例子,仅说明 HHL 算法的过程和正确性。\n", "\n", "为了计算方便,我们选取一个简单的\n", "$A = Z = \\begin{bmatrix}\n", "1 & 0 \\\\\n", "0 & -1\n", "\\end{bmatrix}$,因为其时间演化算子 $e^{i Z \\Delta t} = R_z(- 2 \\Delta t)$ ,比较容易实现。" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "eigenvalues of A:\n", " [ 1. -1.]\n", "eigenvectors of A:\n", " [[1. 0.]\n", " [0. 1.]]\n", "b: [0.6 0.8]\n", "Solution of Ax=b is: [ 0.6 -0.8]\n" ] } ], "source": [ "import numpy as np\n", "\n", "A = np.array([[1, 0], [0, -1]])\n", "es, vs = np.linalg.eig(A)\n", "\n", "print(f\"eigenvalues of A:\\n {es}\")\n", "print(f\"eigenvectors of A:\\n {vs}\")\n", "\n", "b = np.array([0.6, 0.8])\n", "print(f\"b: {b}\")\n", "print(f\"Solution of Ax=b is: {np.linalg.solve(A, b)}\")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "这里导入所需要的各种函数:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from mindquantum.core.gates import H, X, RY, RZ, Measure, Power, BARRIER\n", "from mindquantum.core.circuit import Circuit\n", "from mindquantum.simulator import Simulator" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "对于量子态 $|b\\rangle = \\cos\\theta |0\\rangle + \\sin\\theta |1\\rangle$ 的制备,可以通过一个 $R_y(2 \\theta)$ 实现。\n", "\n", "下面的代码制备了 $|b\\rangle = 0.6 |0\\rangle + 0.8 |1\\rangle$ 并进行了测量。" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "Shots:\n", " 10000 Keys: q0 0.0 0.128 0.255 0.383 0.511 0.639 0 3613 1 6387 probability " ], "text/plain": [ "" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "circ = Circuit()\n", "circ += RY(2 * np.arcsin(0.8)).on(0)\n", "circ += Measure().on(0)\n", "\n", "sim = Simulator(backend=\"mqvector\", n_qubits=1)\n", "res = sim.sampling(circ, shots=10000)\n", "res.svg()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "下面一步步构建出完整的量子线路,取 $t = 4$,$\\Delta t = \\pi / 4$, $C = 0.5$。\n", "\n", "- 首先是 QPE\n", "- 然后是 $CR(k)$,$k = 1, \\ldots, 15$\n", "- 接着是 逆QPE\n", "- 最后测量" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "q0: q1: q2: q3: q4: q5: RY 1.8546 H H H H RZ(-π/2)^1 RZ(-π/2)^2 RZ(-π/2)^4 RZ(-π/2)^8 H PS -π/2 H PS -π/4 PS -π/2 H PS -π/8 PS -π/4 PS -π/2 H X X X RY π X X X X X X RY 1.0472 X X X X X RY 0.6797 X X X X X RY 0.5054 X X X X X RY 0.4027 X X X X RY 0.3349 X X X RY 0.2867 X X X X RY 0.2507 X X X X X RY -0.2867 X X X X RY -0.3349 X X X RY -0.4027 X X X RY -0.5054 X X X RY -0.6797 X X RY -1.0472 X RY H PS π/2 PS π/4 PS π/8 H PS π/2 PS π/4 H PS π/2 H RZ(π/2)^1 RZ(π/2)^2 RZ(π/2)^4 RZ(π/2)^8 H H H H " ], "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from mindquantum.algorithm.library import qft\n", "\n", "# q1 ~ q4 : for QPE\n", "t = 4\n", "qc = [4, 3, 2, 1]\n", "\n", "# store |b> , and store the result |x>\n", "qb = 5\n", "\n", "# use for conditional rotation\n", "qs = 0\n", "\n", "# choose a time evolution\n", "dt = np.pi / 4\n", "\n", "# choose a small C\n", "C = 0.5\n", "\n", "circ = Circuit()\n", "\n", "# prepare b\n", "circ += RY(2 * np.arcsin(0.8)).on(qb)\n", "\n", "# QPE\n", "for i in qc:\n", " circ += H.on(i)\n", "\n", "for (i, q) in enumerate(qc):\n", " circ += Power(RZ(-2 * dt), 2**i).on(qb, q)\n", "\n", "# apply inverse QFT\n", "circ += qft(list(reversed(qc))).hermitian()\n", "\n", "# conditional rotate\n", "circ += BARRIER\n", "for k in range(1, 2**t):\n", " for i in range(t):\n", " if (k & (1 << i)) == 0:\n", " circ += X.on(qc[i])\n", " phi = k / (2**t)\n", " if k > 2**(t-1):\n", " phi -= 1.0\n", " l = 2 * np.pi / dt * phi\n", " circ += RY(2 * np.arcsin(C / l)).on(qs, qc)\n", "\n", " for i in range(t):\n", " if (k & (1 << i)) == 0:\n", " circ += X.on(qc[i])\n", " circ += BARRIER\n", "\n", "# apply inverse QPE\n", "circ += qft(list(reversed(qc)))\n", "\n", "for (i, q) in enumerate(qc):\n", " circ += Power(RZ(2 * dt), 2**i).on(qb, q)\n", "\n", "for i in qc:\n", " circ += H.on(i)\n", "\n", "circ.svg()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "如何验证我们的结果的正确性?测量。因为结果是量子态 $|x\\rangle$,目前 MindQuantum 没有相关的办法能够单独取出某一个量子比特的内部值,我们只能通过测量。" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "Shots:\n", " 100000 Keys: q5 q0 0.0 0.096 0.192 0.288 0.385 0.481 00 26790 01 9057 10 48066 11 16087 probability " ], "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sim = Simulator(backend=\"mqvector\", n_qubits=2 + t)\n", "sim.apply_circuit(circ)\n", "\n", "circ_m = Circuit()\n", "circ_m += Measure().on(qs)\n", "circ_m += Measure().on(qb)\n", "\n", "res = sim.sampling(circ_m, shots=100000)\n", "res.svg()" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.3602052179446389" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "res.data.get(\"01\", 0) / (res.data.get(\"01\", 0) + res.data.get(\"11\", 0))" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "通过统计当 `q0` 是 1 时,`q5` 的测量结果,我们看到 `q5 = 0` 的占比是 `0.36`,和答案的预期一致。\n", "\n", "当然这种方法只能得到振幅大小,想要得到更多信息,例如相位,需要对测量进行调整,这里不去讨论。" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", " \n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", " \n", " \n", "\n", "\n", "
SoftwareVersion
mindquantum0.9.11
scipy1.10.1
numpy1.23.5
SystemInfo
Python3.9.16
OSLinux x86_64
Memory8.3 GB
CPU Max Thread8
DateFri Dec 29 23:43:57 2023
\n" ], "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from mindquantum.utils.show_info import InfoTable\n", "\n", "InfoTable('mindquantum', 'scipy', 'numpy')" ] } ], "metadata": { "kernelspec": { "display_name": "base", "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.9.16" }, "orig_nbformat": 4 }, "nbformat": 4, "nbformat_minor": 2 }