tutorials/quantum_simulation/QSP_and_QSVT_EN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quantum Signal Processing and Quantum Singular Value Transformation\n",
    "*Copyright (c) 2022 Institute for Quantum Computing, Baidu Inc. All Rights Reserved.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Quantum signal processing** (QSP), originally purposed by Guang Hao Low and Issac L. Chuang [[1]](https://quantum-journal.org/papers/q-2019-07-12-163/), is an algorithm for simulation of polynomial transformation of scalars using reflection ($R_X$) and rotation ($R_Z$) operators. QSP states that, for a specific group of polynomials like the Chebyshev polynomials of the first kind, such simulation can be done by a single qubit. Moreover, using a technique called \"linear combination of unitaries\", we can use two more ancilla qubit to simulate all complex polynomials and hence approximate any functions that can be expanded by Taylor series.\n",
    "\n",
    "The idea of QSP is further generalized by paper [[2]](https://doi.org/10.1145/3313276.3316366), so that it can simulate polynomial transformation of matrices using high-dimensional rotation operators. Since the polynomial transformation of a diagonalizable matrix is intrinsically the polynomial transformation of its singular values, this algorithm is called the **quantum singular value transformation** (QSVT).\n",
    "\n",
    "Moreover, a technique called **qubitization** is further employed to substantially reduce the complexity of high-dimensional rotation operators, so that one (ancilla) qubit suffices to perform the rotation in larger space.\n",
    "\n",
    "Based on paper [[2]](https://doi.org/10.1145/3313276.3316366), this tutorial will illustrate quantum signal processing, quantum singular value transformation and qubitization, and how to implement these algorithms in PaddleQuantum. But first of all, we need to present the idea of block encoding."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are some necessary libraries and 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",
    "# PaddleQuantum related packages\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",
    "from paddle_quantum.qsvt import poly_matrix, ScalarQSP, Phi_verification, reflection_based_quantum_signal_processing, QSVT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Block Encoding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Block encoding is a method to encode matrices of data. In quantum computing, for a matrix $A \\in \\mathbb{C}^{2^n \\times 2^n}$, a unitary $U \\in \\mathbb{C}^{2^m \\times 2^m}$ is said to be a **block encoding** of $A$ if there exists two orthogonal projectors $W, V \\in \\mathbb{C}^{2^m \\times 2^m}$ such that\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",
    "This could appear more familiar when $W = V = |0\\rangle^{\\otimes (m - n)}\\langle0|^{\\otimes (m - n)} \\otimes I^{\\otimes n}$, where above equation implies\n",
    "\n",
    "$$\n",
    "U = \\begin{bmatrix}\n",
    "    A & ... \\\\\n",
    "    ... & ...\n",
    "\\end{bmatrix}. \\tag{2}\n",
    "$$\n",
    "\n",
    "That is, $A$ is the left-upper block of matrix $U$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Encoding and Decoding of Matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are numerous ways to find a block encoding unitary $U$ of a matrix $A$. For example, when $A$ is diagonalizable, $U$ can be constructed as\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",
    "when $A$ is a unitary, we can simply select $U = A \\oplus I^{\\otimes (m - n)}$ i.e. $U$ is a controlled $A$. Despite the fact that we do not hold a mature method for doing this in a quantum system, several studies [[3]](http://arxiv.org/abs/2203.10236) [[4]](http://arxiv.org/abs/2206.03505) have been conducted on the quantum implementation of block encodings. In this tutorial, we will assume that there exists an efficient algorithm to block encode a matrix (or at least a Hermitian matrix) in a quantum circuit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Realization in PaddleQuantum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_qubits = 3\n",
    "num_block_qubits = 2\n",
    "\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": [
    "As for the decoding part, the information of matrix $A \\in \\mathbb{C}^{2^n \\times 2^n}$ can be extracted by applying its block encoding unitary $U \\in \\mathbb{C}^{2^m \\times 2^m}$ on $|0\\rangle^{\\otimes (m - n)}\\langle0|^{\\otimes (m - n)} \\otimes \\rho$, where $\\rho$ is an $n$-qubit quantum state. Then the output state is\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",
    "Measuring the first quantum register with outcome $0^{\\otimes (m - n)}$ gives the desired information.\n",
    "\n",
    "![block-decoding](figures/QSVT-fig-block-decoding.png \"Figure 1: Graph Illustration of Block Decoding\")\n",
    "\n",
    "The success probability of decoding is the probability of measuring $0^{\\otimes (m - n)}$, which can be exponentially small as $m - n$ increases. However, this would no longer be a problem if we can control the size of first register. In this tutorial we will assume $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.25072506070137024\n"
     ]
    }
   ],
   "source": [
    "# set density matrix backend\n",
    "paddle_quantum.set_backend('density_matrix')\n",
    "\n",
    "# define input state\n",
    "rho = density_matrix_random(num_block_qubits)\n",
    "zero_state = paddle.eye(2 ** (num_qubits - num_block_qubits), 1) # create a 0 state (in state_vector form)\n",
    "input_state = paddle_quantum.State(paddle.kron(zero_state @ dagger(zero_state), rho))\n",
    "\n",
    "# define auxiliary register\n",
    "aux_register = list(range(num_qubits - num_block_qubits))\n",
    "\n",
    "# construct the circuit\n",
    "cir = Circuit(num_qubits)\n",
    "cir.oracle(U, list(range(num_qubits)))\n",
    "cir.collapse(aux_register, desired_result='0', if_print = True) # call Collapse operator\n",
    "\n",
    "# get output_state and actual output 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": [
    "We can verify that the ``output_rho`` is the same as $\\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": [
      "the difference between input and output is 0.0\n"
     ]
    }
   ],
   "source": [
    "expect_rho = A @ rho @ dagger(A)\n",
    "expect_rho /= paddle.trace(expect_rho)\n",
    "print(f\"the difference between input and output is {paddle.norm(paddle.abs(expect_rho - output_rho)).item()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When $m = 1$ (and hence $n = 0$), $A$ is a scalar. In this case we can calculate the expectance of $U$ in terms of state $|0\\rangle$ to receive $A$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantum Signal Processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Theorem 3-4 in paper [[2]]((https://doi.org/10.1145/3313276.3316366)) state that, for a polynomial $P \\in \\mathbb{C}[x]$ with $\\deg(P) = k$ that satisfies\n",
    "\n",
    "- if $k$ is even/odd, then $P$ is a even/odd polynomial;\n",
    "- $|P(x)| \\leq 1, \\,\\forall x \\in [-1, 1]$;\n",
    "- $|P(x)| \\geq 1, \\,\\forall x \\in (-\\infty, -1] \\cup [1, \\infty)$;\n",
    "- if $k$ is even, then $P(ix)P^*(ix) \\geq 1$,\n",
    "\n",
    "there exists a polynomial $Q \\in \\mathbb{C}[x]$ and a vector of angles $\\Phi = (\\phi_0, ..., \\phi_k) \\in \\mathbb{R}^{k + 1}$ such that $\\forall 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",
    "That is, we can use $k + 1$ rotation ($R_Z$) gates and $k$ reflection ($R_X$) gates to receive a block encoding unitary $W_\\Phi(x)$ of $P(x)$. After block decoding we have completed the simulation of $P(x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are some polynomials that satisfy above requirements. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# simple odd polynomials\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",
    "# Chebyshev polynomials of first kind with degree 10\n",
    "# P = Chebyshev([0 for _ in range(10)] + [1]).convert(kind=Polynomial)\n",
    "\n",
    "# Chebyshev polynomials of first kind with degree 11\n",
    "# P = Chebyshev([0 for _ in range(11)] + [1]).convert(kind=Polynomial)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that such structure of $W_\\Phi$ naturally fits the family of Chebyshev polynomials. Indeed, for a Chebyshev polynomials of first kind with order $k$, its corresponding $\\Phi$ is $(0, \\pi, ..., \\pi)$ (if $k$ is even) or $(\\pi, ..., \\pi)$ (if $k$ is odd). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PaddleQuantum has a built-in class `ScalarQSP` for quantum signal processing and a function `Phi_verification` to verify whether $W_\\Phi$ is a block encoding of $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": [
    "We can also verify the correctness of $\\Phi$ from the corresponding circuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--Rz(-6.28)----Rx(-1.74)----Rz(2.789)----Rx(-1.74)----Rz(-0.88)----Rx(-1.74)----Rz(2.195)----Rx(-1.74)----Rz(-2.19)----Rx(-1.74)----Rz(0.886)----Rx(-1.74)----Rz(-2.78)----Rx(-1.74)----Rz(-6.28)--\n",
      "                                                                                                                                                                                                   \n",
      "the error of simulating P(x) is 6.378721218542117e-08\n"
     ]
    }
   ],
   "source": [
    "x = np.random.rand() * 2 - 1 # random x in [-1, 1]\n",
    "cir = qsp.block_encoding(x)\n",
    "print(cir)\n",
    "print(f\"the error of simulating P(x) is {np.abs(cir.unitary_matrix()[0, 0].item() - P(x))}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Variant of QSP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Moreover, the number of parameters can be further simplified to $k$ by finding a vector of angles $\\Phi \\in \\mathbb{R}^k$ such that\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",
    "where\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",
    "In PaddleQuantum, we can compute such $\\Phi$ by function `reflection_based_quantum_signal_processing`."
   ]
  },
  {
   "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": [
    "> Note: $R_{-\\Phi}(x)$ is a block encoding of $P^*(x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantum Singular Value Transformation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To start with, we need to define the singular value transformation. \n",
    "\n",
    "Let $\\mathcal{X} \\in \\mathbb{C}^{2^m \\times 2^m}$ be a matrix. Suppose there exists a unitary $U$ and two orthogonal projectors $W$ and $V$ such that $\\mathcal{X} = W U V$. Then there exists two orthonormal basis $\\{|\\omega_j\\rangle\\}_{j=1}^{2^m}$ and $\\{|\\nu_j\\rangle\\}_{j=1}^{2^m}$ extended by eigenstates of $W$ and $V$ respectively, such that\n",
    "\n",
    "$$\n",
    "\\mathcal{X} = \\sum_{j=1}^{2^m} \\lambda_j |\\omega_j\\rangle \\langle\\nu_j|, \\tag{8}\n",
    "$$\n",
    "\n",
    "where \n",
    "\n",
    "$$\n",
    "\\lambda_j = \\begin{cases}\n",
    "\\langle\\omega_j| U |\\nu_j\\rangle \\text{  if } j \\leq \\text{rank}(\\mathcal{X}) \\\\\n",
    "0 \\text{  otherwise }\n",
    "\\end{cases} \\tag{9}\n",
    "$$\n",
    "\n",
    "is a singular value of $\\mathcal{X}$. The **singular value transformation** (SVT) of $\\mathcal{X}$ for a function $f$ is defined to be\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{  if } f \\text{ is even} \\\\\n",
    "       |\\omega_j\\rangle\\langle\\nu_j| & \\text{  if } f \\text{ is odd}\n",
    "\\end{cases}. \\tag{10}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Decomposition and QSP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Paper [[2]]((https://doi.org/10.1145/3313276.3316366)) proves that under **appropriate** choices of basis, $U, V$ and $W$ can be decomposed by subspaces of $\\lambda_j$ such that\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",
    "We can decompose $f^{(SV)}$ in a similar manner.\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": [
    "Let's connect above decomposition with quantum signal processing. Suppose function $f$ is a polynomial $P \\in \\mathbb{C}[x]$ with degree $k$ that satisfies the requirements in the last section. \n",
    "Then there exists $\\Phi \\in \\mathbb{R}^k$ such that $R_\\Phi$ is a block encoding of $P$. Moreover, $P$ satisfies $P(1) = e^{i \\sum_{j=1}^k \\phi_j}$ and $P(0) = e^{i \\sum_{j=1}^k (-1)^{j+1} \\phi_j}$.\n",
    "Now Define $U_\\Phi$ as\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{  if } k \\text{ is even} \\\\\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{  if } k \\text{ is odd}\n",
    "\\end{cases}. \\tag{13}\n",
    "$$\n",
    "\n",
    "Then\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{  if } k \\text{ is even} \\\\\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{  if } k \\text{ is odd}\n",
    "\\end{cases} \\tag{14}\n",
    "$$\n",
    "\n",
    "and hence\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{  if } k \\text{ is even} \\\\ \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{  if } k \\text{ is odd}\n",
    "\\end{cases}. \\tag{15}\n",
    "$$\n",
    "\n",
    "If we can realize $V U_\\Phi V$ or $W U_\\Phi V$ by a quantum algorithm, then such transformation is the quantum singular value transformation of $\\mathcal{X}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### QSVT and Block Encoding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose $U$ is a block encoding unitary of $A$ with respect to orthogonal projectors $W$ and $V$. Then $\\mathcal{X} = A \\oplus 0I^{\\otimes (m - n)}$ and $P^{(SV)}(\\mathcal{X}) = P^{(SV)}(A) \\oplus 0I^{\\otimes (m - n)}$. Therefore, $U_\\Phi$ is a block encoding of $P^{(SV)}(A)$.\n",
    "\n",
    "When $W = V$, $\\{|\\omega_j\\rangle \\}_{j=1}^{2^m} = \\{ |\\nu_j\\rangle \\}_{j=1}^{2^m}$. Denote $P(x) := \\sum_{i=0}^k c_i x^i$. We can find that \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",
    "That is, when $W = V$, QSVT maps a block encoding of $A$ to a block encoding of $P(A)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PaddleQuantum has a built-in class ``QSVT`` for quantum singular value transformation when $W = V = |0\\rangle^{\\otimes (m - n)}\\langle0|^{\\otimes (m - n)} \\otimes I^{\\otimes n}$. We can call its member function ``QSVT.block_encoding_matrix()`` to verify the correctness of above theories, from entries to entries and from eigenvalues to eigenvalues."
   ]
  },
  {
   "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": [
    "# find P(A) and its expected eigenvalues, note that they are computed in different ways\n",
    "expect_PX = 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_PX = U_Phi[0:2 ** num_block_qubits, 0:2 ** num_block_qubits]\n",
    "actual_eig_result = np.sort(np.linalg.eigvals(actual_PX))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "error of simulating P(X)\n",
      "     maximum absolute,  1.6997650495437753e-07\n",
      "     percentage,  3.4201093195057237e-07\n"
     ]
    }
   ],
   "source": [
    "print(\"error of simulating P(X)\")\n",
    "print(\"     maximum absolute, \", np.amax(abs(expect_PX - actual_PX)))\n",
    "print(\"     percentage, \", np.linalg.norm(expect_PX - actual_PX) / np.linalg.norm(expect_PX))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "error of eigenvalues of simulating P(X)\n",
      "     maximum absolute,  2.3308557146996258e-07\n",
      "     percentage,  2.2962108806419593e-07\n"
     ]
    }
   ],
   "source": [
    "print(\"error of eigenvalues of simulating P(X)\")\n",
    "print(\"     maximum absolute, \", np.amax(abs(expect_eig_result - actual_eig_result)))\n",
    "print(\"     percentage, \", np.linalg.norm(expect_eig_result - actual_eig_result) / np.linalg.norm(expect_eig_result))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Qubitization: Quantum Realization of $U_\\Phi$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One question is, how to realize $U_\\Phi$ by a quantum circuit? Note that we assume $U$ is implementable. Then the remaining problem is the realization of unitaries $e^{i\\phi (2W - I)}$ and $e^{i\\phi (2V - I)}$, which are essentially $R_Z(-2\\phi)$ operators performed in projection spaces of $W$ and $V$ respectively.\n",
    "\n",
    "To achieve this, Lemma 10 in paper [[1]](https://quantum-journal.org/papers/q-2019-07-12-163/) projects the image space of projector into an ancilla qubit and perform the rotation operation on that qubit. Then by entanglement such rotation is simultaneously applied to the main register. \n",
    "\n",
    "The results are summarized to the following circuit:\n",
    "\n",
    "![U_Phi](figures/QSVT-fig-U_Phi.png \"Figure 2: Quantum Circuit for QSVT, where k is the degree of P\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When $W = V = |0\\rangle^{\\otimes (m - n)}\\langle0|^{\\otimes (m - n)} \\otimes I^{\\otimes n}$, PaddleQuantum can create a circuit in Figure 2 by member function ``QSVT.block_encoding_circuit``. We can show that the constructed quantum circuit can simulate $U_\\Phi$, by comparison of the output state and $(U_\\Phi \\otimes I)|\\psi\\rangle|0\\rangle$ for $|\\psi\\rangle \\in \\mathbb{C}^{2^m}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set state vector backend\n",
    "paddle_quantum.set_backend('state_vector')\n",
    "\n",
    "# arbitrary state such that last qubit is in state |0\\rangle\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": [
      "the different between expected and actual state is 2.383485634049416e-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\"the different between expected and actual state is {np.linalg.norm(expect_state - actual_state)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note: if we add Hadamard gates at the beginning and the end of the ancilla register in Figure 2, then the quantum circuit turns to simulate the real part of the polynomial. Combined with the technique of linear combination of unitaries, theocratically we can simulate all complex polynomials with norm smaller than 1 using quantum circuit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application: Amplitude Amplification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $|\\psi\\rangle$ be a $m$-qubit quantum state such that $|\\psi\\rangle = \\sin(\\theta)|\\psi_{\\text{good}}\\rangle + \\cos(\\theta) |\\psi_{\\text{bad}}\\rangle$. **Amplitude amplification** is a quantum algorithm that can increase the amplitude of $|\\psi_{\\text{good}}\\rangle$ to approximately 1.\n",
    "This algorithm can be viewed in a different prospective. Let $U$ be the unitary and $W, V$ be two orthogonal projectors such that $|\\psi\\rangle = U V |0\\rangle^{\\otimes m}$ and $\\sin(\\theta)|\\psi_{\\text{good}}\\rangle = W |\\psi\\rangle$, implying $\\sin(\\theta)|\\psi_{\\text{good}}\\rangle = W U V |0\\rangle^{\\otimes m}$. Therefore, amplitude amplification is essentially a singular value transformation of $\\mathcal{X} = W U V$. In this section we will show how to use QSVT to do fixed-point (i.e. $\\theta = \\frac{\\pi}{2k}$ for some $k \\in \\mathbb{Z}$) amplitude amplification from Theorem 28 of paper [[2]]((https://doi.org/10.1145/3313276.3316366)). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose $\\sin(\\theta)|\\psi_{\\text{bad}}\\rangle = (|\\psi_{\\text{good}}\\rangle^{\\otimes (m - n)}\\langle0|^{\\otimes (m - n)} \\otimes I^{\\otimes n}) |\\psi\\rangle$ and $\\theta = \\frac{\\pi}{2k}$ for some $k \\in \\mathbb{Z}$. Then $W = V = |0\\rangle^{\\otimes (m - n)}\\langle0|^{\\otimes (m - n)} \\otimes I^{\\otimes n}$ and $\\mathcal{X} = A \\oplus 0I^{\\otimes (m - n)}$, where $A$ is the left-upper block of $U$. \n",
    "\n",
    "Observe that $\\mathcal{X} |\\psi_{\\text{good}\\rangle^{\\otimes m} = \\sin(\\frac{\\pi}{2k}) |\\psi\\rangle}$. Therefore, we aim to find a quantum circuit with unitary $U_\\Phi$ such that $W U_\\Phi V = \\frac{1}{\\sin(\\frac{\\pi}{2k})} \\mathcal{X}$, implying $W U_\\Phi V |\\psi_{\\text{good}}\\rangle^{\\otimes m} = |0\\rangle$. By $\\mathcal{X} = W U V$, the absolute values of singular values of $\\mathcal{X}$ are $\\sin(\\frac{\\pi}{2k})$. Moreover, choose $P(x) = (-1)^k T_k (x)$, where $T_k$ is the Chebyshev polynomial of the first kind of order $k$. Then\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",
    "implies that the QSVT of $\\mathcal{X}$ in terms of polynomial $P$ is a block encoding of $B := \\frac{1}{\\sin(\\frac{\\pi}{2k})} A$, as required."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set $k = 3$ so that $\\sin(\\frac{\\pi}{2k}) = \\frac{1}{2}$, and $U$ is a randomly chosen unitary. We want to show that the left-upper block of $U$ is amplified by 2 after QSVT."
   ]
  },
  {
   "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": [
      "the accuracy of quantum singular value transformation is 3.029211939065135e-07\n"
     ]
    }
   ],
   "source": [
    "B = U_Phi[0:2 ** num_block_qubits, 0:2 ** num_block_qubits].numpy()\n",
    "print(f\"the accuracy of quantum singular value transformation is {np.linalg.norm(B - 2 * A)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As shown above, we successfully amplify $\\mathcal{X}$ by 2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## References\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
}