[GQSP] Negative powers and Hamiltonian simulation

qualtran/bloqs/generalized_qsp.ipynb

@@ -15,10 +15,34 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from qualtran.bloqs.qubitization_walk_operator_test import get_walk_operator_for_1d_Ising_model\n",
+    "import numpy as np\n",
+    "import scipy\n",
     "from qualtran.drawing import show_bloq\n",
     "\n",
-    "from qualtran.bloqs.generalized_qsp import GeneralizedQSP"
+    "import qualtran\n",
+    "from qualtran.bloqs.qubitization_walk_operator_test import get_walk_operator_for_1d_Ising_model\n",
+    "from qualtran.bloqs.qubitization_walk_operator import QubitizationWalkOperator\n",
+    "from qualtran.cirq_interop import CirqGateAsBloq\n",
+    "\n",
+    "from qualtran.bloqs.generalized_qsp import GeneralizedQSP, HamiltonianSimulationByGQSP\n",
+    "from qualtran.bloqs.generalized_qsp_test import RandomGate"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6963c30f339d42de",
+   "metadata": {},
+   "source": [
+    "`GeneralizedQSP` implements the Quantum Eigenvalue Transform on a unitary $U$ using QSP. Given a complex GQSP polynomial $P$ (and its complement $Q$), it implements the unitary:\n",
+    "$$U' = \\begin{bmatrix} P(U) & \\cdot \\\\ Q(U) & \\cdot \\end{bmatrix}$$\n",
+    "\n",
+    "Here, the polynomials $P, Q$ must satisfy the following constraint:\n",
+    "\n",
+    "$$\\left| P(e^{i\\theta}) \\right|^2 + \\left| Q(e^{i\\theta}) \\right|^2 = 1 ~~\\text{for every}~ \\theta \\in [0, 2\\pi]$$\n",
+    "\n",
+    "A polynomial $P$ is said to be a GQSP polynomial iff it satisfies $\\left| P(e^{i\\theta}) \\right|^2 \\le 1$ for every $\\theta \\in [0, 2\\pi]$. \n",
+    "\n",
+    "Reference: https://doi.org/10.48550/arXiv.2308.01501"
    ]
   },
   {
@@ -32,42 +56,111 @@
     "show_bloq(U.decompose_bloq())"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dc7100fd94d58c6b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pU = GeneralizedQSP(U, (0.5, 0.5), (-0.5, 0.5))\n",
+    "show_bloq(pU.decompose_bloq())"
+   ]
+  },
   {
    "cell_type": "markdown",
-   "id": "6963c30f339d42de",
+   "id": "935a03f7-5843-4b11-abe6-5eb9048c0ab5",
    "metadata": {},
    "source": [
-    "`GeneralizedQSP` implements the Quantum Eigenvalue Transform on a unitary $U$ using QSP. Given a complex QSP polynomial $P$ (and its complement $Q$), it implements the unitary:\n",
-    "$$U' = \\begin{bmatrix} P(U) & \\cdot \\\\ Q(U) & \\cdot \\end{bmatrix}$$\n",
+    "There is also a method that directly computes $Q$ from $P$:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "78cd3857297f092b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pU = GeneralizedQSP.from_qsp_polynomial(U, (0.5, 0, 0.5))\n",
+    "show_bloq(pU.decompose_bloq())"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a58f06ba-9287-435d-92a3-256f747024c2",
+   "metadata": {},
+   "source": [
+    "### Negative degree terms\n",
     "\n",
-    "Here, the polynomials $P, Q$ must satisfy the following constraint:\n",
+    "To apply GQSP for a polynomial $P'(z) = z^{-k} P(z)$, we can just pass the polynomial $P$ along with negative power $k$.\n",
+    "The QSP angle sequence is the same for both, and $P'$ can be achieved by running $(U^\\dagger)^k$ at any point in the circuit."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ee60e95f-979f-420d-91b1-b6955b9c5a3b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pU = GeneralizedQSP.from_qsp_polynomial(U, (0.5, 0, 0.5), negative_power=1)\n",
+    "show_bloq(pU)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c84e5443-483c-4294-8a79-1d2bed70e1f5",
+   "metadata": {},
+   "source": [
+    "## Hamiltonian Simulation by GQSP\n",
+    "\n",
+    "Given the Szegedy Quantum Walk Operator for a Hamiltonian $H$, one can construct the walk operator for $e^{-iHt}$ using GQSP (Corollary 8).\n",
+    "\n",
+    "### Recap:\n",
+    "\n",
+    "For a Hamiltonian $H = \\sum_i \\alpha_i U_i$, given the SELECT and PREPARE oracles\n",
+    "$$ \\text{SELECT} = \\sum_j \\ketbra{j}{j} \\otimes U_j $$\n",
+    "$$ \\text{PREPARE} \\ket{0} = \\sum_j \\frac{\\sqrt{\\alpha_j}}{\\|\\alpha\\|_1} \\ket{j} $$\n",
+    "we can implement the [QubitizationWalkOperator](qubitization_walk_operator.ipynb) that encodes the spectrum of $H$ in the eigenphases of the walk operator $W$.\n",
     "\n",
-    "$$\\left\\mid P(e^{i\\theta}) \\right\\mid^2 + \\left\\mid Q(e^{i\\theta}) \\right\\mid^2 = 1 ~~\\text{for every}~ \\theta \\in [0, 2\\pi]$$\n",
+    "### Approximating $\\cos$\n",
+    "We can use the Jacobi-Anger expansion to obtain low-degree polynomial approximations for the $\\cos$ function:\n",
     "\n",
+    "$$e^{it\\cos\\theta} = \\sum_{n = -\\infty}^{\\infty} i^n J_n(t) (e^{i\\theta})^n$$\n",
     "\n",
-    "Reference: https://arxiv.org/abs/2308.01501"
+    "where $J_n$ is the $n$-th [Bessel function of the first kind](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_first_kind), which is provided by `scipy.special.jv`.\n",
+    "\n",
+    "We can cutoff at $d = O(t + \\log(1/\\epsilon) / \\log\\log(1/\\epsilon))$ to get an $\\epsilon$-approximation (Theorem 7):\n",
+    "\n",
+    "$$P[t](z) = \\sum_{n = -d}^d i^n J_n(t) z^n$$\n",
+    "\n",
+    "### Obtaining $e^{iHt}$\n",
+    "\n",
+    "As the eigenphases of the walk operator are $e^{-i\\arccos(E_k / \\|\\alpha\\|_1)}$, we can use the GQSP polynomial with $P = P[\\|\\alpha\\|_1 t]$ (and complementary $Q = 0$) to obtain $P(U) = e^{-iHt}$.\n",
+    "The obtained GQSP operator $W'$ can then be used with the PREPARE oracle to simulate the hamiltonian:\n",
+    "$$(\\langle0| \\otimes \\text{PREPARE}^\\dagger \\otimes I) W' (|0\\rangle \\otimes \\text{PREPARE} \\otimes I) |0\\rangle|\\psi\\rangle = |0\\rangle e^{-iHt}|\\psi\\rangle$$"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "dc7100fd94d58c6b",
+   "id": "c4ccd9c8-92ce-4b53-92bc-722f263a32ad",
    "metadata": {},
    "outputs": [],
    "source": [
-    "pU = GeneralizedQSP(U, (0.5, 0.5))\n",
-    "show_bloq(pU.decompose_bloq())"
+    "W_e_iHt = HamiltonianSimulationByGQSP(U, t=5, alpha=1, precision=1e-3)\n",
+    "show_bloq(W_e_iHt.decompose_bloq())"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "78cd3857297f092b",
+   "id": "c6627cb6-446d-4f76-b727-8ed2f67bc382",
    "metadata": {},
    "outputs": [],
    "source": [
-    "pU = GeneralizedQSP(U, (0.5, 0, 0.5))\n",
-    "show_bloq(pU.decompose_bloq())"
+    "W_e_iHt.call_graph()"
    ]
   }
  ],
@@ -87,7 +180,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.9"
+   "version": "3.11.7"
   }
  },
  "nbformat": 4,