diff --git a/applications/Block_encoding-ND_Laplacian/1D_Lap_BE.png b/applications/Block_encoding-ND_Laplacian/1D_Lap_BE.png new file mode 100644 index 000000000..d9caef196 Binary files /dev/null and b/applications/Block_encoding-ND_Laplacian/1D_Lap_BE.png differ diff --git a/applications/Block_encoding-ND_Laplacian/1D_Periodic_Laplacian_BE.qmod b/applications/Block_encoding-ND_Laplacian/1D_Periodic_Laplacian_BE.qmod new file mode 100644 index 000000000..828f03574 --- /dev/null +++ b/applications/Block_encoding-ND_Laplacian/1D_Periodic_Laplacian_BE.qmod @@ -0,0 +1,25 @@ +qfunc subtractor_mod(a: qnum) { + a += -1; +} + +qfunc adder_mod(a: qnum) { + a += 1; +} + +qfunc main(output j: qnum, output l: qnum) { + allocate(5, j); + allocate(2, l); + apply_to_all(H, l); + apply_to_all(Z, l); + l_arr: qbit[]; + l -> l_arr; + control (l_arr[0] == 0) { + subtractor_mod(j); + } + control (l_arr[1] == 1) { + adder_mod(j); + } + l_arr -> l; + apply_to_all(H, l); + apply_to_all(Z, l); +} diff --git a/applications/Block_encoding-ND_Laplacian/1D_Periodic_Laplacian_BE.synthesis_options.json b/applications/Block_encoding-ND_Laplacian/1D_Periodic_Laplacian_BE.synthesis_options.json new file mode 100644 index 000000000..9ed1d3f02 --- /dev/null +++ b/applications/Block_encoding-ND_Laplacian/1D_Periodic_Laplacian_BE.synthesis_options.json @@ -0,0 +1,45 @@ +{ + "constraints": { + "max_gate_count": {}, + "max_width": 18, + "optimization_parameter": "depth" + }, + "preferences": { + "custom_hardware_settings": { + "basis_gates": [ + "tdg", + "z", + "sx", + "u1", + "cx", + "h", + "p", + "x", + "ry", + "s", + "cy", + "t", + "r", + "cz", + "u", + "rx", + "y", + "id", + "sdg", + "u2", + "rz", + "sxdg" + ], + "is_symmetric_connectivity": true + }, + "debug_mode": true, + "machine_precision": 8, + "optimization_level": 1, + "output_format": ["qasm"], + "pretty_qasm": true, + "random_seed": 185454979, + "synthesize_all_separately": false, + "timeout_seconds": 300, + "transpilation_option": "auto optimize" + } +} diff --git a/applications/Block_encoding-ND_Laplacian/1D_structure.png b/applications/Block_encoding-ND_Laplacian/1D_structure.png new file mode 100644 index 000000000..95e429916 Binary files /dev/null and b/applications/Block_encoding-ND_Laplacian/1D_structure.png differ diff --git a/applications/Block_encoding-ND_Laplacian/1d_discretization.png b/applications/Block_encoding-ND_Laplacian/1d_discretization.png new file mode 100644 index 000000000..8535d4958 Binary files /dev/null and b/applications/Block_encoding-ND_Laplacian/1d_discretization.png differ diff --git a/applications/Block_encoding-ND_Laplacian/2D_Lap_structure.png b/applications/Block_encoding-ND_Laplacian/2D_Lap_structure.png new file mode 100644 index 000000000..e77b626d0 Binary files /dev/null and b/applications/Block_encoding-ND_Laplacian/2D_Lap_structure.png differ diff --git a/applications/Block_encoding-ND_Laplacian/3D_Lap_structure.png b/applications/Block_encoding-ND_Laplacian/3D_Lap_structure.png new file mode 100644 index 000000000..044e47f62 Binary files /dev/null and b/applications/Block_encoding-ND_Laplacian/3D_Lap_structure.png differ diff --git a/applications/Block_encoding-ND_Laplacian/Ddim_Lap_BE.png b/applications/Block_encoding-ND_Laplacian/Ddim_Lap_BE.png new file mode 100644 index 000000000..7bddfcdcc Binary files /dev/null and b/applications/Block_encoding-ND_Laplacian/Ddim_Lap_BE.png differ diff --git a/applications/Block_encoding-ND_Laplacian/ND_Laplacian_BE.ipynb b/applications/Block_encoding-ND_Laplacian/ND_Laplacian_BE.ipynb new file mode 100644 index 000000000..f2735d11a --- /dev/null +++ b/applications/Block_encoding-ND_Laplacian/ND_Laplacian_BE.ipynb @@ -0,0 +1,791 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# A block encoding circuit for a D-dimensional Laplacian with periodic boundary conditions\n", + "\n", + "___\n", + "\n", + "\n", + "This work is an implementation of the paper on **Efficient and Explicit Block Encoding of Finite Difference Discretizations of the Laplacian** (https://arxiv.org/abs/2509.02429)." + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "### **Introduction**\n", + "\n", + "Block encoding is a well known technique in quantum computing used to embed non-unitary matrices on a quantum computer that only allows for unitary evolution. \n", + "\n", + "***Definition:*** \n", + " Let $ a, n, m \\in \\mathbb{N} $ such that $ m = a + n $. A $ m $ -qubit unitary $ U $ is said to be an $(\\alpha, a)$ -block-encoding of an $ n $ -qubit operator $ A $ if \n", + "\n", + "$$\n", + "\\tilde{A} = \\left( \\langle 0 |^{\\otimes a} \\otimes I_n \\right) U \\left( |0 \\rangle^{\\otimes a} \\otimes I_n \\right)\n", + "\\tag{1}\n", + "$$ \n", + "\n", + "where, $A = \\alpha \\tilde{A}$ . The parameters $(\\alpha, a)$ represent the *subnormalization factor* (which adjusts for encoding matrices of any norm), and the *number of ancilla qubits* used in the block-encoding scheme respectively. \n", + "\n", + "Efficiently block encoding arbitrary matrices is a very difficult problem and this task is not trivial even for well structured and sparse matrices. The paper by Sturm et al. 2025 (https://arxiv.org/abs/2509.02429), provides for efficient quantum circuits for block encoding $N-$ dimensional Laplacians, with periodic boundary conditions along all dimensions. This is useful in many applications, especially ones involving problems of linear algebra. Moreover, given an efficient block encoding of a matrix $\\tilde{A}$, its possible to efficiently construct a block encoding of certain polynomials of $\\tilde{A}$ through quantum singular value transformations (QSVT). \n", + "\n", + "### **Notebook contents**\n", + "- ##### Block encoding circuits for D- dimensional Laplacian matrix with periodic boundary conditions along all dimnesions.\n", + "___" + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "## Discretizing the $D$-Dimensional Laplacian with Periodic Boundary Conditions\n", + "\n", + "We consider the Laplacian operator on the $D$-dimensional unit hypercube \n", + "$$\n", + "\\Omega_D = [0,1]^D ,\n", + "$$\n", + "where a point is written as \n", + "$$\n", + "\\mathbf{x} = (x^{(0)}, x^{(1)}, \\ldots, x^{(D-1)}).\n", + "$$\n", + "\n", + "\n", + "The $D$-dimensional Laplacian is the sum of second derivatives along each coordinate axis:\n", + "\\begin{equation}\n", + "L_D = \\sum_{d=0}^{D-1} \\frac{\\partial^2}{\\partial (x^{(d)})^2}.\n", + "\\tag{1}\n", + "\\end{equation}\n", + "\n", + "We impose **periodic boundary conditions**, meaning the function repeats itself across opposite sides of the domain:\n", + "$$\n", + "v(x^{(0)},\\ldots,0,\\ldots,x^{(D-1)}) \n", + "= \n", + "v(x^{(0)},\\ldots,1,\\ldots,x^{(D-1)}).\n", + "$$\n", + "\n", + "This ensures the domain “wraps around” in **every dimension**.\n", + "\n", + "\n", + "To approximate the operator on a computer, we replace each interval with a **uniform grid** (*we assume equidistant grid points as per the paper*). \n", + "\n", + "For one dimension, choose\n", + "$$\n", + "\\Omega_{1,h} = \\{ jh \\mid j=0,1,\\dots,N-1 \\},\n", + "\\qquad h = \\frac{1}{N}.\n", + "$$\n", + "Because of periodicity, the point at $ x=1 $ is identical to $ x=0 $, so we only keep $ N $ samples.\n", + "\n", + "

\n", + " \n", + "

\n", + "\n", + "To build a $D$-dimensional grid, take the Cartesian product:\n", + "$$\n", + "\\Omega_{D,h} = \\Omega_{1,h} \\times \\cdots \\times \\Omega_{1,h}.\n", + "$$\n", + "\n", + "Each grid point is indexed by a $D$-tuple:\n", + "$$\n", + "(j^{(0)}, j^{(1)},\\ldots, j^{(D-1)}),\n", + "\\qquad j^{(d)} \\in \\{0,\\dots,N-1\\}.\n", + "$$\n", + "\n", + "The total number of points is \n", + "$$\n", + "N_D = N^D.\n", + "$$\n", + "\n", + "\n", + "To represent the function values as a vector, we map each integer index $j^{(i)}$ to a $2^n$ dimensional vector $\\ket{j^{(i)}}$ in the Hilbert sapce, where $N=2^n$. Thus, \n", + "$$\n", + "(j^{(0)},j^{(1)},\\ldots,j^{(D-1)})\n", + "\\quad \\longmapsto \\quad\n", + "j = \\ket{j^{(D-1)}}\\otimes \\cdots \\otimes\\ket{j^{(1)}}\\otimes\\ket{j^{(0)}},\n", + "$$\n", + "\n", + "This ordering turns the entire grid $\\Omega_{D,h}$ into a vector of $\\mathbf{\\Omega_{D,h}}$ size $N^D= 2^{nD}$ where,\n", + "\n", + "$$\n", + "\\mathbf{\\Omega_{D,h}}\n", + "= \\left( j^{(0)} h,\\; j^{(1)} h,\\; \\ldots,\\; j^{(D-1)} h \\right)\n", + "\\,\\lvert j^{(D-1)} \\rangle \\cdots \\lvert j^{(1)} \\rangle \\lvert j^{(0)} \\rangle .\n", + "$$\n", + "\n", + "___\n", + "\n", + "#### Building the 1D Finite-Difference Laplacian Matrix (Periodic Boundary Conditions)\n", + "\n", + "To construct the discrete Laplacian on a 1D periodic grid, we approximate the second\n", + "derivative using the standard centered finite-difference stencil.\n", + "\n", + "\n", + "We discretize the interval $[0,1]$ using \n", + "$$\n", + "x_j = jh, \\quad j = 0,1,\\dots, N-1, \\qquad h = \\frac{1}{N}.\n", + "$$\n", + "\n", + "Because the boundary is periodic, the point at $x=1$ is *identified* with $x=0$.\n", + "So the grid has exactly $N$ degrees of freedom.\n", + "\n", + "\n", + "\n", + "For a smooth function $u(x)$, the second derivative at grid point $x_j$ is approximated by\n", + "\n", + "$$\n", + "\\frac{d^2u}{dx^2}(x_j)\n", + "\\approx \n", + "\\frac{u_{j-1} - 2u_j + u_{j+1}}{h^2}.\n", + "$$\n", + "\n", + "This is the classic *three-point centered stencil*.\n", + "\n", + "Without periodicity, $u_{j-1}$ and $u_{j+1}$ would fail at the boundaries. \n", + "With periodicity, the grid “wraps around,” so:\n", + "\n", + "$$\n", + "u_{-1} \\equiv u_{N-1}, \\qquad \n", + "u_{N} \\equiv u_0.\n", + "$$\n", + "\n", + "This single idea gives the Laplacian its circulant structure.\n", + "\n", + "\n", + "For each point $j$, the finite-difference rule says:\n", + "\n", + "- coefficient of $u_{j-1}$ is $+1/h^2$\n", + "- coefficient of $u_{j}$ is $-2/h^2$\n", + "- coefficient of $u_{j+1}$ is $+1/h^2$\n", + "\n", + "Putting these coefficients into matrix form gives an $N \\times N$ matrix where:\n", + "\n", + "- the main diagonal contains $-2/h^2$\n", + "- the first upper and first lower diagonal contain $1/h^2$\n", + "- periodicity adds *two extra wrap-around entries*\n", + "\n", + "Specifically:\n", + "\n", + "- entry $(0, N-1)$ gets $+1/h^2$ (from $u_{-1}=u_{N-1}$)\n", + "- entry $(N-1, 0)$ gets $+1/h^2$ (from $u_N=u_0$)\n", + "\n", + "\n", + "Putting this together, the discrete $1D$ Laplacian is\n", + "\n", + "\\begin{equation}\n", + "L_{1,h} = \\frac{1}{h^2}\n", + "\\begin{pmatrix}\n", + "-2 & 1 & 0 & \\cdots & 0 & 1\\\\\n", + "1 & -2 & 1 & \\cdots & 0 & 0\\\\\n", + "0 & 1 & -2 & \\cdots & 0 & 0\\\\\n", + "\\vdots & \\vdots & \\vdots & \\ddots & \\vdots & \\vdots \\\\\n", + "1 & 0 & 0 & \\cdots & 1 & -2\n", + "\\end{pmatrix}.\n", + "\\tag{2}\n", + "\\end{equation}\n", + "\n", + "This matrix is the building block for higher-dimensional Laplacians via Kronecker sums.\n", + "\n", + "___\n", + "\n", + "\n", + "#### Building the D- dimensional Finite-Difference Laplacian Matrix (Periodic Boundary Conditions)\n", + "\n", + "\n", + "The discrete Laplacian on the full $D$-dimensional grid is a **Kronecker sum** of 1D Laplacians:\n", + "\n", + "\\begin{equation}\n", + "L_{D,h} \n", + "= \n", + "\\sum_{d=0}^{D-1}\n", + "I_N^{\\otimes (D-1-d)} \\otimes L_{1,h} \\otimes I_N^{\\otimes d}.\n", + "\\tag{3}\n", + "\\end{equation}\n", + "\n", + "Each term contributes the second derivative along one dimension, while identities preserve all other coordinates. \n", + "This operator is sparse: every grid point couples only to its $2D$ nearest periodic neighbors.\n", + "\n", + "\n", + "The maximum eigenvalue magnitude of the discrete Laplacian is known analytically:\n", + "$$\n", + "\\lambda_{D, \\max} = \\frac{4D}{h^2}.\n", + "$$\n", + "\n", + "To make the matrix have operator norm 1 (required for several quantum algorithms), define the **scaled Laplacian**:\n", + "$$\n", + "\\widetilde{L}_{D,h}\n", + "= \n", + "\\frac{1}{\\lambda_{D,\\max}} L_{D,h}\n", + "= \n", + "\\frac{1}{4D/h^2} L_{D,h}.\n", + "$$\n", + "\n", + "Plugging this into the Kronecker-sum representation yields\n", + "\n", + "\\begin{equation}\n", + "\\widetilde{L}_{D,h}\n", + "=\n", + "\\frac{1}{D}\n", + "\\sum_{d=0}^{D-1}\n", + "I_N^{\\otimes (D-1-d)} \\otimes \\widetilde{L}_{1,h} \\otimes I_N^{\\otimes d},\n", + "\\tag{4}\n", + "\\end{equation}\n", + "where\n", + "$$\n", + "\\widetilde{L}_{1,h} = \\frac{h^2}{4} L_{1,h}\n", + "$$\n", + "is the 1D Laplacian scaled to have norm 1.\n", + "\n", + "We now provide the block encoding of the scaled Laplacian matrix. \n", + "___" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "#### Block encoding the 1D scaled Laplacian matrix\n", + "\n", + "With periodic boundary conditions, the 1D Laplacian matrix has a circulant structure as shown (for $N= 16$ grid points): \n", + "\n", + "\n", + "

\n", + " \n", + "

\n", + "\n", + "This structure can be block encoded via a simple quantum circuit shown below: \n", + "\n", + "

\n", + " \n", + "

\n", + "\n", + "Here $S_{-}$ and $S_{+}$ are unitary operations that take $\\ket{j}$ to $\\ket{(j+1) \\mod{N}}$ and $\\ket{(j-1) \\mod{N}}$ respectively. \n", + "\n", + "One can easily see that the scaled 1D Laplacian can be written in terms of linear combination of unitaries as:\n", + "\n", + "\\begin{equation}\n", + "\\widetilde L_{1} \\;:=\\; \\frac{1}{4}\\,(S_- - 2I + S_+).\n", + "\\tag{5}\n", + "\\end{equation}\n", + "___\n", + "\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, + "source": [ + "#### Necessary modules and global functions" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "# import all necessary modules\n", + "import builtins\n", + "import math\n", + "import random\n", + "import time\n", + "from typing import List\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "import classiq\n", + "from classiq import *\n", + "from classiq.qmod.symbolic import *\n", + "\n", + "# np.set_printoptions(precision=2, suppress=True, linewidth=120, threshold=10000) # print options for neatness\n", + "\n", + "\n", + "# post-processing function to get reduced state vector after projection on ancilla=0\n", + "def get_projected_state_vector(\n", + " execution_result,\n", + " measured_var: str,\n", + " projections: dict,\n", + ") -> np.ndarray:\n", + " \"\"\"\n", + " This function returns a reduced statevector from execution results.\n", + " measured var: the name of the reduced variable\n", + " projections: on which values of the other variables to project, e.g., {\"anc_M\": 1}\n", + " Note: For this function to work properly all variables, except auxiliary qubits, must be declared as\n", + " output of the model.\n", + " \"\"\"\n", + " projected_size = len(execution_result.output_qubits_map[measured_var])\n", + " proj_statevector = np.zeros(2**projected_size).astype(complex)\n", + " for sample in execution_result.parsed_state_vector:\n", + " if all(\n", + " int(sample.state[key]) == projections[key] for key in projections.keys()\n", + " ):\n", + " value = int(sample.state[measured_var])\n", + " proj_statevector[value] += sample.amplitude\n", + " global_phase = np.angle(proj_statevector[0])\n", + " return np.real(proj_statevector / np.exp(1j * global_phase))\n", + "\n", + "\n", + "def fidelity_with_phase_alignment(psi, phi, align=True):\n", + " \"\"\"\n", + " Compute fidelity between two pure-state vectors (||^2).\n", + " If `align` is True, first remove the global phase of `phi` so the\n", + " inner product with `psi` is real-positive (this also fixes an overall sign).\n", + " Inputs may be unnormalized; they are normalized inside.\n", + " Returns fidelity (float). If return_alignment True returns tuple\n", + " (fidelity, phi_aligned, applied_phase).\n", + " \"\"\"\n", + "\n", + " psi = np.asarray(psi, dtype=complex).ravel()\n", + " phi = np.asarray(phi, dtype=complex).ravel()\n", + " if psi.size != phi.size:\n", + " raise ValueError(\"State vectors must have same dimension\")\n", + "\n", + " npsi = np.linalg.norm(psi)\n", + " nphi = np.linalg.norm(phi)\n", + " if npsi == 0 or nphi == 0:\n", + " raise ValueError(\"Zero vector provided\")\n", + "\n", + " psi = psi / npsi\n", + " phi = phi / nphi\n", + "\n", + " overlap = np.vdot(psi, phi) # \n", + " if align:\n", + " phase = np.angle(overlap)\n", + " phi_aligned = phi * np.exp(-1j * phase)\n", + " overlap_aligned = np.vdot(psi, phi_aligned)\n", + " fidelity = float(np.abs(overlap_aligned) ** 2)\n", + " else:\n", + " phase = 0.0\n", + " phi_aligned = phi\n", + " fidelity = float(np.abs(overlap) ** 2)\n", + "\n", + " return fidelity" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "# classical function to generate scaled 1D Laplacian matrix with periodic boundary conditions\n", + "\n", + "\n", + "def shift_plus(N):\n", + " S = np.zeros((N, N), dtype=complex)\n", + " for j in range(N):\n", + " S[(j + 1) % N, j] = 1.0\n", + " return S\n", + "\n", + "\n", + "def shift_minus(N):\n", + " S = np.zeros((N, N), dtype=complex)\n", + " for j in range(N):\n", + " S[(j - 1) % N, j] = 1.0\n", + " return S\n", + "\n", + "\n", + "def laplacian_1d(N):\n", + " \"\"\"Scaled 1D Laplacian \\tilde L on N grid points (periodic).\"\"\"\n", + " return 0.25 * (shift_minus(N) - 2 * np.eye(N) + shift_plus(N))" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "7", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Circuit Width: 7\n", + "Circuit Depth: 113\n", + "Gate counts: {'u': 111, 'cx': 100}\n", + "The reduced state vector for j when l=0, k=0 is:\n", + "[ 5.00000000e-01 -2.50000000e-01 2.03934950e-16 2.60524259e-16\n", + " -3.00593240e-16 -3.80160209e-17 8.15406091e-17 -4.52150298e-17\n", + " -6.30058944e-17 1.21989910e-16 -4.46287283e-18 5.15617198e-17\n", + " -8.55801898e-17 -1.61909579e-17 2.01956495e-17 -1.05456418e-17\n", + " -1.37383090e-16 9.49128933e-17 -1.17137621e-18 -3.88260210e-17\n", + " 2.36256948e-17 3.94809326e-17 -3.59332343e-17 -9.83622130e-18\n", + " 7.68836822e-17 -2.74072948e-18 -4.98416893e-17 -3.17034970e-16\n", + " 2.84043112e-16 9.53120283e-17 -1.58750885e-16 -2.50000000e-01]\n", + "Expected theoretical state vector:\n", + "[-0.5 +0.j 0.25+0.j 0. +0.j 0. +0.j 0. +0.j 0. +0.j 0. +0.j\n", + " 0. +0.j 0. +0.j 0. +0.j 0. +0.j 0. +0.j 0. +0.j 0. +0.j\n", + " 0. +0.j 0. +0.j 0. +0.j 0. +0.j 0. +0.j 0. +0.j 0. +0.j\n", + " 0. +0.j 0. +0.j 0. +0.j 0. +0.j 0. +0.j 0. +0.j 0. +0.j\n", + " 0. +0.j 0. +0.j 0. +0.j 0.25+0.j]\n", + "********\n" + ] + } + ], + "source": [ + "N = 32\n", + "n = int(math.log2(N))\n", + "lap_scaled = laplacian_1d(N)\n", + "\n", + "\n", + "@qfunc\n", + "def adder_mod(a: QNum):\n", + " a += 1\n", + "\n", + "\n", + "@qfunc\n", + "def subtractor_mod(a: QNum):\n", + " a += -1\n", + "\n", + "\n", + "@qfunc\n", + "def main(j: Output[QNum], l: Output[QNum]):\n", + " allocate(n, j)\n", + " allocate(2, l)\n", + "\n", + " # we chose the input state |0> for testing\n", + "\n", + " apply_to_all(H, l)\n", + " apply_to_all(Z, l)\n", + "\n", + " l_arr = QArray()\n", + " bind(l, l_arr)\n", + "\n", + " control(l_arr[0] == 0, lambda: subtractor_mod(j)), control(\n", + " l_arr[1] == 1, lambda: adder_mod(j)\n", + " )\n", + "\n", + " bind(l_arr, l)\n", + "\n", + " apply_to_all(H, l)\n", + " apply_to_all(Z, l)\n", + "\n", + "\n", + "qmod = create_model(main)\n", + "backend_preferences = ClassiqBackendPreferences(backend_name=\"simulator_statevector\")\n", + "qmod = set_execution_preferences(\n", + " qmod,\n", + " execution_preferences=ExecutionPreferences(\n", + " num_shots=1, backend_preferences=backend_preferences\n", + " ),\n", + ")\n", + "constraints = Constraints(\n", + " max_width=18, optimization_parameter=OptimizationParameter.DEPTH\n", + ")\n", + "qmod = set_constraints(qmod, constraints)\n", + "qprog = synthesize(qmod)\n", + "# show(qprog)\n", + "\n", + "# write_qmod(qmod, \"1D_Periodic_Laplacian_BE\")\n", + "\n", + "print(\"Circuit Width:\", qprog.data.width)\n", + "print(\"Circuit Depth:\", qprog.transpiled_circuit.depth)\n", + "print(\"Gate counts:\", qprog.transpiled_circuit.count_ops)\n", + "\n", + "job = execute(qprog)\n", + "results = job.result()[0].value\n", + "\n", + "# Post processing\n", + "reduced_state = get_projected_state_vector(results, \"j\", {\"l\": 0})\n", + "print(\"The reduced state vector for j when l=0, k=0 is:\")\n", + "print(reduced_state)\n", + "\n", + "# #theoretical reduced state vector\n", + "\n", + "b = [0 for _ in range(N)]\n", + "b[0] = 1\n", + "theoretical_reduced_state = np.matmul(lap_scaled, b)\n", + "print(f\"Expected theoretical state vector:\")\n", + "print(theoretical_reduced_state)\n", + "print(\"********\")" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "8", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Fidelity between theoretical and obtained reduced state vector: 1.0\n" + ] + } + ], + "source": [ + "fidelity = fidelity_with_phase_alignment(reduced_state, theoretical_reduced_state)\n", + "print(\n", + " f\"Fidelity between theoretical and obtained reduced state vector: {np.round(fidelity, 4)}\"\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "#### **Note:** The state vector results should match upto a overall global phase. Hence the above theoretical and simulated results match exactly. \n", + "\n", + "___" + ] + }, + { + "cell_type": "markdown", + "id": "10", + "metadata": {}, + "source": [ + "#### Block encoding D-dimensional scaled Laplacian matrix\n", + "\n", + "\n", + "The scaled D-dimensional laplacian from equation $4$ is as follows:\n", + "\n", + "$$\n", + "\\widetilde{L}_{D,h}\n", + "=\n", + "\\frac{1}{D}\n", + "\\sum_{d=0}^{D-1}\n", + "I_N^{\\otimes (D-1-d)} \\otimes \\widetilde{L}_{1,h} \\otimes I_N^{\\otimes d},\n", + "$$\n", + "\n", + "where we can use the expression from equation $5$ for $\\widetilde{L}_{1,h}$.\n", + "\n", + "The structure of the matrices for $D=2$ and $D=3$ is shown below:\n", + "\n", + "
\n", + " \"D-dim\n", + " \"D-dim\n", + "
\n", + "\n", + "\n", + "\n", + "Due to the above decomposition, we can easily extend the idea of 1D Laplacian block encoding to a D-dimensional laplacian by introducing an additional control qubit register that controls on the dimension number, as follows: \n", + "\n", + "

\n", + " \n", + "

\n", + "\n", + "here $\\hat{d}= \\lceil{\\log D}\\rceil$\n", + "\n", + "___" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "# classical function to generate scaled 1D Laplacian matrix with periodic boundary conditions\n", + "\n", + "\n", + "def laplacian_multiD(Ns):\n", + " \"\"\"\n", + " Build D-dimensional Laplacian with periodic BCs.\n", + " Ns = [N0, N1, ..., N_{D-1}], where Ni = 2**n_i.\n", + " \"\"\"\n", + " D = len(Ns)\n", + " N_total = np.prod(Ns)\n", + " L_total = np.zeros((N_total, N_total), dtype=complex)\n", + "\n", + " for i, Ni in enumerate(Ns):\n", + " Li = laplacian_1d(Ni)\n", + " # Build tensor with reversed order (so last entry in Ns = qubit 0)\n", + " kron_terms = []\n", + " for j, Nj in enumerate(Ns):\n", + " if j == i:\n", + " kron_terms.append(Li)\n", + " else:\n", + " kron_terms.append(np.eye(Nj, dtype=complex))\n", + " # reverse list for little-endian convention\n", + " kron_terms = kron_terms[::-1]\n", + " Li_full = kron_terms[0]\n", + " for term in kron_terms[1:]:\n", + " Li_full = np.kron(Li_full, term)\n", + " L_total += Li_full\n", + "\n", + " d = math.ceil(math.log2(builtins.max(1, D)))\n", + "\n", + " return L_total * (1 / (2**d)) # scale by 1/D" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "12", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Circuit Width: 6\n", + "Circuit Depth: 78\n", + "Gate counts: {'u': 56, 'cx': 48}\n", + "The reduced state vector for j when l=0, k=0 is:\n", + "[ 5.00000000e-01 -2.50000000e-01 -1.25000000e-01 -3.61302947e-17\n", + " -1.71584428e-16 3.26277430e-17 -1.25000000e-01 5.27811857e-17]\n", + "Expected theoretical state vector:\n", + "[-0.5 +0.j 0.25 +0.j 0.125+0.j 0. +0.j 0. +0.j 0. +0.j\n", + " 0.125+0.j 0. +0.j]\n", + "********\n" + ] + } + ], + "source": [ + "# Classical Inputs: N_i, D\n", + "D = 2 # number of dimensions\n", + "d = math.ceil(math.log2(builtins.max(1, D))) # d hat of the paper\n", + "N = [2, 4] # grid size in each dimension is 4, 16\n", + "n = [int(math.log2(i)) for i in N] # number of qubits in each dimension\n", + "total_system_size = np.sum(n)\n", + "lap_scaled = laplacian_multiD(N)\n", + "\n", + "\n", + "@qfunc\n", + "def adder_mod(a: QNum):\n", + " a += 1\n", + "\n", + "\n", + "@qfunc\n", + "def subtractor_mod(a: QNum):\n", + " a += -1\n", + "\n", + "\n", + "def control_adder_mod(j: List[QNum], l: QArray, control_bit: QNum):\n", + " for i in range(D):\n", + " control(\n", + " control_bit == i,\n", + " lambda: [\n", + " control(l[0] == 0, lambda: subtractor_mod(j[i])),\n", + " control(l[1] == 1, lambda: adder_mod(j[i])),\n", + " ],\n", + " )\n", + "\n", + "\n", + "@qfunc\n", + "def main(j: Output[QNum], l: Output[QNum], k: Output[QNum]):\n", + " allocate(total_system_size, j)\n", + " allocate(2, l)\n", + " allocate(d, k)\n", + "\n", + " # we chose the input state |0> for testing\n", + " apply_to_all(H, l)\n", + " apply_to_all(Z, l)\n", + " apply_to_all(H, k)\n", + "\n", + " j_regs = [QNum(f\"j{i}\", n[i], False, 0) for i in range(len(n))]\n", + " l_arr = QArray()\n", + "\n", + " bind(j, j_regs)\n", + " bind(l, l_arr)\n", + "\n", + " control_adder_mod(j_regs, l_arr, k)\n", + "\n", + " bind(l_arr, l)\n", + " bind(j_regs, j)\n", + "\n", + " apply_to_all(H, l)\n", + " apply_to_all(H, k)\n", + "\n", + "\n", + "qmod = create_model(main)\n", + "qmod = set_execution_preferences(\n", + " qmod,\n", + " execution_preferences=ExecutionPreferences(\n", + " num_shots=1, backend_preferences=backend_preferences\n", + " ),\n", + ")\n", + "qmod = set_constraints(qmod, constraints)\n", + "qprog = synthesize(qmod)\n", + "# show(qprog)\n", + "\n", + "write_qmod(qmod, \"ND_Periodic_Laplacian_BE\")\n", + "\n", + "print(\"Circuit Width:\", qprog.data.width)\n", + "print(\"Circuit Depth:\", qprog.transpiled_circuit.depth)\n", + "print(\"Gate counts:\", qprog.transpiled_circuit.count_ops)\n", + "\n", + "job = execute(qprog)\n", + "results = job.result()[0].value\n", + "\n", + "# Post processing\n", + "reduced_state = get_projected_state_vector(results, \"j\", {\"l\": 0, \"k\": 0})\n", + "print(\"The reduced state vector for j when l=0, k=0 is:\")\n", + "print(reduced_state)\n", + "\n", + "# #theoretical reduced state vector\n", + "b = [0 for i in range(N[0] * N[1])]\n", + "b[0] = 1\n", + "theoretical_reduced_state = np.matmul(lap_scaled, b)\n", + "print(f\"Expected theoretical state vector:\")\n", + "print(theoretical_reduced_state)\n", + "print(\"********\")" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "13", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Fidelity between theoretical and obtained reduced state vector: 1.0\n" + ] + } + ], + "source": [ + "fidelity = fidelity_with_phase_alignment(reduced_state, theoretical_reduced_state)\n", + "print(\n", + " f\"Fidelity between theoretical and obtained reduced state vector: {np.round(fidelity, 4)}\"\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "14", + "metadata": {}, + "source": [ + "### Reference\n", + "\n", + "- Sturm, A., & Schillo, N. (2025). Efficient and Explicit Block Encoding of Finite Difference Discretizations of the Laplacian. arXiv preprint https://arxiv.org/abs/2509.02429" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "qiskit-fresh", + "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.10.18" + } + }, + "nbformat": 4, + "nbformat_minor": 9 +} diff --git a/applications/Block_encoding-ND_Laplacian/ND_Periodic_Laplacian_BE.qmod b/applications/Block_encoding-ND_Laplacian/ND_Periodic_Laplacian_BE.qmod new file mode 100644 index 000000000..c6354d224 --- /dev/null +++ b/applications/Block_encoding-ND_Laplacian/ND_Periodic_Laplacian_BE.qmod @@ -0,0 +1,41 @@ +qfunc subtractor_mod(a: qnum) { + a += -1; +} + +qfunc adder_mod(a: qnum) { + a += 1; +} + +qfunc main(output j: qnum, output l: qnum, output k: qnum) { + allocate(3, j); + allocate(2, l); + allocate(1, k); + apply_to_all(H, l); + apply_to_all(Z, l); + apply_to_all(H, k); + j0: qnum<1, False, 0>; + j1: qnum<2, False, 0>; + l_arr: qbit[]; + j -> {j0, j1}; + l -> l_arr; + control (k == 0) { + control (l_arr[0] == 0) { + subtractor_mod(j0); + } + control (l_arr[1] == 1) { + adder_mod(j0); + } + } + control (k == 1) { + control (l_arr[0] == 0) { + subtractor_mod(j1); + } + control (l_arr[1] == 1) { + adder_mod(j1); + } + } + l_arr -> l; + {j0, j1} -> j; + apply_to_all(H, l); + apply_to_all(H, k); +} diff --git a/applications/Block_encoding-ND_Laplacian/ND_Periodic_Laplacian_BE.synthesis_options.json b/applications/Block_encoding-ND_Laplacian/ND_Periodic_Laplacian_BE.synthesis_options.json new file mode 100644 index 000000000..b5abe982c --- /dev/null +++ b/applications/Block_encoding-ND_Laplacian/ND_Periodic_Laplacian_BE.synthesis_options.json @@ -0,0 +1,45 @@ +{ + "constraints": { + "max_gate_count": {}, + "max_width": 18, + "optimization_parameter": "depth" + }, + "preferences": { + "custom_hardware_settings": { + "basis_gates": [ + "id", + "u1", + "u", + "sx", + "sxdg", + "sdg", + "y", + "u2", + "r", + "z", + "h", + "s", + "cy", + "tdg", + "cz", + "cx", + "rx", + "x", + "t", + "ry", + "rz", + "p" + ], + "is_symmetric_connectivity": true + }, + "debug_mode": true, + "machine_precision": 8, + "optimization_level": 1, + "output_format": ["qasm"], + "pretty_qasm": true, + "random_seed": 3612357240, + "synthesize_all_separately": false, + "timeout_seconds": 300, + "transpilation_option": "auto optimize" + } +} diff --git a/applications/block_encoding/2D_Lap.png b/applications/block_encoding/2D_Lap.png new file mode 100644 index 000000000..7dbdf6139 Binary files /dev/null and b/applications/block_encoding/2D_Lap.png differ diff --git a/applications/block_encoding/2D_Laplacian_circuit.png b/applications/block_encoding/2D_Laplacian_circuit.png new file mode 100644 index 000000000..3d49bf376 Binary files /dev/null and b/applications/block_encoding/2D_Laplacian_circuit.png differ diff --git a/applications/block_encoding/General_BE.png b/applications/block_encoding/General_BE.png new file mode 100644 index 000000000..6399413cf Binary files /dev/null and b/applications/block_encoding/General_BE.png differ diff --git a/applications/block_encoding/checkerboard.png b/applications/block_encoding/checkerboard.png new file mode 100644 index 000000000..5c0894704 Binary files /dev/null and b/applications/block_encoding/checkerboard.png differ diff --git a/applications/block_encoding/checkerboard_BE_N8.qmod b/applications/block_encoding/checkerboard_BE_N8.qmod new file mode 100644 index 000000000..0152906ff --- /dev/null +++ b/applications/block_encoding/checkerboard_BE_N8.qmod @@ -0,0 +1,45 @@ +qfunc array_swap(j: qnum, k: qnum) { + z1: qbit[]; + z2: qbit[]; + j -> z1; + k -> z2; + repeat (i: z1.len) { + SWAP(z1[i], z2[i]); + } + z1 -> j; + z2 -> k; +} + +qfunc oracle_d(s: qnum, data: qnum, TM: real[]) { + repeat (i: TM::len) { + control (s == i) { + RX(2 * acos(TM[i]), data); + } + } +} + +qfunc checkerboard_BE(data: qnum, s: qnum, j: qnum) { + s0: qnum<1, False, 0>; + s1: qnum<2, False, 0>; + j0: qnum<2, False, 0>; + j1: qnum<1, False, 0>; + s -> {s0, s1}; + j -> {j0, j1}; + array_swap(s1, j0); + CX(s0, j1); + oracle_d(s0, data, [(-0.69), (-0.34)]); + {s0, s1} -> s; + {j0, j1} -> j; +} + +qfunc main(output data: qnum, output s: qnum, output j: qnum) { + allocate(3, s); + allocate(3, j); + allocate(1, data); + apply_to_all(H, j); + within { + apply_to_all(H, s); + } apply { + checkerboard_BE(data, s, j); + } +} diff --git a/applications/block_encoding/checkerboard_BE_N8.synthesis_options.json b/applications/block_encoding/checkerboard_BE_N8.synthesis_options.json new file mode 100644 index 000000000..0538cc4e6 --- /dev/null +++ b/applications/block_encoding/checkerboard_BE_N8.synthesis_options.json @@ -0,0 +1,45 @@ +{ + "constraints": { + "max_gate_count": {}, + "max_width": 10, + "optimization_parameter": "depth" + }, + "preferences": { + "custom_hardware_settings": { + "basis_gates": [ + "z", + "rz", + "u", + "u2", + "sx", + "cx", + "cz", + "cy", + "y", + "tdg", + "r", + "rx", + "sdg", + "t", + "x", + "s", + "u1", + "ry", + "p", + "id", + "sxdg", + "h" + ], + "is_symmetric_connectivity": true + }, + "debug_mode": true, + "machine_precision": 8, + "optimization_level": 3, + "output_format": ["qasm"], + "pretty_qasm": true, + "random_seed": 4262529357, + "synthesize_all_separately": false, + "timeout_seconds": 1500, + "transpilation_option": "auto optimize" + } +} diff --git a/applications/block_encoding/checkerboard_circuit.png b/applications/block_encoding/checkerboard_circuit.png new file mode 100644 index 000000000..be9e14961 Binary files /dev/null and b/applications/block_encoding/checkerboard_circuit.png differ diff --git a/applications/block_encoding/laplacian_BE_Nx4_Ny4.qmod b/applications/block_encoding/laplacian_BE_Nx4_Ny4.qmod new file mode 100644 index 000000000..ef4fe785f --- /dev/null +++ b/applications/block_encoding/laplacian_BE_Nx4_Ny4.qmod @@ -0,0 +1,101 @@ +qfunc eq_sup_lap(s: qnum) { + inplace_prepare_state([0.2, 0, 0.2, 0.2, 0.2, 0.2, 0, 0], 0, s); +} + +qfunc oracle_org(s: qbit[], j: qnum, dlt: qbit) { + j1: qnum<2.0, False, 0>; + j2: qnum<2.0, False, 0>; + s0: qnum<1, False, 0>; + s1: qnum<1, False, 0>; + s2: qnum<1, False, 0>; + s -> {s0, s1, s2}; + j -> {j1, j2}; + control (s0 == 1) { + control (s1 == 0) { + control (s2 == 0) { + X(dlt); + } + } + } + control (s1 == 1) { + control (j1 == 0) { + X(dlt); + } + } + control (s2 == 1) { + control (j2 == 0) { + X(dlt); + } + } + {j1, j2} -> j; + {s0, s1, s2} -> s; +} + +qfunc oracle_d(s: qnum, data: qnum, TM: real[]) { + repeat (i: TM::len) { + control (s == i) { + RX(2 * acos(TM[i]), data); + } + } +} + +qfunc laplacian_BE(s: qnum, j: qnum, data: qnum, dlt: qnum) { + j1: qnum<2, False, 0>; + j2: qnum<2, False, 0>; + s0: qnum<1, False, 0>; + s1: qnum<1, False, 0>; + s2: qnum<1, False, 0>; + s -> {s0, s1, s2}; + j -> {j1, j2}; + control (s0 == 0) { + control (s2 == 1) { + j2 += 1; + } + } + {j1, j2} -> j; + control (s1 == 1) { + control (s0 == 0) { + j += 1; + } + } + {s0, s1, s2} -> s; + oracle_org(s, j, dlt); + s -> {s0, s1, s2}; + d: qnum<2, False, 0>; + {s1, s2} -> d; + Z(data); + oracle_d(d, data, [(-1.0), 0.25, 0.25]); + d -> {s1, s2}; + control (s1 == 1) { + X(s0); + } + control (s2 == 1) { + X(s0); + } + control (s1 == 1) { + control (s0 == 0) { + j += -1; + } + } + j -> {j1, j2}; + control (s2 == 1) { + control (s0 == 0) { + j2 += -1; + } + } + {j1, j2} -> j; + {s0, s1, s2} -> s; +} + +qfunc main(output s: qnum, output j: qnum, output data: qnum, output dlt: qbit) { + allocate(3, s); + allocate(1, dlt); + allocate(1, data); + allocate(4, j); + apply_to_all(H, j); + within { + eq_sup_lap(s); + } apply { + laplacian_BE(s, j, data, dlt); + } +} diff --git a/applications/block_encoding/laplacian_BE_Nx4_Ny4.synthesis_options.json b/applications/block_encoding/laplacian_BE_Nx4_Ny4.synthesis_options.json new file mode 100644 index 000000000..627cafd7f --- /dev/null +++ b/applications/block_encoding/laplacian_BE_Nx4_Ny4.synthesis_options.json @@ -0,0 +1,45 @@ +{ + "constraints": { + "max_gate_count": {}, + "max_width": 10, + "optimization_parameter": "depth" + }, + "preferences": { + "custom_hardware_settings": { + "basis_gates": [ + "z", + "rz", + "u", + "u2", + "sx", + "cx", + "cz", + "cy", + "y", + "tdg", + "r", + "rx", + "sdg", + "t", + "x", + "s", + "u1", + "ry", + "p", + "id", + "sxdg", + "h" + ], + "is_symmetric_connectivity": true + }, + "debug_mode": true, + "machine_precision": 8, + "optimization_level": 1, + "output_format": ["qasm"], + "pretty_qasm": true, + "random_seed": 1506744138, + "synthesize_all_separately": false, + "timeout_seconds": 300, + "transpilation_option": "auto optimize" + } +} diff --git a/applications/block_encoding/select_structures_BE.ipynb b/applications/block_encoding/select_structures_BE.ipynb new file mode 100644 index 000000000..d86c9bcb8 --- /dev/null +++ b/applications/block_encoding/select_structures_BE.ipynb @@ -0,0 +1,1374 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Block encoding of select structured matrices\n", + "___\n", + "\n", + "This work is an implementation of the paper on **Block-encoding structured matrices for data input in quantum computing** (https://quantum-journal.org/papers/q-2024-01-11-1226/)." + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "### **Introduction**\n", + "\n", + "Block encoding is a well known technique in quantum computing used to embed non-unitary matrices on a quantum computer that only allows for unitary evolution. \n", + "\n", + "***Definition:*** \n", + " Let $ a, n, m \\in \\mathbb{N} $ such that $ m = a + n $. A $ m $ -qubit unitary $ U $ is said to be an $(\\alpha, a)$ -block-encoding of an $ n $ -qubit operator $ A $ if \n", + "\n", + "$$\n", + "\\tilde{A} = \\left( \\langle 0 |^{\\otimes a} \\otimes I_n \\right) U \\left( |0 \\rangle^{\\otimes a} \\otimes I_n \\right)\n", + "\\tag{1}\n", + "$$ \n", + "\n", + "where, $A = \\alpha \\tilde{A}$ . The parameters $(\\alpha, a)$ represent the *subnormalization factor* (which adjusts for encoding matrices of any norm), and the *number of ancilla qubits* used in the block-encoding scheme respectively. \n", + "\n", + "Efficiently block encoding arbitrary matrices is a very difficult problem and this task is not trivial even for well structured and sparse matrices. The paper by Sunderhauf et al. 2024 (https://quantum-journal.org/papers/q-2024-01-11-1226/), provides for efficient quantum circuits for block encoding arithmetically structured matrices. This is useful in many applications, especially ones involving problems of linear algebra. Moreover, given an efficient block encoding of a matrix $\\tilde{A}$, its possible to efficiently construct a block encoding of certain polynomials of $\\tilde{A}$ through quantum singular value transformations (QSVT). \n", + "\n", + "### **Notebook contents**\n", + "- ##### Block encoding circuits for Checkerboard matrix, Toeplitz matrix, Tridiagonal symmetric matrix and 2D Laplacian matrix.\n", + "___" + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "### **Notation for Sparse Matrix Structures**\n", + "\n", + "We consider an $N \\times N$ sparse matrix $A$ whose nonzero entries are drawn from a fixed set of $D$ distinct data values:\n", + "$$\n", + "A_d \\quad (0 \\le d < D)\n", + "$$\n", + "Multiple positions in the matrix may contain the same data value.\n", + "\n", + "We denote:\n", + "- $M$: maximum number of times any of the $D$ values appears in the matrix (multiplicity of a data value)\n", + "- $S_c$: maximum number of nonzero entries in any column \n", + "- $S_r$: maximum number of nonzero entries in any row \n", + "- $S = \\max(S_c, S_r)$: overall sparsity level (smaller $S$ means a sparser matrix)\n", + "\n", + "We assume the following to hold for this work:\n", + "- all $D$ data items have the same multiplicity $M$,\n", + "- all rows and columns have equal sparsity: $S_c = S_r = S$,\n", + "- all quantities are powers of two.\n", + "\n", + "Then, the total number of nonzero entries satisfies\n", + "$$\n", + "MD = NS_c = NS_r = \\#\\text{nonzero}. \n", + "$$\n", + "\n", + "**Note:** If the above equality does not hold, we pad $S$ and/or $D$ till the equality is satisfied. \n", + "\n", + "Each nonzero entry can then be equivalently uniquely labeled from three perspectives:\n", + "1. **By data value and repetition index:** $(d, m)$ \n", + "2. **By column index and sparsity index:** $(j, s_c)$ \n", + "3. **By row index and sparsity index:** $(i, s_r)$ \n", + "\n", + "These perspectives describe the same sparsity pattern using different coordinate systems.\n", + "\n", + "\n", + "### **Block encoding ingredients**\n", + "\n", + "We first label each matrix entry by a data label, say $(d,m). Then, we would require the following ingredients to block encode the matrix given above.\n", + "\n", + "- Oracle $ O_c$: $$ O_c |d, m\\rangle \\rightarrow\\left|j, s_c\\right\\rangle $$ where $s_c$ is the column sparsity index, $0 \\le s_c \\le S_c$.\n", + "- Oracle $ O_r$: $$ O_r |d, m\\rangle \\rightarrow\\left|i, s_r\\right\\rangle $$ where $s_r$ is the row sparsity index, $0 \\le s_r \\le S_r$.\n", + "\n", + "The above two oracles are solely dependent on the structure of the matrix in question and can be constructed through appropriate quantum arithmetic. \n", + "\n", + "- Apart from these, a dataloading oracle encoding the values of the $D$ data items is required:\n", + "$$\n", + "O_{\\text {data }}=\\sum_{d=0}^{D-1} R_X\\left(2 \\arccos A_d /\\|A\\|_{\\max }\\right) \\otimes|d\\rangle\\langle d|\n", + "$$\n", + "\n", + "The factor $\\|A\\|_{\\max }=\\max _d\\left|A_d\\right|$ is needed to ensure all values are in-range of the arccosine. Provided the first qubit is initialised and postselected as $|0\\rangle$, the effect of $O_{\\text {data }}$ is loading the correct values:\n", + "\n", + "$$\n", + "O_{\\text {data }}|0\\rangle \\otimes|d\\rangle=\\left(A_d /\\|A\\|_{\\max }|0\\rangle-i \\sqrt{1-\\left(\\left|A_d\\right| /\\|A\\|_{\\max }\\right)^2}|1\\rangle\\right) \\otimes|d\\rangle .\n", + "$$\n", + "\n", + "### **General Block encoding circuit**\n", + "\n", + "

\n", + " \n", + "

\n", + "\n", + "\n", + "#### **Important Note**\n", + "- In the below illustrations we scale the matrix entires of $A$ by dividing the unscaled $A$ by $||A||_{max}$ so that entries of scaled $A_s$ are between $-1$ and $+1$. We block encode the scaled matrix $A_s$ (where $||A_{s}||_{max}= 1$). \n", + "\n", + "___" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "#### Necessary modules and global functions" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "# import all necessary modules\n", + "import builtins\n", + "import math\n", + "import random\n", + "import time\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "np.set_printoptions(\n", + " precision=2, suppress=True, linewidth=120, threshold=10000\n", + ") # print options for neatness\n", + "\n", + "import classiq\n", + "from classiq import *\n", + "from classiq.qmod.symbolic import *\n", + "\n", + "\n", + "# function to extract unique values from structured matrices\n", + "def _unique_values_core(matrix, mode, Nx=None):\n", + " \"\"\"\n", + " helper to extract unique values (as needed for O_data oracle) from structured matrices.\n", + "\n", + " Parameters\n", + " ----------\n", + " matrix : array-like or sparse-like\n", + " Input matrix (Laplacian, tridiagonal, Toeplitz).\n", + " mode : {'laplacian', 'tridiagonal', 'toeplitz'}\n", + " Extraction mode.\n", + " Nx : int, optional\n", + " Grid size in x-direction (only needed for 'laplacian').\n", + "\n", + " Returns\n", + " -------\n", + " list\n", + " List of unique values\n", + " \"\"\"\n", + "\n", + " mat = matrix\n", + "\n", + " if mode == \"laplacian\":\n", + " if Nx is None:\n", + " raise ValueError(\"Nx must be provided for mode='laplacian'.\")\n", + " # 2D Laplacian: pick the main diag, 1-step, and Nx-step diagonals\n", + " vals = [\n", + " mat.diagonal(0)[0], # A0\n", + " mat.diagonal(1)[0], # A1\n", + " mat.diagonal(Nx)[0], # A2\n", + " ]\n", + " return vals\n", + "\n", + " elif mode == \"tridiagonal\":\n", + " size = mat.shape[0]\n", + " vals = []\n", + "\n", + " for i in range(size):\n", + " for j in range(size):\n", + " v = mat[i, j]\n", + " if v not in vals:\n", + " vals.append(v)\n", + "\n", + " return vals\n", + "\n", + " elif mode == \"toeplitz\":\n", + " N = mat.shape[0]\n", + " vals = []\n", + " for offset in range(-N + 1, N):\n", + " diagonal = np.diag(mat, k=offset)\n", + " if diagonal.size > 0:\n", + " v = diagonal[0]\n", + " if v != 0 and v not in vals:\n", + " vals.append(v)\n", + " return vals[::-1]\n", + "\n", + " else:\n", + " raise ValueError(f\"Unknown mode '{mode}'.\")\n", + "\n", + "\n", + "# post-processing function to get reduced state vector after projection on ancilla=0\n", + "def get_projected_state_vector(\n", + " execution_result,\n", + " measured_var: str,\n", + " projections: dict,\n", + ") -> np.ndarray:\n", + " \"\"\"\n", + " This function returns a reduced statevector from execution results.\n", + " measured var: the name of the reduced variable\n", + " projections: on which values of the other variables to project, e.g., {\"anc_M\": 1}\n", + " \"\"\"\n", + " projected_size = len(execution_result.output_qubits_map[measured_var])\n", + " proj_statevector = np.zeros(2**projected_size).astype(complex)\n", + " for sample in execution_result.parsed_state_vector:\n", + " if all(\n", + " int(sample.state[key]) == projections[key] for key in projections.keys()\n", + " ):\n", + " value = int(sample.state[measured_var])\n", + " proj_statevector[value] += sample.amplitude\n", + " global_phase = np.angle(proj_statevector[0])\n", + " return np.real(proj_statevector / np.exp(1j * global_phase))\n", + "\n", + "\n", + "# general oracle to load data values from a classical array TM into data register using angle encoding\n", + "@qfunc\n", + "def oracle_d(s: QNum, data: QNum, TM: CArray[CReal]):\n", + " repeat(\n", + " TM.len,\n", + " lambda i: control(\n", + " s == i, lambda: (a_d := TM[i], theta := 2 * acos(a_d), RX(theta, data))\n", + " ),\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "### 1. Checkerboard Matrix" + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "\n", + "\n", + "A checkerboard matrix $A\\in\\mathbb{R}^{N\\times N}$ is any matrix whose entries alternate between two scalar values depending on the parity of the index sum.\n", + "\n", + "Given two values $a_0,a_1\\in[-1,1]$, we define\n", + "$$\n", + "\n", + "A_{ij} = \n", + "\\begin{cases}\n", + "a_0, & (i+j)\\ \\text{even},\\\\[4pt]\n", + "a_1, & (i+j)\\ \\text{odd},\n", + "\\end{cases} \\qquad i,j=0,\\dots,N-1.\n", + "\n", + "$$\n", + "Thus the matrix contains only two distinct values and a simple periodic pattern\n", + "that is useful for testing structure-aware block-encoding. An example heatmap showing the structure is given below:\n", + "\n", + "

\n", + " \n", + "

\n", + "\n", + "For convenience, the index $m\\in\\{0,\\dots,N^2/2-1\\}$ is decomposed into three parts,\n", + "$$\n", + "m \\;=\\; N\\,m^{\\mathrm{hi}} \\;+\\; \\frac{N}{2}\\,m^{\\mathrm{mid}} \\;+\\; m^{\\mathrm{lo}},\n", + "$$\n", + "corresponding respectively to the high, mid, and low $\\log(N/2)$-bit segments.\n", + "Given a data index $d\\in\\{0,1\\}$ and multiplicity index $m$, the row and column\n", + "indices $(i,j)$ of the checkerboard matrix are obtained as\n", + "\\begin{align}\n", + "i(d,m) &= \\left\\lfloor \\frac{m}{N/2} \\right\\rfloor\n", + " \\;=\\; 2\\,m^{\\mathrm{hi}} + m^{\\mathrm{mid}}, \\\\[6pt]\n", + "j(d,m) &= 2\\,\\big(m \\bmod (N/2)\\big) \\;+\\; \\big(d + i(d,m)\\bmod 2\\big)\n", + " \\;=\\; 2\\,m^{\\mathrm{lo}} + \\big(d + m^{\\mathrm{mid}} \\bmod 2\\big).\n", + "\\end{align}\n", + "These relations give a compact arithmetic mapping from $(d,m)$ labels to the corresponding matrix entry $(i,j)$ in the $N\\times N$ checkerboard matrix. One can deduce that $O_r$ would simply be the identity operator and the block encoding circuit is provided as below. \n", + "\n", + "\n", + "

\n", + " \n", + "

\n" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "Here the subnormalization factor is $N$ and $n+1$ ancilla qubits are used, where $n= log_2N$. Greater the subnormalization and ancilla qubits greater is the block encoding cost, so for efficient block encodings one tries to minimize these resources. " + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "8", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Checkerboard matrix:\n", + "\n", + "[[-0.08 -0.84 -0.08 -0.84 -0.08 -0.84 -0.08 -0.84]\n", + " [-0.84 -0.08 -0.84 -0.08 -0.84 -0.08 -0.84 -0.08]\n", + " [-0.08 -0.84 -0.08 -0.84 -0.08 -0.84 -0.08 -0.84]\n", + " [-0.84 -0.08 -0.84 -0.08 -0.84 -0.08 -0.84 -0.08]\n", + " [-0.08 -0.84 -0.08 -0.84 -0.08 -0.84 -0.08 -0.84]\n", + " [-0.84 -0.08 -0.84 -0.08 -0.84 -0.08 -0.84 -0.08]\n", + " [-0.08 -0.84 -0.08 -0.84 -0.08 -0.84 -0.08 -0.84]\n", + " [-0.84 -0.08 -0.84 -0.08 -0.84 -0.08 -0.84 -0.08]]\n", + "Values (val_even, val_odd): [-0.08, -0.84]\n" + ] + } + ], + "source": [ + "def checkerboard_matrix(N, val_even=None, val_odd=None, rng=None, decimals=2):\n", + " \"\"\"\n", + " Generate an N×N checkerboard matrix with alternating values.\n", + "\n", + " Parameters\n", + " ----------\n", + " N : int\n", + " Matrix dimension.\n", + "\n", + " val_even : float, optional\n", + " Value at positions where (i + j) is even. Random in [-1,1] if None.\n", + "\n", + " val_odd : float, optional\n", + " Value at positions where (i + j) is odd. Random in [-1,1] if None.\n", + "\n", + " rng : np.random.Generator, optional\n", + " Random generator for reproducibility.\n", + "\n", + " decimals : int, optional\n", + " Decimal precision for the values.\n", + "\n", + " Returns\n", + " -------\n", + " np.ndarray, list\n", + " The checkerboard matrix and the list [val_even, val_odd].\n", + " \"\"\"\n", + "\n", + " if rng is None:\n", + " rng = np.random.default_rng()\n", + "\n", + " # Randomly choose the two values if not provided and round them\n", + " if val_even is None:\n", + " val_even = float(rng.uniform(-1, 1))\n", + " else:\n", + " val_even = float(val_even)\n", + " if val_odd is None:\n", + " val_odd = float(rng.uniform(-1, 1))\n", + " else:\n", + " val_odd = float(val_odd)\n", + "\n", + " val_even = round(val_even, decimals)\n", + " val_odd = round(val_odd, decimals)\n", + "\n", + " # checkerboard pattern\n", + " i = np.arange(N).reshape(-1, 1)\n", + " j = np.arange(N).reshape(1, -1)\n", + " parity = (i + j) % 2 # 0 for even, 1 for odd positions\n", + "\n", + " M = np.where(parity == 0, val_even, val_odd).astype(float)\n", + "\n", + " # Round matrix entries to the requested number of decimals\n", + " M = np.round(M, decimals)\n", + "\n", + " return M, [val_even, val_odd]\n", + "\n", + "\n", + "# Test case\n", + "if __name__ == \"__main__\":\n", + " N = 8\n", + " M, values = checkerboard_matrix(N)\n", + " print(\"Checkerboard matrix:\\n\")\n", + " print(M)\n", + " print(\"Values (val_even, val_odd):\", values)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "9", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The circuit width is: 7\n", + "The circuit depth is: 11\n", + "The reduced state vector for j when s=0, data=0 is:\n", + "[1.3 1.3 1.3 1.3 1.3 1.3 1.3 1.3]\n", + "Theoretical reduced state vector for j when l=0 and data=0 is:\n", + "[-1.3 -1.3 -1.3 -1.3 -1.3 -1.3 -1.3 -1.3]\n", + "The reduced state vector for j when s=0, data=0 is:\n", + "[1.3 1.3 1.3 1.3 1.3 1.3 1.3 1.3]\n", + "Theoretical reduced state vector for j when l=0 and data=0 is:\n", + "[-1.3 -1.3 -1.3 -1.3 -1.3 -1.3 -1.3 -1.3]\n" + ] + } + ], + "source": [ + "@qfunc\n", + "def array_swap(j: QNum, k: QNum):\n", + " z_arr1 = QArray(\"z1\")\n", + " z_arr2 = QArray(\"z2\")\n", + " bind(j, z_arr1)\n", + " bind(k, z_arr2)\n", + " repeat(z_arr1.len, lambda i: SWAP(z_arr1[i], z_arr2[i]))\n", + " bind(z_arr1, j)\n", + " bind(z_arr2, k)\n", + "\n", + "\n", + "@qfunc\n", + "def checkerboard_BE(data: QNum, s: QNum, j: QNum):\n", + " n = int(math.log2(N / 2))\n", + "\n", + " s0 = QNum(\"s0\", 1, False, 0)\n", + " s1 = QNum(\"s1\", n, False, 0)\n", + "\n", + " j0 = QNum(\"j0\", n, False, 0)\n", + " j1 = QNum(\"j1\", 1, False, 0)\n", + "\n", + " bind(s, [s0, s1])\n", + " bind(j, [j0, j1])\n", + "\n", + " array_swap(s1, j0)\n", + " CX(s0, j1)\n", + "\n", + " oracle_d(\n", + " s0, data, values\n", + " ) # refers to O_data which is common across all block encodings in this notebook\n", + "\n", + " bind([s0, s1], s)\n", + " bind([j0, j1], j)\n", + "\n", + "\n", + "@qfunc\n", + "def main(data: Output[QNum], s: Output[QNum], j: Output[QNum]):\n", + " n = int(math.log2(N))\n", + "\n", + " allocate(n, s)\n", + " allocate(n, j)\n", + " allocate(1, data)\n", + "\n", + " apply_to_all(H, j)\n", + "\n", + " within_apply(lambda: apply_to_all(H, s), lambda: checkerboard_BE(data, s, j))\n", + "\n", + "\n", + "preferences = Preferences(timeout_seconds=1500, optimization_level=3)\n", + "constraints = Constraints(\n", + " optimization_parameter=OptimizationParameter.DEPTH,\n", + " max_width=10, # change max_width as needed, minimum width is 2n+1\n", + ")\n", + "qmod = create_model(main, preferences=preferences)\n", + "qmod = set_constraints(qmod, constraints)\n", + "\n", + "# write_qmod(\n", + "# qmod,\n", + "# 'checkerboard_BE_N{}'.format(N)\n", + "# )\n", + "\n", + "# Synthesize the quantum program\n", + "backend_preferences = ClassiqBackendPreferences(backend_name=\"simulator_statevector\")\n", + "qmod = set_execution_preferences(\n", + " qmod,\n", + " execution_preferences=ExecutionPreferences(\n", + " num_shots=1, backend_preferences=backend_preferences\n", + " ),\n", + ")\n", + "qprog = synthesize(qmod)\n", + "# show(qprog)\n", + "circuit_width = qprog.data.width\n", + "circuit_depth = qprog.transpiled_circuit.depth\n", + "print(\"The circuit width is:\", circuit_width)\n", + "print(\"The circuit depth is:\", circuit_depth)\n", + "\n", + "job = execute(qprog)\n", + "results = job.result()[0].value\n", + "\n", + "\n", + "reduced_state = get_projected_state_vector(results, \"j\", {\"s\": 0, \"data\": 0})\n", + "reduced_state = reduced_state * N # rescale by subnormalization factor N\n", + "print(\"The reduced state vector for j when s=0, data=0 is:\")\n", + "print(reduced_state)\n", + "\n", + "\n", + "# theoretical reduced state vector\n", + "theoretical_reduced_state = np.matmul(M, [1 / np.sqrt(N) for i in range(N)])\n", + "print(\"Theoretical reduced state vector for j when l=0 and data=0 is:\")\n", + "print(theoretical_reduced_state)" + ] + }, + { + "cell_type": "markdown", + "id": "10", + "metadata": {}, + "source": [ + "#### **Note:** The state vector results should match upto a overall global phase. Hence the above theoretical and simulated results match exactly. " + ] + }, + { + "cell_type": "markdown", + "id": "11", + "metadata": {}, + "source": [ + "____" + ] + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "### 2. Toeplitz Matrix" + ] + }, + { + "cell_type": "markdown", + "id": "13", + "metadata": {}, + "source": [ + "Consider a $N$ x $N$ Toeplitz matrix with $D$ diagonals (or $D$ values), offset from the main diagonal by $k$ :\n", + "\n", + "$$\n", + "\\left(\\begin{array}{cccccccc}\n", + "A_k & A_{k-1} & \\cdots & A_0 & & & & \\\\\n", + "A_{k+1} & A_k & A_{k-1} & \\cdots & A_0 & & & \\\\\n", + "\\vdots & A_{k+1} & A_k & A_{k-1} & \\cdots & A_0 & & \\\\\n", + "A_{D-1} & \\vdots & A_{k+1} & A_k & A_{k-1} & \\cdots & A_0 & \\\\\n", + "& A_{D-1} & \\vdots & A_{k+1} & A_k & A_{k-1} & \\cdots & A_0 \\\\\n", + "& & A_{D-1} & \\vdots & A_{k+1} & A_k & A_{k-1} & \\cdots \\\\\n", + "& & & A_{D-1} & \\vdots & A_{k+1} & A_k & A_{k-1} \\\\\n", + "& & & & A_{D-1} & \\vdots & A_{k+1} & A_k\n", + "\\end{array}\\right)\n", + "$$\n", + "\n", + "For the distinct values, we choose $d$ as equal to the subscript of the matrix elements above. For $m$, we simply choose the column index. Arithmetically, the mapping to row ($i$) and column indices ($j$) is then:\n", + "\n", + "$$\n", + "i(d, m)=d-k+m, \\quad j(d, m)=m\n", + "$$\n", + "\n", + "and out-of-range $(d, m)$ pairs are those where an overflow or underflow in the calculation of $i$ occurs, that is when\n", + "\n", + "$$\n", + "d-k+m<0 \\text { or } d-k+m \\geq N \\text {. }\n", + "$$\n", + "\n", + "Note that the range of the data labels are as follows:\n", + "\n", + "$i,j,m: \\{0,1,.., N-1\\} , d: \\{0,1,..., D-1\\}$\n", + "\n", + "\n", + "In the case of Topelitz matrix, elements in any row and any column has a unique data label $d$ hence we do not require a separate arithmetic for $s_c$ and $s_r$, they can be set equal to $d$ as their range also match. \n", + "\n", + "#### Block encoding circuit:\n", + "\n", + "

\n", + " \n", + "

\n", + "\n", + "Here the subnormalization factor is $D$ and $2+log_2D$ ancilla qubits are used." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "14", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Toeplitz Matrix:\n", + " [[0.45 0.54 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]\n", + " [0.86 0.45 0.54 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]\n", + " [0.93 0.86 0.45 0.54 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]\n", + " [0. 0.93 0.86 0.45 0.54 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]\n", + " [0. 0. 0.93 0.86 0.45 0.54 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]\n", + " [0. 0. 0. 0.93 0.86 0.45 0.54 0. 0. 0. 0. 0. 0. 0. 0. 0. ]\n", + " [0. 0. 0. 0. 0.93 0.86 0.45 0.54 0. 0. 0. 0. 0. 0. 0. 0. ]\n", + " [0. 0. 0. 0. 0. 0.93 0.86 0.45 0.54 0. 0. 0. 0. 0. 0. 0. ]\n", + " [0. 0. 0. 0. 0. 0. 0.93 0.86 0.45 0.54 0. 0. 0. 0. 0. 0. ]\n", + " [0. 0. 0. 0. 0. 0. 0. 0.93 0.86 0.45 0.54 0. 0. 0. 0. 0. ]\n", + " [0. 0. 0. 0. 0. 0. 0. 0. 0.93 0.86 0.45 0.54 0. 0. 0. 0. ]\n", + " [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.93 0.86 0.45 0.54 0. 0. 0. ]\n", + " [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.93 0.86 0.45 0.54 0. 0. ]\n", + " [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.93 0.86 0.45 0.54 0. ]\n", + " [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.93 0.86 0.45 0.54]\n", + " [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.93 0.86 0.45]]\n", + "Unique values: [np.float64(0.5430706158006551), np.float64(0.4451945190809161), np.float64(0.8588474400968182), np.float64(0.9258545815408299)]\n" + ] + } + ], + "source": [ + "# Toeplitz matrix\n", + "def generate_toeplitz_matrix(D, N, k):\n", + " \"\"\"\n", + " Generate a Toeplitz matrix with D diagonals and size NxN.\n", + "\n", + " Parameters:\n", + " D (int): Number of diagonals, must be a power of 2.\n", + " N (int): Size of the matrix, must be a power of 2 and N >= D.\n", + " k (int): Offset from the main diagonal.\n", + "\n", + " Returns:\n", + " numpy.ndarray: NxN Toeplitz matrix.\n", + " \"\"\"\n", + " # Check if D and N are powers of 2\n", + " if not (D & (D - 1) == 0 and D > 0):\n", + " raise ValueError(\"D must be a power of 2.\")\n", + " if not (N & (N - 1) == 0 and N > 0):\n", + " raise ValueError(\"N must be a power of 2.\")\n", + " if N < D:\n", + " raise ValueError(\"N must be greater than or equal to D.\")\n", + " if k < 0 or k >= D:\n", + " raise ValueError(\"k must be between 0 and D-1.\")\n", + "\n", + " # Generate random diagonal values between -1 and 1\n", + " diagonals = np.random.uniform(-1, 1, D)\n", + "\n", + " # Generate the Toeplitz matrix\n", + " toeplitz_matrix = np.zeros((N, N))\n", + " for offset in range(D):\n", + " diag_index = k - offset\n", + " if diag_index >= 0:\n", + " np.fill_diagonal(toeplitz_matrix[:, diag_index:], diagonals[offset])\n", + " if diag_index < 0:\n", + " np.fill_diagonal(toeplitz_matrix[-diag_index:, :], diagonals[offset])\n", + "\n", + " return toeplitz_matrix\n", + "\n", + "\n", + "def get_unique_toeplitz(toeplitz_matrix):\n", + " \"\"\"\n", + " Get the list of unique values used in the Toeplitz matrix,\n", + " ordered from 0th to (D-1)th diagonal\n", + " \"\"\"\n", + " return _unique_values_core(toeplitz_matrix, mode=\"toeplitz\")\n", + "\n", + "\n", + "# Test case\n", + "D = 4 # Number of diagonals (must be a power of 2, for now)\n", + "N = 16 # Size of the matrix (must be a power of 2 and >= D)\n", + "k = 1 # Offset from the main diagonal, k must be between 0 and D-1\n", + "\n", + "toeplitz_matrix = generate_toeplitz_matrix(D, N, k)\n", + "print(\"Toeplitz Matrix:\\n\", toeplitz_matrix)\n", + "\n", + "unique_values_toeplitz = get_unique_toeplitz(\n", + " toeplitz_matrix\n", + ") # unique values in the matrix (each has a unique label- d)\n", + "print(\"Unique values:\", unique_values_toeplitz)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "15", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/opt/anaconda3/envs/qiskit-fresh/lib/python3.10/site-packages/classiq/synthesis.py:60: ClassiqDeprecationWarning: The computation sequence of the local variable 'tmp' includes a non-permutation operation\n", + "\t\tat ipynb cell line 13 in function 'Toep_BE'\n", + " quantum_program.raise_warnings()\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The circuit width is: 18\n", + "The circuit depth is: 225\n", + "The reduced state vector for j when s=0, dlt=1 and data=0 is:\n", + "[0.25 0.46 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.56]\n", + "Theoretical reduced state vector for j when l=0 and data=0 is:\n", + "[0.25 0.46 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.56]\n", + "The reduced state vector for j when s=0, dlt=1 and data=0 is:\n", + "[0.25 0.46 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.56]\n", + "Theoretical reduced state vector for j when l=0 and data=0 is:\n", + "[0.25 0.46 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.56]\n" + ] + } + ], + "source": [ + "@qfunc\n", + "def oracle_r(j: QNum, s: QNum, k: CInt):\n", + " inplace_add(s - k, j)\n", + "\n", + "\n", + "@qfunc\n", + "def Toep_BE(data: QNum, s: QNum, j: QNum, dlt: QBit):\n", + " T_M = unique_values_toeplitz\n", + " oracle_d(s, data, T_M)\n", + "\n", + " tmp = QNum(\"tmp\")\n", + " tmp |= s - k + j\n", + "\n", + " oracle_r(j, s, k)\n", + " control(tmp < 0, lambda: X(dlt))\n", + " control(tmp >= N, lambda: X(dlt))\n", + "\n", + "\n", + "@qfunc\n", + "def main(data: Output[QNum], s: Output[QNum], j: Output[QNum], dlt: Output[QBit]):\n", + " d = int(math.log2(D))\n", + " n = int(math.log2(N))\n", + "\n", + " allocate(d, s)\n", + " allocate(n, j)\n", + " allocate(1, data)\n", + " allocate(1, dlt)\n", + "\n", + " apply_to_all(H, j)\n", + "\n", + " within_apply(lambda: apply_to_all(H, s), lambda: Toep_BE(data, s, j, dlt))\n", + "\n", + "\n", + "preferences = Preferences(timeout_seconds=1500, optimization_level=3)\n", + "constraints = Constraints(\n", + " optimization_parameter=OptimizationParameter.DEPTH,\n", + " max_width=18, # change max_width as needed, minimum width as per algo is n+2+ log2(D)- but more might be used by synthesizer to do the arithmetics\n", + ")\n", + "qmod = create_model(main, preferences=preferences)\n", + "qmod = set_constraints(qmod, constraints)\n", + "\n", + "# write_qmod(\n", + "# qmod,\n", + "# 'toeplitz_BE_N{}'.format(N)\n", + "# )\n", + "\n", + "# Synthesize the quantum program\n", + "backend_preferences = ClassiqBackendPreferences(backend_name=\"simulator_statevector\")\n", + "qmod = set_execution_preferences(\n", + " qmod,\n", + " execution_preferences=ExecutionPreferences(\n", + " num_shots=1, backend_preferences=backend_preferences\n", + " ),\n", + ")\n", + "qprog = synthesize(qmod)\n", + "# show(qprog)\n", + "circuit_width = qprog.data.width\n", + "circuit_depth = qprog.transpiled_circuit.depth\n", + "print(\"The circuit width is:\", circuit_width)\n", + "print(\"The circuit depth is:\", circuit_depth)\n", + "\n", + "\n", + "job = execute(qprog)\n", + "results = job.result()[0].value\n", + "\n", + "\n", + "reduced_state = get_projected_state_vector(results, \"j\", {\"s\": 0, \"data\": 0, \"dlt\": 0})\n", + "reduced_state = reduced_state * D # D is the subnormalization factor\n", + "print(\"The reduced state vector for j when s=0, dlt=1 and data=0 is:\")\n", + "print(reduced_state)\n", + "\n", + "\n", + "# theoretical reduced state vector\n", + "theoretical_reduced_state = np.matmul(\n", + " toeplitz_matrix, [1 / np.sqrt(N) for i in range(N)]\n", + ")\n", + "print(\"Theoretical reduced state vector for j when l=0 and data=0 is:\")\n", + "print(theoretical_reduced_state)" + ] + }, + { + "cell_type": "markdown", + "id": "16", + "metadata": {}, + "source": [ + "____" + ] + }, + { + "cell_type": "markdown", + "id": "17", + "metadata": {}, + "source": [ + "### 3. Tridiagonal symmetric Matrix" + ] + }, + { + "cell_type": "markdown", + "id": "18", + "metadata": {}, + "source": [ + "We now consider a tridiagonal matrix which is also symmetric. This matrix has the below structure and sparsity pattern.\n", + "\n", + "\n", + "

\n", + " \n", + "

" + ] + }, + { + "cell_type": "markdown", + "id": "19", + "metadata": {}, + "source": [ + "The block encoding technique is similar to the one for Toeplitz matrix. Oracles $O_c$ and $O_r$ need to be appropriately constructed using quantum arithmetics. The oracle $O_{data}$ for loading the unique data values would be identical as above. The below is the full block encoding circuit:\n", + "\n", + "\n", + "

\n", + " \n", + "

\n", + "\n", + "Here the subnormalization factor is $3$ and $3$ ancilla qubits are used irrespective of matrix size." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "20", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Tridiagonal symmetric system [[-0.67 0.91 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]\n", + " [ 0.91 -0.01 0.48 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]\n", + " [ 0. 0.48 -0.85 0.02 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]\n", + " [ 0. 0. 0.02 -0.18 -0.08 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]\n", + " [ 0. 0. 0. -0.08 0.99 -0.2 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]\n", + " [ 0. 0. 0. 0. -0.2 -0.31 -0.61 0. 0. 0. 0. 0. 0. 0. 0. 0. ]\n", + " [ 0. 0. 0. 0. 0. -0.61 0.68 -0.1 0. 0. 0. 0. 0. 0. 0. 0. ]\n", + " [ 0. 0. 0. 0. 0. 0. -0.1 -0.05 -0.63 0. 0. 0. 0. 0. 0. 0. ]\n", + " [ 0. 0. 0. 0. 0. 0. 0. -0.63 -0.5 0.93 0. 0. 0. 0. 0. 0. ]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0. 0.93 -0.65 0.58 0. 0. 0. 0. 0. ]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.58 0.2 -0.31 0. 0. 0. 0. ]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -0.31 -0.48 -0.68 0. 0. 0. ]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -0.68 -0.38 0.99 0. 0. ]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.99 -0.47 -0.15 0. ]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -0.15 -0.44 -0.33]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -0.33 -0.97]]\n", + "Unique values in the matrix: [np.float64(-0.6669), np.float64(0.907), np.float64(-0.01), np.float64(0.4769), np.float64(-0.8479), np.float64(0.0245), np.float64(-0.1823), np.float64(-0.0761), np.float64(0.9882), np.float64(-0.1966), np.float64(-0.3112), np.float64(-0.611), np.float64(0.6835), np.float64(-0.0973), np.float64(-0.0452), np.float64(-0.6281), np.float64(-0.4978), np.float64(0.931), np.float64(-0.6478), np.float64(0.5827), np.float64(0.2046), np.float64(-0.3149), np.float64(-0.4843), np.float64(-0.6751), np.float64(-0.384), np.float64(0.9934), np.float64(-0.4661), np.float64(-0.1534), np.float64(-0.4401), np.float64(-0.3313), np.float64(-0.9661), 0]\n" + ] + } + ], + "source": [ + "# Tridiagonal symmetric matrix\n", + "def generate_tridiagonal_symmetric_matrix(size):\n", + " \"\"\"\n", + " Generate a symmetric tridiagonal matrix of given size.\n", + " The size must be a power of 2, and diagonal values are distinct random values between -1 and 1.\n", + "\n", + " Parameters:\n", + " size (int): Size of the matrix (must be a power of 2).\n", + "\n", + " Returns:\n", + " np.ndarray: A symmetric tridiagonal matrix.\n", + "\n", + " Raises:\n", + " ValueError: If the size is not a power of 2.\n", + " \"\"\"\n", + " if size <= 0 or (size & (size - 1)) != 0:\n", + " raise ValueError(\"Size must be a power of 2.\")\n", + "\n", + " # Generate distinct random diagonal values\n", + " diagonal_values = np.random.uniform(-1, 1, size)\n", + " while len(set(diagonal_values)) < size:\n", + " diagonal_values = np.random.uniform(-1, 1, size) # Ensure uniqueness\n", + "\n", + " # Generate random off-diagonal values between -1 and 1 for tridiagonal matrix\n", + " off_diagonal_values = np.random.uniform(-1, 1, size - 1)\n", + "\n", + " # Construct the tridiagonal symmetric matrix\n", + " matrix = np.zeros((size, size))\n", + " np.fill_diagonal(matrix, diagonal_values)\n", + " np.fill_diagonal(matrix[1:], off_diagonal_values)\n", + " np.fill_diagonal(matrix[:, 1:], off_diagonal_values)\n", + "\n", + " return matrix\n", + "\n", + "\n", + "def unique_val_tridiagonal(matrix):\n", + " \"\"\"\n", + " Extract unique values from a (tri)diagonal matrix in a specific order.\n", + " The value 0 is always placed at the end of the list.\n", + " \"\"\"\n", + " vals = _unique_values_core(matrix, mode=\"tridiagonal\")\n", + " non_zero = [v for v in vals if v != 0]\n", + " return non_zero + [0]\n", + "\n", + "\n", + "# Test case for tridiagonal symmetric matrix, let's take a 16x16 matrix\n", + "size = 16 # Must be a power of 2\n", + "matrix = generate_tridiagonal_symmetric_matrix(size)\n", + "matrix = np.round(matrix, 4)\n", + "print(\"Tridiagonal symmetric system\", matrix)\n", + "unique_values = unique_val_tridiagonal(\n", + " matrix\n", + ") # unique values in the matrix (each has a unique label- d)\n", + "\n", + "print(\"Unique values in the matrix:\", unique_values)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "21", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The circuit width is: 10\n", + "The circuit depth is: 1382\n", + "The reduced state vector for j when s=0, data=0 is:\n", + "[ 0.02 0.11 -0.03 -0.02 0.06 -0.09 -0. -0.06 -0.02 0.07 0.04 -0.12 -0.01 0.03 -0.08 -0.11]\n", + "Expected theoretical state vector:\n", + "[ 0.02 0.11 -0.03 -0.02 0.06 -0.09 -0. -0.06 -0.02 0.07 0.04 -0.12 -0.01 0.03 -0.08 -0.11]\n", + "********\n", + "The reduced state vector for j when s=0, data=0 is:\n", + "[ 0.02 0.11 -0.03 -0.02 0.06 -0.09 -0. -0.06 -0.02 0.07 0.04 -0.12 -0.01 0.03 -0.08 -0.11]\n", + "Expected theoretical state vector:\n", + "[ 0.02 0.11 -0.03 -0.02 0.06 -0.09 -0. -0.06 -0.02 0.07 0.04 -0.12 -0.01 0.03 -0.08 -0.11]\n", + "********\n" + ] + } + ], + "source": [ + "@qfunc\n", + "def eq_sup(s: QNum):\n", + " inplace_prepare_amplitudes(\n", + " [np.sqrt(1 / 3), np.sqrt(1 / 3), 0, np.sqrt(1 / 3)], 0, s\n", + " )\n", + "\n", + "\n", + "@qfunc\n", + "def main(s: Output[QNum], j: Output[QNum], data: Output[QNum]):\n", + " n = int(math.log2(size))\n", + " allocate(2, s)\n", + " allocate(n, j)\n", + " allocate(1, data)\n", + "\n", + " apply_to_all(H, j)\n", + " eq_sup(s)\n", + " T_M = unique_values\n", + "\n", + " d = QNum(\"d\", n + 1, False, 0)\n", + " s0 = QNum(\"s0\", 1, False, 0)\n", + " s1 = QNum(\"s1\", 1, False, 0)\n", + " bind(s, [s0, s1])\n", + " control(s0 == 1, lambda: control(s1 == 0, lambda: inplace_add(-1, j)))\n", + "\n", + " bind([s0, j], d)\n", + " oracle_d(d, data, T_M)\n", + " bind(d, [s0, j])\n", + "\n", + " control(s1 == 1, lambda: inplace_add(1, j))\n", + " bind([s0, s1], s)\n", + "\n", + " invert(lambda: eq_sup(s))\n", + "\n", + "\n", + "backend_preferences = ClassiqBackendPreferences(backend_name=\"simulator_statevector\")\n", + "constraints = Constraints(\n", + " optimization_parameter=OptimizationParameter.DEPTH,\n", + " max_width=10, # change max_width as needed, minimum width is n+3\n", + ")\n", + "preferences = Preferences(timeout_seconds=900)\n", + "\n", + "qmod = create_model(main)\n", + "qmod = set_execution_preferences(\n", + " qmod,\n", + " execution_preferences=ExecutionPreferences(\n", + " num_shots=100000, backend_preferences=backend_preferences\n", + " ),\n", + ")\n", + "qmod = set_preferences(qmod, preferences)\n", + "qmod = set_constraints(qmod, constraints)\n", + "# write_qmod(\n", + "# qmod,\n", + "# 'tridiagonal_BE_size{}'.format(size)\n", + "# )\n", + "\n", + "qprog = synthesize(qmod)\n", + "# show(qprog)\n", + "circuit_width = qprog.data.width\n", + "circuit_depth = qprog.transpiled_circuit.depth\n", + "print(\"The circuit width is:\", circuit_width)\n", + "print(\"The circuit depth is:\", circuit_depth)\n", + "\n", + "job = execute(qprog)\n", + "results = job.result()[0].value\n", + "\n", + "reduced_state = get_projected_state_vector(results, \"j\", {\"s\": 0, \"data\": 0})\n", + "print(\"The reduced state vector for j when s=0, data=0 is:\")\n", + "print(reduced_state)\n", + "\n", + "theoretical_reduced_state = (\n", + " np.matmul(matrix, [1 / np.sqrt(size) for i in range(size)]) / 3\n", + ") # here 3 is the subnormalization factor\n", + "print(f\"Expected theoretical state vector:\")\n", + "print(theoretical_reduced_state)\n", + "print(\"********\")" + ] + }, + { + "cell_type": "markdown", + "id": "22", + "metadata": {}, + "source": [ + "___" + ] + }, + { + "cell_type": "markdown", + "id": "23", + "metadata": {}, + "source": [ + "### 4. Two-dimensional Laplacian\n", + "\n", + "The discrete Laplacian in one dimension is obtained from the second–order finite difference stencil \n", + "$$\n", + "\\Delta f(x_i) = \\frac{f(x_{i-1}) - 2 f(x_i) + f(x_{i+1})}{(\\Delta x)^2},\n", + "$$\n", + "which corresponds to a tridiagonal Toeplitz matrix with $-2/(\\Delta x)^2$ on the diagonal and $1/(\\Delta x)^2$ on the off-diagonals.\n", + "\n", + "To generalize this to two dimensions, consider a rectangular grid of size \n", + "$$\n", + "N_x \\times N_y,\n", + "$$\n", + "with grid spacings $\\Delta x$ and $\\Delta y$. \n", + "A grid point is indexed as $(a,b)$ with \n", + "$$\n", + "a = 0,\\ldots,N_x-1, \\qquad b = 0,\\ldots,N_y-1.\n", + "$$\n", + "\n", + "The 2D finite-difference Laplacian uses the standard 5-point stencil:\n", + "$$\n", + "\\Delta f(x_a, y_b)\n", + "=\n", + "\\frac{f(x_{a-1},y_b) - 2 f(x_a,y_b) + f(x_{a+1},y_b)}{(\\Delta x)^2}\n", + "+\n", + "\\frac{f(x_a,y_{b-1}) - 2 f(x_a,y_b) + f(x_a,y_{b+1})}{(\\Delta y)^2}.\n", + "$$\n", + "\n", + "To write this as a matrix, we flatten the 2D grid into a vector of length \n", + "$$\n", + "N = N_x N_y\n", + "$$\n", + "using **row-major ordering**:\n", + "$$\n", + "f_{a + b N_x} = f(x_a, y_b).\n", + "$$\n", + "\n", + "With this ordering, the Laplacian becomes an $N \\times N$ sparse matrix $A$. \n", + "Each row corresponds to a grid point and couples only to its nearest neighbors. \n", + "The matrix entries take only a few values:\n", + "$$\n", + "A_{a_1 + b_1 N_x,\\; a_2 + b_2 N_x} =\n", + "\\begin{cases}\n", + "A_0 = -2\\!\\left(\\frac{1}{(\\Delta x)^2} + \\frac{1}{(\\Delta y)^2}\\right), & a_1=a_2,\\; b_1=b_2, \\\\[6pt]\n", + "A_1 = \\frac{1}{(\\Delta x)^2}, & |a_1-a_2| = 1,\\; b_1=b_2, \\\\[6pt]\n", + "A_2 = \\frac{1}{(\\Delta y)^2}, & |b_1-b_2| = 1,\\; a_1=a_2, \\\\[6pt]\n", + "0, & \\text{otherwise}.\n", + "\\end{cases}\n", + "$$\n", + "\n", + "This produces the familiar **five-diagonal block structure**: \n", + "- the main diagonal contains $A_0$, \n", + "- horizontal neighbors contribute $A_1$, \n", + "- vertical neighbors contribute $A_2$, \n", + "- all other entries are zero.\n", + "\n", + "For example, when $N_x = N_y = 4$, we get the below structure:\n", + "\n", + "\n", + "

\n", + " \n", + "

\n", + "\n", + "\n", + "Below is the full block encoding circuit (which is also Hermitian in this case):\n", + "\n", + "

\n", + " \n", + "

\n", + "\n", + "Here the subnormalization factor is $5$ and $5$ ancilla qubits are used, irrespective of matrix size. \n", + "\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "24", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2D Laplacian Matrix:\n", + " [[-4. 1. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [ 1. -4. 1. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [ 0. 1. -4. 1. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [ 0. 0. 1. -4. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [ 1. 0. 0. 0. -4. 1. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0.]\n", + " [ 0. 1. 0. 0. 1. -4. 1. 0. 0. 1. 0. 0. 0. 0. 0. 0.]\n", + " [ 0. 0. 1. 0. 0. 1. -4. 1. 0. 0. 1. 0. 0. 0. 0. 0.]\n", + " [ 0. 0. 0. 1. 0. 0. 1. -4. 0. 0. 0. 1. 0. 0. 0. 0.]\n", + " [ 0. 0. 0. 0. 1. 0. 0. 0. -4. 1. 0. 0. 1. 0. 0. 0.]\n", + " [ 0. 0. 0. 0. 0. 1. 0. 0. 1. -4. 1. 0. 0. 1. 0. 0.]\n", + " [ 0. 0. 0. 0. 0. 0. 1. 0. 0. 1. -4. 1. 0. 0. 1. 0.]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 1. -4. 0. 0. 0. 1.]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. -4. 1. 0. 0.]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 1. -4. 1. 0.]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 1. -4. 1.]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 1. -4.]]\n" + ] + } + ], + "source": [ + "# Laplacian Matrix\n", + "def generate_laplacian_matrix(Nx, Ny, delta_x, delta_y):\n", + " N = Nx * Ny\n", + " A = np.zeros((N, N))\n", + "\n", + " A0 = -2 / (delta_x**2) - 2 / (delta_y**2)\n", + " A1 = 1 / (delta_x**2)\n", + " A2 = 1 / (delta_y**2)\n", + "\n", + " for a1 in range(Nx):\n", + " for b1 in range(Ny):\n", + " index = a1 + b1 * Nx\n", + "\n", + " A[index, index] = A0\n", + "\n", + " # off-diagonal in x-direction\n", + " if a1 > 0:\n", + " A[index, index - 1] = A1\n", + " if a1 < Nx - 1:\n", + " A[index, index + 1] = A1\n", + "\n", + " # off-diagonal in y-direction\n", + " if b1 > 0:\n", + " A[index, index - Nx] = A2\n", + " if b1 < Ny - 1:\n", + " A[index, index + Nx] = A2\n", + "\n", + " return A\n", + "\n", + "\n", + "def unique_val_lap(laplacian, Nx):\n", + " \"\"\"Extract unique diagonal values from a 2D Laplacian matrix.\"\"\"\n", + " return _unique_values_core(laplacian, mode=\"laplacian\", Nx=Nx)\n", + "\n", + "\n", + "# Test case for 2D Laplacian matrix\n", + "Nx = 4 # Number of grid points in x-direction (must be a power of 2)\n", + "Ny = 4 # Number of grid points in y-direction (must be a power of 2)\n", + "delta_x = 1.0 # Grid spacing in x-direction\n", + "delta_y = 1.0 # Grid spacing in y-direction\n", + "laplacian = generate_laplacian_matrix(Nx, Ny, delta_x, delta_y)\n", + "print(\"2D Laplacian Matrix:\\n\", laplacian)\n", + "lap_scaled = laplacian / np.max(np.abs(laplacian))\n", + "\n", + "unique_values_laplacian = unique_val_lap(lap_scaled, Nx)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "25", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The circuit width is: 10\n", + "The circuit depth is: 324\n", + "The reduced state vector for j when s=0, data=0 is:\n", + "[0.03 0.01 0.01 0.03 0.01 0. 0. 0.01 0.01 0. 0. 0.01 0.03 0.01 0.01 0.03]\n", + "Expected theoretical state vector:\n", + "[-0.03 -0.01 -0.01 -0.03 -0.01 0. 0. -0.01 -0.01 0. 0. -0.01 -0.03 -0.01 -0.01 -0.03]\n", + "********\n" + ] + } + ], + "source": [ + "j1 = int(math.log2(Nx))\n", + "j2 = int(math.log2(Ny))\n", + "\n", + "\n", + "@qfunc\n", + "def eq_sup_lap(s: QNum):\n", + " inplace_prepare_state([1 / 5, 0, 1 / 5, 1 / 5, 1 / 5, 1 / 5, 0, 0], 0, s)\n", + "\n", + "\n", + "@qfunc\n", + "def oracle_org(s: QArray, j: QNum, dlt: QBit):\n", + " j1 = math.log2(Nx)\n", + " j2 = math.log2(Ny)\n", + " j1reg = QNum(\"j1\", j1, False, 0)\n", + " j2reg = QNum(\"j2\", j2, False, 0)\n", + " s0 = QNum(\"s0\", 1, False, 0)\n", + " s1 = QNum(\"s1\", 1, False, 0)\n", + " s2 = QNum(\"s2\", 1, False, 0)\n", + " bind(s, [s0, s1, s2])\n", + " bind(j, [j1reg, j2reg])\n", + " control(s0 == 1, lambda: control(s1 == 0, lambda: control(s2 == 0, lambda: X(dlt))))\n", + " control(s1 == 1, lambda: control(j1reg == 0, lambda: X(dlt)))\n", + " control(s2 == 1, lambda: control(j2reg == 0, lambda: X(dlt)))\n", + " bind([j1reg, j2reg], j)\n", + " bind([s0, s1, s2], s)\n", + "\n", + "\n", + "@qfunc\n", + "def laplacian_BE(s: QNum, j: QNum, data: QNum, dlt: QNum):\n", + " j1reg = QNum(\"j1\", j1, False, 0)\n", + " j2reg = QNum(\"j2\", j2, False, 0)\n", + " s0 = QNum(\"s0\", 1, False, 0)\n", + " s1 = QNum(\"s1\", 1, False, 0)\n", + " s2 = QNum(\"s2\", 1, False, 0)\n", + " bind(s, [s0, s1, s2])\n", + " bind(j, [j1reg, j2reg])\n", + " control(s0 == 0, lambda: control(s2 == 1, lambda: inplace_add(1, j2reg)))\n", + " bind([j1reg, j2reg], j)\n", + " control(s1 == 1, lambda: control(s0 == 0, lambda: inplace_add(1, j)))\n", + " bind([s0, s1, s2], s)\n", + "\n", + " oracle_org(s, j, dlt)\n", + "\n", + " bind(s, [s0, s1, s2])\n", + " d = QNum(\"d\", 2, False, 0)\n", + " bind([s1, s2], d)\n", + " Z(data)\n", + " oracle_d(d, data, unique_values_laplacian)\n", + " bind(d, [s1, s2])\n", + "\n", + " control(s1 == 1, lambda: X(s0))\n", + " control(s2 == 1, lambda: X(s0))\n", + "\n", + " control(\n", + " s1 == 1, lambda: control(s0 == 0, lambda: inplace_add(-1, j))\n", + " ) # is j=0 here causing an issue? problem here when superposition used\n", + " bind(j, [j1reg, j2reg])\n", + " control(\n", + " s2 == 1, lambda: control(s0 == 0, lambda: inplace_add(-1, j2reg))\n", + " ) # is j=0 here causing an issue?\n", + " bind([j1reg, j2reg], j)\n", + " bind([s0, s1, s2], s)\n", + "\n", + "\n", + "@qfunc\n", + "def main(s: Output[QNum], j: Output[QNum], data: Output[QNum], dlt: Output[QBit]):\n", + " allocate(3, s)\n", + " allocate(1, dlt)\n", + " allocate(1, data)\n", + " allocate(j1 + j2, j)\n", + " apply_to_all(H, j)\n", + "\n", + " within_apply(lambda: eq_sup_lap(s), lambda: laplacian_BE(s, j, data, dlt))\n", + "\n", + "\n", + "# print(states)\n", + "qmod = create_model(main)\n", + "backend_preferences = ClassiqBackendPreferences(backend_name=\"simulator_statevector\")\n", + "qmod = set_execution_preferences(\n", + " qmod,\n", + " execution_preferences=ExecutionPreferences(\n", + " num_shots=1000000, backend_preferences=backend_preferences\n", + " ),\n", + ")\n", + "constraints = Constraints(\n", + " optimization_parameter=OptimizationParameter.DEPTH,\n", + " max_width=10, # change max_width as needed, minimum width is j1+j2+5\n", + ")\n", + "qmod = set_constraints(qmod, constraints)\n", + "\n", + "# write_qmod(\n", + "# qmod,\n", + "# 'laplacian_BE_Nx{}_Ny{}'.format(Nx, Ny)\n", + "# )\n", + "\n", + "qprog = synthesize(qmod)\n", + "# show(qprog)\n", + "circuit_width = qprog.data.width\n", + "circuit_depth = qprog.transpiled_circuit.depth\n", + "print(\"The circuit width is:\", circuit_width)\n", + "print(\"The circuit depth is:\", circuit_depth)\n", + "\n", + "job = execute(qprog)\n", + "results = job.result()[0].value\n", + "\n", + "# Post processing\n", + "reduced_state = get_projected_state_vector(results, \"j\", {\"s\": 0, \"data\": 0, \"dlt\": 0})\n", + "print(\"The reduced state vector for j when s=0, data=0 is:\")\n", + "print(reduced_state)\n", + "\n", + "# theoretical reduced state vector\n", + "theoretical_reduced_state = (\n", + " np.matmul(lap_scaled, [1 / np.sqrt(Nx * Ny) for i in range(Nx * Ny)]) / 5\n", + ")\n", + "print(f\"Expected theoretical state vector:\")\n", + "print(theoretical_reduced_state)\n", + "print(\"********\")" + ] + }, + { + "cell_type": "markdown", + "id": "26", + "metadata": {}, + "source": [ + "____" + ] + }, + { + "cell_type": "markdown", + "id": "27", + "metadata": {}, + "source": [ + "### Reference\n", + "\n", + "- Sunderhauf, C., Campbell, E., Camps, J.: Block-encoding structured matrices fordata input in quantum computing. Quantum 8, 1226 (2024) https://doi.org/10.22331/q-2024-01-11-1226." + ] + }, + { + "cell_type": "markdown", + "id": "28", + "metadata": {}, + "source": [ + "___" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "qiskit-fresh", + "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.10.18" + } + }, + "nbformat": 4, + "nbformat_minor": 9 +} diff --git a/applications/block_encoding/symmetric_tridiagonal.png b/applications/block_encoding/symmetric_tridiagonal.png new file mode 100644 index 000000000..4fc96529b Binary files /dev/null and b/applications/block_encoding/symmetric_tridiagonal.png differ diff --git a/applications/block_encoding/toeplitz_BE_N16.qmod b/applications/block_encoding/toeplitz_BE_N16.qmod new file mode 100644 index 000000000..431b60594 --- /dev/null +++ b/applications/block_encoding/toeplitz_BE_N16.qmod @@ -0,0 +1,42 @@ +qfunc oracle_d(s: qnum, data: qnum, TM: real[]) { + repeat (i: TM::len) { + control (s == i) { + RX(2 * acos(TM[i]), data); + } + } +} + +qfunc oracle_r(j: qnum, s: qnum, k: int) { + j += s - k; +} + +qfunc Toep_BE(data: qnum, s: qnum, j: qnum, dlt: qbit) { + oracle_d(s, data, [ + 0.6406, + 0.6246, + (-0.3873), + 0.952 + ]); + tmp: qnum; + tmp = (s - 1) + j; + oracle_r(j, s, 1); + control (tmp < 0) { + X(dlt); + } + control (tmp >= 16) { + X(dlt); + } +} + +qfunc main(output data: qnum, output s: qnum, output j: qnum, output dlt: qbit) { + allocate(2, s); + allocate(4, j); + allocate(1, data); + allocate(1, dlt); + apply_to_all(H, j); + within { + apply_to_all(H, s); + } apply { + Toep_BE(data, s, j, dlt); + } +} diff --git a/applications/block_encoding/toeplitz_BE_N16.synthesis_options.json b/applications/block_encoding/toeplitz_BE_N16.synthesis_options.json new file mode 100644 index 000000000..1815d2e1c --- /dev/null +++ b/applications/block_encoding/toeplitz_BE_N16.synthesis_options.json @@ -0,0 +1,45 @@ +{ + "constraints": { + "max_gate_count": {}, + "max_width": 18, + "optimization_parameter": "depth" + }, + "preferences": { + "custom_hardware_settings": { + "basis_gates": [ + "z", + "rz", + "u", + "u2", + "sx", + "cx", + "cz", + "cy", + "y", + "tdg", + "r", + "rx", + "sdg", + "t", + "x", + "s", + "u1", + "ry", + "p", + "id", + "sxdg", + "h" + ], + "is_symmetric_connectivity": true + }, + "debug_mode": true, + "machine_precision": 8, + "optimization_level": 3, + "output_format": ["qasm"], + "pretty_qasm": true, + "random_seed": 2157362051, + "synthesize_all_separately": false, + "timeout_seconds": 1500, + "transpilation_option": "auto optimize" + } +} diff --git a/applications/block_encoding/toeplitz_circuit.png b/applications/block_encoding/toeplitz_circuit.png new file mode 100644 index 000000000..95c005c1f Binary files /dev/null and b/applications/block_encoding/toeplitz_circuit.png differ diff --git a/applications/block_encoding/tridiagonal_BE_size16.qmod b/applications/block_encoding/tridiagonal_BE_size16.qmod new file mode 100644 index 000000000..b9020adb4 --- /dev/null +++ b/applications/block_encoding/tridiagonal_BE_size16.qmod @@ -0,0 +1,71 @@ +qfunc eq_sup(s: qnum) { + inplace_prepare_amplitudes([0.5774, 0.5774, 0, 0.5774], 0, s); +} + +qfunc oracle_d(s: qnum, data: qnum, TM: real[]) { + repeat (i: TM::len) { + control (s == i) { + RX(2 * acos(TM[i]), data); + } + } +} + +qfunc main(output s: qnum, output j: qnum, output data: qnum) { + allocate(2, s); + allocate(4, j); + allocate(1, data); + apply_to_all(H, j); + eq_sup(s); + d: qnum<5, False, 0>; + s0: qnum<1, False, 0>; + s1: qnum<1, False, 0>; + s -> {s0, s1}; + control (s0 == 1) { + control (s1 == 0) { + j += -1; + } + } + {s0, j} -> d; + oracle_d(d, data, [ + 0.8015, + (-0.8087), + 0.6652, + (-0.1181), + (-0.2983), + 0.6675, + 0.8014, + 0.6437, + 0.7776, + 0.3898, + (-0.7502), + (-0.2429), + 0.683, + 0.2305, + 0.1147, + 0.6802, + 0.7784, + 0.3091, + (-0.2962), + (-0.2641), + 0.7889, + 0.3176, + 0.1564, + (-0.4584), + (-0.3464), + 0.5023, + (-0.1425), + (-0.8012), + 0.811, + 0.0217, + 0.465, + 0 + ]); + d -> {s0, j}; + control (s1 == 1) { + j += 1; + } + {s0, s1} -> s; + invert { + eq_sup(s); + } +} diff --git a/applications/block_encoding/tridiagonal_BE_size16.synthesis_options.json b/applications/block_encoding/tridiagonal_BE_size16.synthesis_options.json new file mode 100644 index 000000000..8fc9f1360 --- /dev/null +++ b/applications/block_encoding/tridiagonal_BE_size16.synthesis_options.json @@ -0,0 +1,45 @@ +{ + "constraints": { + "max_gate_count": {}, + "max_width": 10, + "optimization_parameter": "depth" + }, + "preferences": { + "custom_hardware_settings": { + "basis_gates": [ + "z", + "rz", + "u", + "u2", + "sx", + "cx", + "cz", + "cy", + "y", + "tdg", + "r", + "rx", + "sdg", + "t", + "x", + "s", + "u1", + "ry", + "p", + "id", + "sxdg", + "h" + ], + "is_symmetric_connectivity": true + }, + "debug_mode": true, + "machine_precision": 8, + "optimization_level": 1, + "output_format": ["qasm"], + "pretty_qasm": true, + "random_seed": 2292430670, + "synthesize_all_separately": false, + "timeout_seconds": 900, + "transpilation_option": "auto optimize" + } +} diff --git a/applications/block_encoding/tridiagonal_symmetric_circuit.png b/applications/block_encoding/tridiagonal_symmetric_circuit.png new file mode 100644 index 000000000..e71c0a072 Binary files /dev/null and b/applications/block_encoding/tridiagonal_symmetric_circuit.png differ diff --git a/tests/notebooks/test_ND_Laplacian_BE.py b/tests/notebooks/test_ND_Laplacian_BE.py new file mode 100644 index 000000000..bf96ff3fd --- /dev/null +++ b/tests/notebooks/test_ND_Laplacian_BE.py @@ -0,0 +1,21 @@ +from tests.utils_for_testbook import ( + validate_quantum_program_size, + validate_quantum_model, + wrap_testbook, +) +from testbook.client import TestbookNotebookClient + + +@wrap_testbook("ND_Laplacian_BE", timeout_seconds=300) +def test_notebook(tb: TestbookNotebookClient) -> None: + # test models + validate_quantum_model(tb.ref("qmod")) + # test quantum programs + validate_quantum_program_size( + tb.ref_pydantic("qprog"), + expected_width=None, + expected_depth=None, + ) + + # test notebook content + pass # Todo diff --git a/tests/notebooks/test_select_structures_BE.py b/tests/notebooks/test_select_structures_BE.py new file mode 100644 index 000000000..9e0f31441 --- /dev/null +++ b/tests/notebooks/test_select_structures_BE.py @@ -0,0 +1,21 @@ +from tests.utils_for_testbook import ( + validate_quantum_program_size, + validate_quantum_model, + wrap_testbook, +) +from testbook.client import TestbookNotebookClient + + +@wrap_testbook("select_structures_BE", timeout_seconds=600) +def test_notebook(tb: TestbookNotebookClient) -> None: + # test models + validate_quantum_model(tb.ref("qmod")) + # test quantum programs + validate_quantum_program_size( + tb.ref_pydantic("qprog"), + expected_width=None, + expected_depth=None, + ) + + # test notebook content + pass # Todo diff --git a/tests/resources/timeouts.yaml b/tests/resources/timeouts.yaml index 85a239109..4252d75e7 100644 --- a/tests/resources/timeouts.yaml +++ b/tests/resources/timeouts.yaml @@ -1,8 +1,11 @@ +1D_Periodic_Laplacian_BE.qmod: 30 3_sat_grover_large.qmod: 60 3_sat_grover_small.qmod: 90 3sat_oracles.ipynb: 1800 CRX.qmod: 10 CX.qmod: 10 +ND_Laplacian_BE.ipynb: 300 +ND_Periodic_Laplacian_BE.qmod: 10 Option_Pricing_Workshop.ipynb: 400 PHASE.qmod: 10 Qmod_tutorial_part1.ipynb: 500 @@ -53,6 +56,7 @@ bpde.qmod: 10 cheb_lcu_solver_banded_be.qmod: 400 cheb_lcu_solver_pauli_be.qmod: 400 chebyshev_approximation.ipynb: 20 +checkerboard_BE_N8.qmod: 200 classiq_chemistry_application.ipynb: 80 classiq_discrete_quantum_walk.ipynb: 300 classiq_iQuHack_2025_final.ipynb: 200 @@ -154,6 +158,7 @@ ising_model.ipynb: 300 ising_model.qmod: 300 kidney_exchange_problem.ipynb: 300 kidney_exchange_problem.qmod: 300 +laplacian_BE_Nx4_Ny4.qmod: 200 learning_optimization.ipynb: 80 linear_combination_of_unitaries.ipynb: 20 linear_combination_of_unitaries.qmod: 10 @@ -304,6 +309,7 @@ rectangles_packing.qmod: 1800 rectangles_packing_grid.ipynb: 1800 second_quantized_hamiltonian.ipynb: 44 second_quantized_hamiltonian.qmod: 40 +select_structures_BE.ipynb: 20 set_cover.ipynb: 1440 set_cover.qmod: 1216 set_partition.ipynb: 350 @@ -337,8 +343,10 @@ time_marching.ipynb: 500 time_marching.qmod: 400 tket_discrete_quantum_walk.ipynb: 300 tket_qsvt_example.ipynb: 300 +toeplitz_BE_N16.qmod: 200 traveling_salesman_problem.ipynb: 1068 traveling_salesman_problem.qmod: 1068 +tridiagonal_BE_size16.qmod: 200 unitary.ipynb: 20 unitary.qmod: 10 variational_data_encoding.ipynb: 20