tutorials/quantum_simulation/QSP_and_QSVT_CN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 量子信号处理与量子奇异值变换\n",
    "\n",
    "*Copyright (c) 2022 Institute for Quantum Computing, Baidu Inc. All Rights Reserved.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 概览"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**量子信号处理**（quantum signal processing, QSP）是一个首次由 Guang Hao Low 和 Issac L. Chuang [[1]](https://quantum-journal.org/papers/q-2019-07-12-163/) 提出的量子算法，其使用反射（$R_X$）和旋转（$R_Z$）算子来模拟标量的多项式变换。对于满足特定条件的多项式（比如切比雪夫多项式），QSP 表明这些多项式变换可以通过单比特上的量子门模拟。使用名为“酉算子线性组合”（linear combination of unitaries, LCU）的技巧，我们只需增加两个额外辅助比特，就能模拟任意复多项式变换。简而言之，QSP可拟合任意能被泰勒级数逼近的函数。\n",
    "\n",
    "论文 [[2]](https://doi.org/10.1145/3313276.3316366)进一步发展了 QSP 的思想，并提出使用高维旋转算子的量子算法来模拟矩阵的多项式变换。由于矩阵的多项式变换本质上是它的奇异值的多项式变换，因此该算法被称为**量子奇异值变换**（quantum singular value transformation, QSVT）。另外，我们还介绍名为**单比特化**（qubitization）的技术，该技术可使用一个辅助比特就能模拟高维旋转算子，大幅减少了模拟高维度旋转算子的资源消耗，使得我们能够用量子电路实现矩阵的多项式变换。\n",
    "\n",
    "基于论文 [[2]](https://doi.org/10.1145/3313276.3316366)，该教程会向大家介绍 QSP、QSVT 和单比特化的概念及其在量桨上的实现。\n",
    "\n",
    "接下来，我们先介绍块编码（block-encoding）的概念。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以下是必要的 libraries 和 packages。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.polynomial.polynomial import Polynomial\n",
    "from numpy.polynomial import Chebyshev\n",
    "import paddle\n",
    "\n",
    "# 量桨的通用函数\n",
    "import paddle_quantum\n",
    "from paddle_quantum.ansatz import Circuit\n",
    "from paddle_quantum.qinfo import dagger\n",
    "from paddle_quantum.linalg import unitary_random_with_hermitian_block, density_matrix_random\n",
    "\n",
    "# 量桨的 QSVT 模块函数\n",
    "from paddle_quantum.qsvt import poly_matrix, ScalarQSP, Phi_verification, reflection_based_quantum_signal_processing, QSVT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 块编码"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "（量子）块编码是一种给矩阵数据编码的方式。在量子计算中，对于一个给定的矩阵 $A \\in \\mathbb{C}^{2^n \\times 2^n}$，我们能找到一个酉矩阵 $U \\in \\mathbb{C}^{2^m \\times 2^m}$ 满足如下等式\n",
    "\n",
    "$$\n",
    "W U V = \n",
    "\\begin{bmatrix}\n",
    "    A & 0 \\\\\n",
    "    0 & 0\n",
    "\\end{bmatrix}\n",
    "= A \\oplus 0I^{\\otimes (m - n)}, \\tag{1}\n",
    "$$\n",
    "\n",
    "其中，$W, V \\in \\mathbb{C}^{2^m \\times 2^m}$ 是两个正交投影算子，这样的酉矩阵 $U$ 被称为 $A$ 的**块编码**（block encoding）。\n",
    "\n",
    "一般地，投影算子通常是 $W = V = |0\\rangle^{\\otimes (m - n)}\\langle0|^{\\otimes (m - n)} \\otimes I^{\\otimes n}$，那么酉矩阵 $U$ 在标准正交基下可表示为\n",
    "\n",
    "$$\n",
    "U = \\begin{bmatrix}\n",
    "    A & ... \\\\\n",
    "    ... & ...\n",
    "\\end{bmatrix}, \\tag{2}\n",
    "$$\n",
    "\n",
    "即 $A$ 位于 $U$ 的左上角。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 编码与解码"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "一般来说，块编码的表达形式并非是唯一的。比如，当 $A$ 可以对角化时，$U$ 可以被构建为\n",
    "\n",
    "$$\n",
    "U = \\begin{bmatrix}\n",
    "    A & i\\sqrt{I^{\\otimes n} - A^2} \\\\\n",
    "    i\\sqrt{I^{\\otimes n} - A^2} & A\n",
    "\\end{bmatrix}; \\tag{3}\n",
    "$$\n",
    "\n",
    "当 $A$ 是酉矩阵时，我们则可以选择 $U = A \\oplus I^{\\otimes (m - n)}$，即 $U$ 是受控的 $A$。特别地，文献[[3]](http://arxiv.org/abs/2203.10236) [[4]](http://arxiv.org/abs/2206.03505)给出了一些使用量子电路实现块编码的方案。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_qubits = 3\n",
    "num_block_qubits = 2\n",
    "\n",
    "# 创建一个 2 比特大小的厄密矩阵 A 的 3 比特大小的块编码 U\n",
    "# create a 3-qubit block encoding unitary U of a random 2-qubit Hermitian matrix A\n",
    "U = unitary_random_with_hermitian_block(num_qubits)\n",
    "A = U[0: 2 ** num_block_qubits, 0: 2 ** num_block_qubits]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "另一方面，我们也可以从 矩阵 $A \\in \\mathbb{C}^{2^n \\times 2^n}$ 的块编码 $U$ 中解码出 $A$ 的信息。用$|0\\rangle^{\\otimes (m - n)}\\langle0|^{\\otimes (m - n)} \\otimes \\rho$表示输入量子态，块编码 $U$ 作为量子电路，那么输出量子态为\n",
    "\n",
    "$$\n",
    "U (|0\\rangle^{\\otimes (m - n)}\\langle0|^{\\otimes (m - n)} \\otimes \\rho) U^\\dagger\n",
    "= \\begin{bmatrix}\n",
    "    A \\rho A^\\dagger & ... \\\\\n",
    "    ... & ...\n",
    "\\end{bmatrix}\n",
    "= |0\\rangle^{\\otimes (m - n)}\\langle0|^{\\otimes (m - n)} \\otimes A \\rho A^\\dagger + (I^{\\otimes (m - n)} - |0\\rangle^{\\otimes (m - n)}\\langle0|^{\\otimes (m - n)}) \\otimes ... \\tag{4}\n",
    "$$\n",
    "\n",
    "对 $|0\\rangle^{\\otimes (m - n)}$ 所在的寄存器作测量，若结果为 $0^{\\otimes (m - n)}$，那么我们成功地提取了 $A$ 的信息\n",
    "<!-- 测量第一个量子寄存器。若结果为 $0^{\\otimes (m - n)}$，则信息提取成功。 -->\n",
    "\n",
    "![block-decoding](figures/QSVT-fig-block-decoding.png \"图 1：块解码的量子实现\")\n",
    "\n",
    "解码的成功概率取决于测量到 $0^{\\otimes (m - n)}$ 的几率，其会随着 $m - n$ 的增长而指数性下降。然而，若我们可以控制 $m - n$ 的大小，那么这就不再是一个问题。在本教程中我们会假设 $m = n + 1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "qubits [0] collapse to the state |0> with probability 0.24548491835594177\n"
     ]
    }
   ],
   "source": [
    "# 设置密度矩阵后端\n",
    "paddle_quantum.set_backend('density_matrix')\n",
    "\n",
    "# 定义输入态\n",
    "rho = density_matrix_random(num_block_qubits)\n",
    "zero_state = paddle.eye(2 ** (num_qubits - num_block_qubits), 1) # 创建 0 态的矢量\n",
    "input_state = paddle_quantum.State(paddle.kron(zero_state @ dagger(zero_state), rho))\n",
    "\n",
    "# 定义辅助寄存器\n",
    "aux_register = list(range(num_qubits - num_block_qubits))\n",
    "\n",
    "# 创建线路\n",
    "cir = Circuit(num_qubits)\n",
    "cir.oracle(U, list(range(num_qubits)))\n",
    "cir.collapse(aux_register, desired_result='0', if_print = True) # 调用塌缩算子\n",
    "\n",
    "# 获取输出态及输出 rho\n",
    "output_state = cir(input_state)\n",
    "output_rho = output_state.data[0: 2 ** num_block_qubits, 0: 2 ** num_block_qubits]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们可以验证 ``output_rho`` 与 $\\frac{A \\rho A^\\dagger}{\\text{Tr}(A \\rho A^\\dagger)}$ 相同。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "期望 rho 和 输出 rho 的误差是 0.0\n"
     ]
    }
   ],
   "source": [
    "expect_rho = A @ rho @ dagger(A)\n",
    "expect_rho /= paddle.trace(expect_rho)\n",
    "print(f\"期望 rho 和 输出 rho 的误差是 {paddle.norm(paddle.abs(expect_rho - output_rho)).item()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "当 $m = 1$（因此 $n = 0$）时，$A$ 是一个标量。在这种情况下我们则可以算出 $U$ 对于量子态 $|0\\rangle$ 的期望值来获取 $A$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 量子信号处理"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "论文 [[2]]((https://doi.org/10.1145/3313276.3316366)) 中的 Theorem 3-4 证明，若一个次数为 $k$ 的复数多项式 $P \\in \\mathbb{C}[x]$ 满足\n",
    "\n",
    "- 若 $k$ 为偶/奇数，则 $P$ 是一个偶/奇多项式；\n",
    "- 对于所有的 $x \\in [-1, 1]$，$|P(x)| \\leq 1$;\n",
    "- 对于所有的 $x \\in (-\\infty, -1] \\cup [1, \\infty)$，$|P(x)| \\geq 1$;\n",
    "- 若 $k$ 为偶，则 $P(ix)P^*(ix) \\geq 1$,\n",
    "\n",
    "则存在多项式 $Q \\in \\mathbb{C}[x]$ 和角度矢量 $\\Phi = (\\phi_0, ..., \\phi_k) \\in \\mathbb{R}^{k + 1}$ 使得对于所有的 $x \\in [-1, 1]$,\n",
    "\n",
    "$$\n",
    "W_\\Phi(x) := R_Z(-2\\phi_0) \\prod_{j = 1}^k R_X(\\arccos(x)) R_Z(-2\\phi_j)\n",
    "= \\begin{bmatrix}\n",
    "   P(x) & i Q(x)\\sqrt{1 - x^2} \\\\\n",
    "   iQ^*(x)\\sqrt{1 - x^2} & P^*(x) \n",
    "\\end{bmatrix}. \\tag{5}\n",
    "$$\n",
    "\n",
    "换句话说，我们可以使用 $k + 1$ 个旋转（$R_Z$）门和 $k$ 个反射 ($R_X$) 门去获得一个 $P(x)$ 的块编码酉矩阵 $W_\\Phi(x)$。块解码后 $P(x)$ 的模拟就完成了。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以下是一些满足上述条件的例子。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 简单的奇多项式\n",
    "# P = Polynomial([0, 0, 0, 0, 0, 1])\n",
    "# P = Polynomial([0, 0.5, 0, 0, 0, 0.5])\n",
    "P = Polynomial([0, 1 / 3, 0, 0, 0, 1 / 3, 0, 1 / 3])\n",
    "\n",
    "# 次数为10的第一类切比雪夫多项式\n",
    "# P = Chebyshev([0 for _ in range(10)] + [1]).convert(kind=Polynomial)\n",
    "\n",
    "# 次数为11的第一类切比雪夫多项式\n",
    "# P = Chebyshev([0 for _ in range(11)] + [1]).convert(kind=Polynomial)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "值得注意的是，$W_\\Phi$ 可以模拟切比雪夫多项式变换。对于次数为 $k$ 的第一类切比雪夫多项式，其对应的角度矢量为 $(0, \\pi, ..., \\pi)$ （如 $k$ 为偶数）或 $(\\pi, ..., \\pi)$ (如 $k$ 为奇数)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "量桨的 QSVT 模块有一个内置类 `ScalarQSP` 用于完成量子信号处理，并可以调用函数 `Phi_verification` 来验证是否 $W_\\Phi$ 是 $P$ 的块编码。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tensor(shape=[8], dtype=float64, place=Place(cpu), stop_gradient=True,\n",
      "       [ 3.14159265,  1.39460474, -0.44313783,  1.09757975, -1.09757975,\n",
      "         0.44313783, -1.39460474,  3.14159265])\n"
     ]
    }
   ],
   "source": [
    "qsp = ScalarQSP(P)\n",
    "print(qsp.Phi)\n",
    "assert Phi_verification(qsp.Phi.numpy(), P)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们也可以通过量子线路来验证 $\\Phi$ 的正确性。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--Rz(-6.28)----Rx(-2.73)----Rz(2.789)----Rx(-2.73)----Rz(-0.88)----Rx(-2.73)----Rz(2.195)----Rx(-2.73)----Rz(-2.19)----Rx(-2.73)----Rz(0.886)----Rx(-2.73)----Rz(-2.78)----Rx(-2.73)----Rz(-6.28)--\n",
      "                                                                                                                                                                                                   \n",
      "模拟 P(x) 的误差是 1.7893723625681406e-08\n"
     ]
    }
   ],
   "source": [
    "x = np.random.rand() * 2 - 1 # 随机在 [-1, 1] 间选取 x\n",
    "cir = qsp.block_encoding(x)\n",
    "print(cir)\n",
    "print(f\"模拟 P(x) 的误差是 {np.abs(cir.unitary_matrix()[0, 0].item() - P(x))}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### QSP 变种"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "另外，参数数量可以被进一步减少至 $k$，通过找到一个角度矢量 $\\Phi \\in \\mathbb{R}^k$ 使得\n",
    "\n",
    "$$\n",
    "R_\\Phi(x):=\\prod_{i = 1}^{k} e^{i \\phi_i \\sigma_z} R(x)\n",
    "= \\begin{bmatrix}\n",
    "    P(x) & ... \\\\\n",
    "    ... & ...\n",
    "\\end{bmatrix}, \\tag{6}\n",
    "$$\n",
    "\n",
    "这里\n",
    "\n",
    "$$\n",
    "R(x) = \\begin{bmatrix}\n",
    "    x & \\sqrt{1 - x^2} \\\\\n",
    "    \\sqrt{1 - x^2} & -x\n",
    "\\end{bmatrix}. \\tag{7}\n",
    "$$\n",
    "\n",
    "在量桨里我们通过调用 `reflection_based_quantum_signal_processing` 来计算 $\\Phi$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[15.70796327 -0.17619159 -2.01393416 -0.47321657 -2.66837608 -1.12765849\n",
      " -2.96540107]\n"
     ]
    }
   ],
   "source": [
    "print(reflection_based_quantum_signal_processing(P))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> 注: $R_{-\\Phi}(x)$ 是 $P^*(x)$ 的块编码。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 量子奇异值变换"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来，我们介绍量子奇异值变换的概念。\n",
    "\n",
    "对于 $\\mathcal{X} \\in \\mathbb{C}^{2^m \\times 2^m}$ 为一个矩阵，当它可以被一个酉矩阵 $U$ 和两个正交投影算子 $ W, V$ 表示时，即 $\\mathcal{X} = W U V$，那么矩阵 $\\mathcal{X}$ 有如下表达式：\n",
    "\n",
    "$$\n",
    "\\mathcal{X} = \\sum_{j=1}^{2^m} \\lambda_j |\\omega_j\\rangle \\langle\\nu_j|, \\tag{8}\n",
    "$$\n",
    "\n",
    "其中 $\\{|\\omega_j\\rangle\\}_{j=1}^{2^m}$ 和 $\\{|\\nu_j\\rangle\\}_{j=1}^{2^m}$ 是由 $W$ 和 $V$ 的本征态扩展生成的标准正交基，且 $\\lambda_j$ 是 $\\mathcal{X}$ 的奇异值，\n",
    "\n",
    "$$\n",
    "\\lambda_j = \\begin{cases}\n",
    "\\langle\\omega_j| U |\\nu_j\\rangle \\text{ 如 } j \\leq \\text{rank}(\\mathcal{X}) \\\\\n",
    "0 \\text{  其他情况 }\n",
    "\\end{cases} \\tag{9}\n",
    "$$\n",
    "\n",
    "对于给定的周期性函数 $f$，矩阵$\\mathcal{X}$ 的**奇异值变换**(SVT) 定义为\n",
    "\n",
    "$$\n",
    "f^{(SV)}(\\mathcal{X}) := \\sum_{j=1}^{2^m} f(\\lambda_j)\n",
    "\\begin{cases}\n",
    "       |\\nu_j\\rangle\\langle\\nu_j| & \\text{  如 } f \\text{ 是偶函数} \\\\\n",
    "       |\\omega_j\\rangle\\langle\\nu_j| & \\text{  如 } f \\text{ 是奇函数}\n",
    "\\end{cases}. \\tag{10}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 分解与 QSP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "论文 [[2]]((https://doi.org/10.1145/3313276.3316366)) 证明在**合适**的基下, 与 $\\lambda_j$ 相关的子空间可以将 $U, V$ 和 $W$ 分解为\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "U & = \\bigoplus_{j: \\lambda_j \\in [0, 1)} R(\\lambda_j) \\oplus\n",
    "        \\bigoplus_{j: \\lambda_j = 1} [1] \\oplus\n",
    "        \\bigoplus_{j: V|\\nu_j\\rangle = |\\nu_j\\rangle,\\, WU|\\nu_j\\rangle = 0} [1] \\oplus\n",
    "        \\bigoplus_{j: W|\\omega_j\\rangle = |\\omega_j\\rangle,\\, VU^\\dagger|\\omega_j\\rangle = 0} [1] \\oplus \n",
    "        \\bigoplus_{j: V|\\nu_j\\rangle = W|\\omega_j\\rangle = 0} [...], \\\\\n",
    "V & = \\bigoplus_{j: \\lambda_j \\in [0, 1)} |0\\rangle\\langle0| \\oplus\n",
    "        \\bigoplus_{j: \\lambda_j = 1} [1] \\oplus\n",
    "        \\bigoplus_{j: V|\\nu_j\\rangle = |\\nu_j\\rangle,\\, WU|\\nu_j\\rangle = 0} [1] \\oplus\n",
    "        \\bigoplus_{j: W|\\omega_j\\rangle = |\\omega_j\\rangle,\\, VU^\\dagger|\\omega_j\\rangle = 0} [0]  \\oplus \n",
    "        \\bigoplus_{j: V|\\nu_j\\rangle = W|\\omega_j\\rangle = 0} [0] \\text{ and}\\\\\n",
    "W & = \\bigoplus_{j: \\lambda_j \\in [0, 1)} |0\\rangle\\langle0| \\oplus\n",
    "        \\bigoplus_{j: \\lambda_j = 1} [1] \\oplus\n",
    "        \\bigoplus_{j: V|\\nu_j\\rangle = |\\nu_j\\rangle,\\, WU|\\nu_j\\rangle = 0} [0] \\oplus\n",
    "        \\bigoplus_{j: W|\\omega_j\\rangle = |\\omega_j\\rangle,\\, VU^\\dagger|\\omega_j\\rangle = 0} [1]  \\oplus \n",
    "        \\bigoplus_{j: V|\\nu_j\\rangle = W|\\omega_j\\rangle = 0} [0] .\\\\\n",
    "\\end{split} \\tag{11}\n",
    "$$\n",
    "\n",
    "$f^{(SV)}$ 也可以被分解为\n",
    "\n",
    "$$\n",
    "f^{(SV)}(\\mathcal{X}) = \\sum_{j: \\lambda_j \\in [0, 1)} f(\\lambda_j)\\,... + \\sum_{j: \\lambda_j = 1} f(1)\\,... + \\sum_{j: \\lambda_j = 0, V |\\nu_j\\rangle \\neq 0} f(0)\\,... + \\sum_{j: \\lambda_j = 0, W |\\omega_j\\rangle \\neq 0} f(0)\\,... + \\sum_{j: \\lambda_j = 0, \\text{otherwise}} f(0)\\,...\\,, \\tag{12}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在我们将上述分解与量子信号处理联系到一起。\n",
    "\n",
    "假设函数 $P$ 是一个次数为 $k$ 并满足上节所述条件，则存在 $\\Phi \\in \\mathbb{R}^k$ 使得 $R_\\Phi$ 是 $P$ 的块编码，并且 $P$ 满足 $P(1) = e^{i \\sum_{j=1}^k \\phi_j}$ 和 $P(0) = e^{i \\sum_{j=1}^k (-1)^{j+1} \\phi_j}$。\n",
    "\n",
    "定义 $U_\\Phi$ 为\n",
    "\n",
    "$$\n",
    "U_\\Phi :=  \\begin{cases}\n",
    "& \\prod_{j = 1}^{k / 2} e^{i\\phi_{2j - 1} (2V - I)} U^\\dagger e^{i\\phi_{2j} (2W - I)} U \\text{  如 } k \\text{ 是偶数} \\\\\n",
    "e^{i\\phi_1 (2W - I)} U & \\prod_{j = 1}^{(k - 1) / 2} e^{i\\phi_{2j} (2V - I)} U^\\dagger e^{i\\phi_{2j+1} (2W - I)} U \\text{  如 } k \\text{ 是奇数}\n",
    "\\end{cases}. \\tag{13}\n",
    "$$\n",
    "\n",
    "那么\n",
    "\n",
    "$$\n",
    "U_\\Phi =  \\begin{cases}\n",
    "        \\bigoplus R_\\Phi(\\lambda_j) \\oplus\n",
    "        \\bigoplus [e^{i \\sum_{j=1}^k \\phi_j}] \\oplus\n",
    "        \\bigoplus [e^{i \\sum_{j=1}^k (-1)^{j+1} \\phi_j}] \\oplus\n",
    "        \\bigoplus [e^{i \\sum_{j=1}^k (-1)^{j+1} \\phi_j}] \\oplus \n",
    "        \\bigoplus [...]\n",
    "\\text{  如 } k \\text{ 是偶数} \\\\\n",
    "        \\bigoplus R_\\Phi(\\lambda_j) \\oplus\n",
    "        \\bigoplus [e^{i \\sum_{j=1}^k \\phi_j}] \\oplus\n",
    "        \\bigoplus [e^{i \\sum_{j=1}^k (-1)^{j+1} \\phi_j}] \\oplus\n",
    "        \\bigoplus [e^{i \\sum_{j=1}^k (-1)^{j+1} \\phi_j}] \\oplus \n",
    "        \\bigoplus [...]\n",
    "\\text{  如 } k \\text{ 是奇数}\n",
    "\\end{cases}, \\tag{14}\n",
    "$$\n",
    "\n",
    "因此\n",
    "\n",
    "$$\n",
    "P^{(SV)}(\\mathcal{X}) = \\begin{cases} \n",
    "        \\bigoplus \\begin{bmatrix}\n",
    "                P(\\lambda_j) & 0 \\\\\n",
    "                0 & 0 \\\\\n",
    "        \\end{bmatrix} \\oplus\n",
    "        \\bigoplus [P(1)] \\oplus\n",
    "        \\bigoplus [P(0)] \\oplus\n",
    "        \\bigoplus [0] \\oplus \n",
    "        \\bigoplus [0] = V U_\\Phi V\n",
    "\\text{  如 } k \\text{ 是偶数} \\\\ \n",
    "        \\bigoplus \\begin{bmatrix}\n",
    "                P(\\lambda_j) & 0 \\\\\n",
    "                0 & 0 \\\\\n",
    "        \\end{bmatrix} \\oplus\n",
    "        \\bigoplus [P(1)] \\oplus\n",
    "        \\bigoplus [0] \\oplus\n",
    "        \\bigoplus [0] \\oplus \n",
    "        \\bigoplus [0] = W U_\\Phi V \n",
    "\\text{  如 } k \\text{ 是奇数}\n",
    "\\end{cases}. \\tag{15}\n",
    "$$\n",
    "\n",
    "若是我们可以使用量子算法实现 $V U_\\Phi V$ 或者 $W U_\\Phi V$，那么这样的变换就是 $\\mathcal{X}$ 的量子奇异值变换。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### QSVT 和块编码"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "假设相对于正交投影算子 $W$ 和 $V$， $U$ 是一个 $A$ 的块编码酉矩阵。那么 $\\mathcal{X} = A \\oplus 0I^{\\otimes (m - n)}$ 且 $P^{(SV)}(\\mathcal{X}) = P^{(SV)}(A) \\oplus 0I^{\\otimes (m - n)}$. 因此，$U_\\Phi$ 是 $P^{(SV)}(A)$ 的块编码.\n",
    "\n",
    "当 $W = V$ 时，$\\{|\\omega_j\\rangle \\}_{j=1}^{2^m} = \\{ |\\nu_j\\rangle \\}_{j=1}^{2^m}$。记 $P(x) := \\sum_{i=0}^k c_i x^i$，我们就可以找到\n",
    "\n",
    "$$\n",
    "P^{(SV)}(\\mathcal{X}) = \\sum_{j=1}^{2^m} P(\\lambda_j) |\\nu_j\\rangle \\langle\\nu_j| = P(X) = P(A) \\oplus 0I^{\\otimes (m - n)}. \\tag{16}\n",
    "$$\n",
    "\n",
    "换句话说，当 $W = V$ 时，QSVT 将一个 $A$ 的块编码转为了一个 $P(A)$ 的块编码。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "当 $W = V = |0\\rangle^{\\otimes (m - n)}\\langle0|^{\\otimes (m - n)} \\otimes I^{\\otimes n}$ 时，量桨有一个内置类 `QSVT` 来完成量子奇异值变换。我们可以调用其成员函数 `QSVT.block_encoding_matrix()` 来验证上述理论的正确性。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "qsvt = QSVT(P, U, num_block_qubits)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 找到 P(A) 及\n",
    "# find P(A) and its expected eigenvalues, note that they are computed in different ways\n",
    "expect_PA = poly_matrix(P, A).numpy()\n",
    "expect_eig_result = np.sort(list(map(lambda x: P(x), np.linalg.eigvals(A.numpy()))))\n",
    "\n",
    "# Calculate U_\\Phi and extract eigenvalues of block encoded matrix\n",
    "U_Phi = qsvt.block_encoding_matrix().numpy()\n",
    "actual_PA = U_Phi[0:2 ** num_block_qubits, 0:2 ** num_block_qubits]\n",
    "actual_eig_result = np.sort(np.linalg.eigvals(actual_PA))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "模拟 P(X) 的误差\n",
      "     最大绝对值,  1.450928305889582e-07\n",
      "     百分比,  3.370838085602552e-07\n"
     ]
    }
   ],
   "source": [
    "print(\"模拟 P(X) 的误差\")\n",
    "print(\"     最大绝对值, \", np.amax(abs(expect_PA - actual_PA)))\n",
    "print(\"     百分比, \", np.linalg.norm(expect_PA - actual_PA) / np.linalg.norm(expect_PA))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "模拟 P(X) 的本征值误差\n",
      "     最大绝对值,  1.8375562631798232e-07\n",
      "     百分比,  2.776073800450616e-07\n"
     ]
    }
   ],
   "source": [
    "print(\"模拟 P(X) 的本征值误差\")\n",
    "print(\"     最大绝对值, \", np.amax(abs(expect_eig_result - actual_eig_result)))\n",
    "print(\"     百分比, \", np.linalg.norm(expect_eig_result - actual_eig_result) / np.linalg.norm(expect_eig_result))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 单比特化: $U_\\Phi$ 的量子实现"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在的问题是，如何用量子线路去实现 $U_\\Phi$ ？通常我们认为块编码 $U$ 已经通过量子电路实现，那么剩下来只需实现 $e^{i\\phi (2W - I)}$ 和 $e^{i\\phi (2V - I)}$。值得注意的是，这两个算子本质上是分别在 $W$ 和 $V$ 的投影空间下施加的 $R_Z(-2\\phi)$ 算子。\n",
    "\n",
    "论文 [[1]](https://quantum-journal.org/papers/q-2019-07-12-163/) 中的 Lemma 10 认为，应该将投影算子的映射空间投射到一个辅助比特上再进行旋转，那么通过纠缠的性质这个旋转也会同时应用于主寄存器上。\n",
    "\n",
    "这些理论被总结至下图：\n",
    "\n",
    "![U_Phi](figures/QSVT-fig-U_Phi.png \"图 2: QSVT 的量子实现，这里 k 是 P 的多项式次数\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "当 $W = V = |0\\rangle^{\\otimes (m - n)}\\langle0|^{\\otimes (m - n)} \\otimes I^{\\otimes n}$ 时，量桨可以调用成员函数 `QSVT.block_encoding_circuit` 来创建一个图 2 所示的量子线路。我们可以通过比较该线路的输出态和量子态 $(U_\\Phi \\otimes I)|\\psi\\rangle|0\\rangle$（$|\\psi\\rangle \\in \\mathbb{C}^{2^m}$）来证明所建立的线路是可以模拟 $U_\\Phi$ 的。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 设置态矢量后端\n",
    "paddle_quantum.set_backend('state_vector')\n",
    "\n",
    "# 随即量子态满足最后一个量子比特固定为 0 态\n",
    "ket_0 = [1.0, 0.0]\n",
    "psi = np.array([np.random.rand() + 1j * np.random.rand() for _ in range(2 ** num_qubits)])\n",
    "psi = np.kron(psi / np.linalg.norm(psi), ket_0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "期望量子态和实际量子态的差距是 3.0184617834444857e-07\n"
     ]
    }
   ],
   "source": [
    "expect_state = np.kron(U_Phi, np.eye(2)) @ psi\n",
    "\n",
    "cir = qsvt.block_encoding_circuit()\n",
    "actual_state = cir(paddle_quantum.State(psi, dtype='complex64')).data.numpy()\n",
    "\n",
    "print(f\"期望量子态和实际量子态的差距是 {np.linalg.norm(expect_state - actual_state)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> 注：如果我们在图 2 中的辅助寄存器两边加上 Hadamard 门，那么该量子线路就会转为模拟其多项式的实数变化。再结合 LCU 的技巧，理论上我们只需要额外 1 个量子比特来量子模拟任意范数小于 1 的复数多项式。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 应用: 振幅放大"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "设 $|\\psi\\rangle$ 为 $m$ 比特大小的量子态并满足 $|\\psi\\rangle = \\sin(\\theta)|\\psi_{\\text{good}}\\rangle + \\cos(\\theta) |\\psi_{\\text{bad}}\\rangle$。**振幅放大**是一个用于增大 $|\\psi_{\\text{good}}\\rangle$ 至 1 的量子算法。该算法可以从另一个角度来解释。设定 $U$ 为一个酉矩阵且 $W$ 和 $V$ 为两个正交投影算子，使得 $|\\psi\\rangle = U V |0\\rangle^{\\otimes m}$ 且 $\\sin(\\theta)|\\psi_{\\text{good}}\\rangle = W |\\psi\\rangle$，意味着 $\\sin(\\theta)|\\psi_{\\text{good}}\\rangle = W U V |0\\rangle^{\\otimes m}$。因此，振幅放大本质上是 $\\mathcal{X} = W U V$ 的奇异值变换。\n",
    "在这一章节我们会根据论文 [[2]]((https://doi.org/10.1145/3313276.3316366)) 的 Theorem 28，来展示如何用 QSVT 来做定点（即 $\\theta$ 满足 $\\theta = \\frac{\\pi}{2k}$，$k \\in \\mathbb{Z}$）振幅放大算法。\n",
    "\n",
    "假设 $\\sin(\\theta)|\\psi_{\\text{good}}\\rangle = (|0\\rangle^{\\otimes (m - n)}\\langle0|^{\\otimes (m - n)} \\otimes I^{\\otimes n}) |\\psi\\rangle$ 且 $\\theta = \\frac{\\pi}{2k}$（$k \\in \\mathbb{Z}$）。那么就能推导出 $W = V = |0\\rangle^{\\otimes (m - n)}\\langle0|^{\\otimes (m - n)} \\otimes I^{\\otimes n}$ 且 $\\mathcal{X} = A \\oplus 0I^{\\otimes (m - n)}$，这里 $A$ 位于 $U$ 的左上角。 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可见 $\\mathcal{X} |0\\rangle^{\\otimes m} = \\sin(\\frac{\\pi}{2k}) |\\psi_\\text{good}\\rangle$。因此，我们旨在找到一个酉矩阵为 $U_\\Phi$ 的量子线路使得 $W U_\\Phi V = \\sin^{-1}(\\theta) \\mathcal{X}$， 从而达到 $W U_\\Phi V |0\\rangle^{\\otimes m} = |\\psi_{\\text{good}}\\rangle$ 的目的。因为 $\\mathcal{X} = W U V$，$\\mathcal{X}$ 的奇异值的绝对值是 $\\sin(\\frac{\\pi}{2k})$。我们进一步选择 $P(x) = (-1)^k T_k (x)$，这里 $T_k$ 是次数为 $k$ 的第一类切比雪夫多项式。最后通过\n",
    "\n",
    "$$\n",
    "P(\\sin(\\frac{\\pi}{2k})) = (-1)^k T_k (\\sin(\\frac{\\pi}{2k})) = (-1)^k \\cos(\\frac{k - 1}{2}\\pi) = 1, \\tag{17}\n",
    "$$\n",
    "\n",
    "我们可得，对于多项式 $P$ 的 $\\mathcal{X}$ 的量子奇异值变换是 $B := \\frac{1}{\\sin(\\frac{\\pi}{2k})} A$ 的块编码。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "设定 $k = 3$ 使得 $\\sin(\\frac{\\pi}{2k}) = \\frac{1}{2}$，同时随机选择酉矩阵 $U$。其左上角为 $A$，我们需要证明 $U$ 的 QSVT 的左上角为 $B = 2A$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "P = -1 * Chebyshev([0 for _ in range(3)] + [1]).convert(kind=Polynomial)\n",
    "U = unitary_random_with_hermitian_block(num_qubits, is_unitary=True)\n",
    "A = U[0:2 ** num_block_qubits, 0:2 ** num_block_qubits].numpy()\n",
    "\n",
    "amplifier = QSVT(P, U, num_block_qubits)\n",
    "U_Phi = amplifier.block_encoding_unitary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "QSVT 的精确值为 3.1799856969882967e-07\n"
     ]
    }
   ],
   "source": [
    "B = U_Phi[0:2 ** num_block_qubits, 0:2 ** num_block_qubits].numpy()\n",
    "print(f\"QSVT 的精确值为 {np.linalg.norm(B - 2 * A)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由上可见，我们成功将 $A$ 的大小翻倍。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "## 参考文献\n",
    "\n",
    "[1] Low, Guang Hao, and Isaac L. Chuang. \"Hamiltonian simulation by qubitization.\" [Quantum 3 (2019): 163.](https://doi.org/10.22331/q-2019-07-12-163)\n",
    "\n",
    "[2] Gilyén, András, et al. \"Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics.\" [Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing. 2019.](https://doi.org/10.1145/3313276.3316366)\n",
    "\n",
    "[3] Camps, Daan, et al. \"Explicit Quantum Circuits for Block Encodings of Certain Sparse Matrices.\" [arXiv preprint arXiv:2203.10236 (2022).](http://arxiv.org/abs/2203.10236)\n",
    "\n",
    "[4] Clader, B. David, et al. \"Quantum Resources Required to Block-Encode a Matrix of Classical Data.\" [arXiv preprint arXiv:2206.03505 (2022).](http://arxiv.org/abs/2206.03505)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.13 ('pq')",
   "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.8.13"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "08942b1340a5932ff3a93f52933a99b0e263568f3aace1d262ffa4d9a0f2da31"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}