![](\"images/quantum_enhanced_optimization_LABS/radar.png\")
\n",
+ "\n",
+ "\n",
+ "So, how do you select a good signal? This is where LABS comes in, defining these questions as a binary optimization problem. Given a binary sequence of length $N$, $(s_1 \\cdots s_N) \\in {\\pm 1}^N$, the goal is to minimize the following objective function.\n",
+ "\n",
+ "$$ E(s) = \\sum_{k=1}^{N-1} C_k^2 $$\n",
+ "\n",
+ "Where $C_k$ is defined as. \n",
+ "\n",
+ " $$C_k= \\sum_{i=1}^{N-k} s_is_{i+k}$$\n",
+ "\n",
+ "\n",
+ "So, each $C_k$ computes how similar the original signal is to the shifted one for each offset value $k$. To explore this more, try the interactive widget linked [here](https://nvidia.github.io/cuda-q-academic/interactive_widgets/labs_visualization.html). See if you can select a very good and very poor sequence and see the difference for the case of $N=7$."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "84fa6dff-0fee-4251-a006-76d6ad426116",
+ "metadata": {},
+ "source": [
+ "## Classical Solution of the LABS problem\n",
+ "\n",
+ "The LABS problem is tricky to solve for a few reasons. First, the configuration space grows exponentially. Second, underlying symmetries of the problem result in many degeneracies in the optimization landscape severely inhibiting local search methods. \n",
+ "\n",
+ "![](\"images/quantum_enhanced_optimization_LABS/mts_algorithm.png\")
\n",
+ "\n",
+ "Such an approach is fast, parallelizable, and allows for exploration with improved searching of the solution landscape. \n",
+ "\n",
+ "![](\"images/quantum_enhanced_optimization_LABS/combine_mutate.png\")
\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "id": "341c9a3f-565c-49ab-9903-0aa9a508722c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "#TODO - Write code to perform MTS\n",
+ "\n",
+ "import cudaq\n",
+ "import numpy as np\n",
+ "from math import floor\n",
+ "import auxiliary_files.labs_utils as utils\n",
+ "\n",
+ "\n",
+ "import random \n",
+ "\n",
+ "#creates bitstrings of length N\n",
+ "def bitstring(N):\n",
+ " string = []\n",
+ " for bit in range(N):\n",
+ " string.append(random.choice([0,1])) \n",
+ " return string\n",
+ "\n",
+ "\n",
+ "\n",
+ "#creates s_max bitstrings of lenth N\n",
+ "def populate(s_max,N):\n",
+ " pop = []\n",
+ " for i in range(s_max):\n",
+ " pop.append(bitstring(N))\n",
+ " return pop\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "#makes a array from 0 to N\n",
+ "def cut(N):\n",
+ " k_list = []\n",
+ " for bit in range(N):\n",
+ " k_list.append(bit) \n",
+ " return k_list\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "#cuts two parent bitstrings p1 and p2 of length N at index k and combines the halves into a new bitstring child of length N\n",
+ "def combine(p1,p2):\n",
+ " if len(p1)==len(p2):\n",
+ " k = random.choice(cut(len(p1)))\n",
+ " child = p1[0:k]\n",
+ " for element in p2[k:len(p2)]:\n",
+ " child.append(element)\n",
+ " return child\n",
+ "\n",
+ "\n",
+ "\n",
+ "#pm probability of flippinng each bit in a bitstring s \n",
+ "def mutate(s,pm):\n",
+ " for i in range(len(s)):\n",
+ " if random.random()\n",
+ "
\n",
+ "\n",
+ "Exercise 3:
\n", + "\n", + "At first glance, this equation might looks quite complicated. However, observe the structure and note two \"blocks\" of terms. Can you spot them? \n", + "\n", + "They are 2 qubit terms that look like $R_{YZ}(\\theta)$ or $R_{ZY}(\\theta)$.\n", + "\n", + "As well as 4 qubit terms that look like $R_{YZZZ}(\\theta)$, $R_{ZYZZ}(\\theta)$, $R_{ZZYZ}(\\theta)$, or $R_{ZZZY}(\\theta)$.\n", + "\n", + "Thankfully the authors derive a pair of circuit implementations for the two and four qubit terms respectively, shown in the figures below.\n", + "\n", + "Using CUDA-Q, write a kernel for each which will be used later to construct the full implementation.\n", + "\n", + "* Hint: Remember that the adjoint of a rotation gate is the same as rotating in the opposite direction. \n", + "\n", + "* Hint: You may also want to define a CUDA-Q kernel for an R$_{ZZ}$ gate.\n", + "\n", + "* Hint: Implementing a circuit from a paper is a great place where AI can help accelerate your work. If you have access to a coding assistant, feel free to use it here.\n", + "
![](\"images/quantum_enhanced_optimization_LABS/kernels.png\")
\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "id": "c91bbaae-62a1-41e7-a285-af885627a942",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ " ╭────────────╮╭───╮ ╭────────────╮╭───╮╭───────────╮ ╭───╮ »\n",
+ "q0 : ┤ rx(-1.571) ├┤ h ├──●──┤ rz(-1.571) ├┤ x ├┤ rx(1.571) ├─────┤ h ├─────»\n",
+ " ├───────────┬╯├───┤╭─┴─╮├────────────┤╰─┬─╯├───────────┴╮╭───┴───┴────╮»\n",
+ "q1 : ┤ ry(1.571) ├─┤ h ├┤ x ├┤ rz(-1.571) ├──●──┤ ry(-1.571) ├┤ rx(-1.571) ├»\n",
+ " ├───────────┴╮├───┤╰───╯├────────────┤ ├───────────┬╯├───────────┬╯»\n",
+ "q2 : ┤ ry(-1.571) ├┤ h ├──●──┤ rz(-1.571) ├──●──┤ ry(1.571) ├─┤ rx(1.571) ├─»\n",
+ " ╰───┬───┬────╯╰───╯╭─┴─╮├────────────┤╭─┴─╮├───────────┴╮╰───┬───┬───╯ »\n",
+ "q3 : ────┤ h ├──────────┤ x ├┤ rz(-1.571) ├┤ x ├┤ rx(-1.571) ├────┤ h ├─────»\n",
+ " ╰───╯ ╰───╯╰────────────╯╰───╯╰────────────╯ ╰───╯ »\n",
+ "\n",
+ "################################################################################\n",
+ "\n",
+ " ╭───────────╮╭───╮»\n",
+ "─────────────────────────────────────────────────────────●──┤ rz(1.571) ├┤ x ├»\n",
+ "╭───╮ ╭───────╮╭───╮╭───────────╮╭───────────╮╭───╮╭─┴─╮├───────────┤╰─┬─╯»\n",
+ "┤ h ├──●──┤ rz(0) ├┤ x ├┤ rx(1.571) ├┤ ry(1.571) ├┤ h ├┤ x ├┤ rz(1.571) ├──●──»\n",
+ "├───┤╭─┴─╮├───────┤╰─┬─╯├───────────┤╰───┬───┬───╯╰───╯╰───╯╰───────────╯ »\n",
+ "┤ h ├┤ x ├┤ rz(0) ├──●──┤ rx(3.142) ├────┤ h ├────────────────────────────────»\n",
+ "╰───╯╰───╯╰───────╯ ╰───────────╯ ╰───╯ »\n",
+ "──────────────────────────────────────────────────────────────────────────────»\n",
+ " »\n",
+ "\n",
+ "################################################################################\n",
+ "\n",
+ "╭───────────╮ ╭───╮ »\n",
+ "┤ rx(1.571) ├─┤ h ├─────────────────────────────────────────────────────────»\n",
+ "├───────────┴╮├───┤ ╭───────╮╭───╮╭────────────╮╭───╮ »\n",
+ "┤ ry(-1.571) ├┤ h ├──●──┤ rz(0) ├┤ x ├┤ rx(-3.142) ├┤ h ├───────────────────»\n",
+ "╰────────────╯╰───╯╭─┴─╮├───────┤╰─┬─╯├───────────┬╯├───┤ ╭────────────╮»\n",
+ "───────────────────┤ x ├┤ rz(0) ├──●──┤ ry(1.571) ├─┤ h ├──●──┤ rz(-1.571) ├»\n",
+ " ╰───╯╰───────╯ ╰───────────╯ ╰───╯╭─┴─╮├────────────┤»\n",
+ "─────────────────────────────────────────────────────────┤ x ├┤ rz(-1.571) ├»\n",
+ " ╰───╯╰────────────╯»\n",
+ "\n",
+ "################################################################################\n",
+ "\n",
+ " »\n",
+ "──────────────────────────────────────────────────────────────────────»\n",
+ " ╭───────╮╭───╮╭───────────╮»\n",
+ "────────────────────────────────────────●──┤ rz(0) ├┤ x ├┤ rx(1.571) ├»\n",
+ "╭───╮╭────────────╮╭────────────╮╭───╮╭─┴─╮├───────┤╰─┬─╯├───────────┤»\n",
+ "┤ x ├┤ ry(-1.571) ├┤ rx(-1.571) ├┤ h ├┤ x ├┤ rz(0) ├──●──┤ rx(1.571) ├»\n",
+ "╰─┬─╯├────────────┤╰───┬───┬────╯╰───╯╰───╯╰───────╯ ╰───────────╯»\n",
+ "──●──┤ rx(-1.571) ├────┤ h ├──────────────────────────────────────────»\n",
+ " ╰────────────╯ ╰───╯ »\n",
+ "\n",
+ "################################################################################\n",
+ "\n",
+ " ╭───────────╮╭───╮ \n",
+ "─────────────────────●──┤ rz(1.571) ├┤ x ├──────────────\n",
+ "╭────────────╮╭───╮╭─┴─╮├───────────┤╰─┬─╯╭───────────╮ \n",
+ "┤ ry(-1.571) ├┤ h ├┤ x ├┤ rz(1.571) ├──●──┤ ry(1.571) ├─\n",
+ "├───────────┬╯├───┤╰───╯├───────────┤╭───╮├───────────┴╮\n",
+ "┤ ry(1.571) ├─┤ h ├──●──┤ rz(1.571) ├┤ x ├┤ ry(-1.571) ├\n",
+ "╰───────────╯ ╰───╯╭─┴─╮├───────────┤╰─┬─╯├───────────┬╯\n",
+ "───────────────────┤ x ├┤ rz(1.571) ├──●──┤ rx(1.571) ├─\n",
+ " ╰───╯╰───────────╯ ╰───────────╯ \n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "# the following code creates the circuit highlighted in excercise 3, which models the U(0,T) trotterizer function\n",
+ "# this function can be subdivided into two parts, into two qubit and four qubit rotations\n",
+ "# the two qubit block rotated around RYZ(theta) RZY(theta) \n",
+ "# this requires 2 entangling gates (4 gates each) and 4 single qubit gates\n",
+ "# the second block is the quad-qubit rotational block, rotated around RYZZZ(theta),RZYZZ(Theta),RZZYZ(theta\n",
+ "# and RZZZY(theta). \n",
+ "# the circuit breaks down the technique in quantum computing used to simulate the time evolution of a quantum system. \n",
+ "# It breaks down complex, large-scale Hamiltonian simulations into a series of smaller, manageable steps\n",
+ "\n",
+ "import cudaq\n",
+ "import numpy as np\n",
+ "import math as mat\n",
+ "\n",
+ "theta = 0 #in degrees, change this var after finding in part 4\n",
+ "pie_2 = mat.pi / 2\n",
+ "pie = mat.pi\n",
+ "\n",
+ "@cudaq.kernel\n",
+ "def dual_qubit_rotation():\n",
+ " q = cudaq.qvector(2)\n",
+ " rx(pie_2, q[1])\n",
+ " h(q[1])\n",
+ " h(q[0])\n",
+ " x.ctrl(q[0], q[1])\n",
+ " rz(theta, q[0])\n",
+ " rz(theta, q[1])\n",
+ " x.ctrl(q[1], q[0])\n",
+ " rx(pie_2, q[0])\n",
+ " rx(pie_2, q[1])\n",
+ " rx(pie_2, q[0])\n",
+ "\n",
+ "@cudaq.kernel\n",
+ "def quad_qubit_rotation():\n",
+ " q = cudaq.qvector(4)\n",
+ " ry(-pie_2, q[2])\n",
+ " ry(pie_2, q[1])\n",
+ " h(q[3])\n",
+ " rx(-pie_2, q[0])\n",
+ " h(q[2])\n",
+ " h(q[1])\n",
+ " h(q[0])\n",
+ " x.ctrl(q[0], q[1])\n",
+ " x.ctrl(q[2], q[3])\n",
+ " rz(-pie_2, q[0])\n",
+ " rz(-pie_2, q[1])\n",
+ " rz(-pie_2, q[2])\n",
+ " rz(-pie_2, q[3])\n",
+ " x.ctrl(q[1], q[0])\n",
+ " x.ctrl(q[2], q[3])\n",
+ " rx(pie_2, q[0])\n",
+ " rx(-pie_2, q[3])\n",
+ " ry(pie_2, q[2])\n",
+ " ry(-pie_2, q[1])\n",
+ " rx(-pie_2, q[1])\n",
+ " rx(pie_2, q[2])\n",
+ " h(q[1])\n",
+ " h(q[2])\n",
+ " x.ctrl(q[1], q[2])\n",
+ " rz(theta, q[1])\n",
+ " rz(theta, q[2])\n",
+ " x.ctrl(q[2], q[1])\n",
+ " rx(pie_2, q[1])\n",
+ " rx(pie, q[2])\n",
+ " ry(pie_2, q[1])\n",
+ " h(q[0])\n",
+ " h(q[1])\n",
+ " x.ctrl(q[0], q[1])\n",
+ " rz(pie_2, q[0])\n",
+ " rz(pie_2, q[1])\n",
+ " x.ctrl(q[1], q[0])\n",
+ " rx(pie_2, q[0])\n",
+ " ry(-pie_2, q[1])\n",
+ " h(q[1])\n",
+ " h(q[2])\n",
+ " x.ctrl(q[1], q[2])\n",
+ " rz(theta, q[1])\n",
+ " rz(theta, q[2])\n",
+ " x.ctrl(q[2], q[1])\n",
+ " rx(-pie, q[1])\n",
+ " ry(pie_2, q[2])\n",
+ " h(q[2])\n",
+ " h(q[3])\n",
+ " x.ctrl(q[2], q[3])\n",
+ " rz(-pie_2, q[2])\n",
+ " rz(-pie_2, q[3])\n",
+ " x.ctrl(q[3], q[2])\n",
+ " ry(-pie_2, q[2])\n",
+ " rx(-pie_2, q[3])\n",
+ " rx(-pie_2, q[2])\n",
+ " h(q[1])\n",
+ " h(q[2])\n",
+ " x.ctrl(q[1], q[2])\n",
+ " rz(theta, q[1])\n",
+ " rz(theta, q[2])\n",
+ " x.ctrl(q[2], q[1])\n",
+ " rx(pie_2, q[1])\n",
+ " rx(pie_2, q[2])\n",
+ " ry(-pie_2, q[1])\n",
+ " ry(pie_2, q[2])\n",
+ " h(q[0])\n",
+ " h(q[1])\n",
+ " h(q[2])\n",
+ " h(q[3])\n",
+ " x.ctrl(q[0], q[1])\n",
+ " x.ctrl(q[2], q[3])\n",
+ " rz(pie_2, q[0])\n",
+ " rz(pie_2, q[1])\n",
+ " rz(pie_2, q[2])\n",
+ " rz(pie_2, q[3])\n",
+ " x.ctrl(q[1], q[0])\n",
+ " x.ctrl(q[3], q[2])\n",
+ " ry(pie_2, q[1])\n",
+ " ry(-pie_2, q[2])\n",
+ " rx(pie_2, q[3])\n",
+ "\n",
+ "print(cudaq.draw(quad_qubit_rotation))\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f113f324-6f58-4b5e-ab93-339d70f88c46",
+ "metadata": {},
+ "source": [
+ "There are a few additional items we need to consider before completing the final implementation of the entire circuit. One simplification we can make is that for our problem the $h_i^x$ terms are all 1 and any $h_b^x$ terms are 0, and are only there for generalizations of this model. \n",
+ "\n",
+ "The remaining challenge is derivation of the angles that are used to apply each of the circuits you defined above. These are obtained from two terms $\\lambda(t)$ and $\\alpha(t)$. \n",
+ "\n",
+ "\n",
+ "The $\\lambda(t)$ defines an annealing schedule and is generally a Sin function which slowly \"turns on\" the problem Hamiltonian. For computing our angles, we need the derivative of $\\lambda(t)$.\n",
+ "\n",
+ "The $\\alpha$ term is a bit trickier and is the solution to a set of differential equations which minimize the distance between $H^1_{CD}(\\lambda)$ and $A(\\lambda)$. The result is \n",
+ "\n",
+ "$$\\alpha(t) = \\frac{-\\Gamma_1(t)}{\\Gamma_2(t)} $$\n",
+ "\n",
+ "Where $\\Gamma_1(t)$ and $\\Gamma_2(t)$ are defined in equations 16 and 17 of the paper and essentially depend on the structure of the optimization problem. Curious learners can look at the functions in `labs_utils.py` to see how these are computed, based on the problem size and specific time step in the Trotter process. \n",
+ "\n",
+ "\n",
+ "\n",
+ "
\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "id": "abf34168-42f9-4dbf-86e1-99976232ad7e",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "([(1, 2), (2, 3)], [(1, 3, 2, 4), (1, 4, 2, 5)])"
+ ]
+ },
+ "execution_count": 20,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "\n",
+ "\n",
+ "def two_body_pairs(N):\n",
+ " # two-body terms\n",
+ " \n",
+ " for i in range(1, N - 2):\n",
+ " max_k = (N - i)//2\n",
+ " \n",
+ " for k in range(1, max_k):\n",
+ " j = i + k\n",
+ " yield(i,j)\n",
+ "\n",
+ "\n",
+ "def four_body_pairs(N):\n",
+ " # four-body terms\n",
+ "\n",
+ " for i in range(1, N - 3):\n",
+ " max_t = (N - i - 1)//2\n",
+ "\n",
+ " for t in range(1, max_t):\n",
+ "\n",
+ " for k in range(t+1, N - i - t):\n",
+ " j = i + k\n",
+ " l = i + t\n",
+ " m = i + k + t\n",
+ " yield(i, j, l, m)\n",
+ "\n",
+ "\n",
+ "def get_interactions(N):\n",
+ " G2 = []\n",
+ " G4 = []\n",
+ " for i,j in two_body_pairs(N):\n",
+ " G2.append((i,j))\n",
+ "\n",
+ " for i,j,l,m in four_body_pairs(N):\n",
+ " G4.append((i,j,l,m))\n",
+ " \n",
+ " return G2, G4\n",
+ "\n",
+ "\n",
+ "\n",
+ "get_interactions(6)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "3450f200-b191-41f5-88d2-974052ee45ac",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "Exercise 4:
\n", + "\n", + "The `compute_theta` function, called in the following cells, requires all indices of the two and four body terms. These will be used again in our main kernel to apply the respective gates. Use the products in the formula below to finish the function in the cell below. Save them as `G2` and `G4` where each is a list of lists of indices defining the two and four term interactions. As you are translating an equation to a set of loops, this is a great opportunity to use an AI coding assistant.\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "\\begin{aligned}\n", + "U(0, T) = \\prod_{n=1}^{n_{\\text{trot}}} & \\left[ \\prod_{i=1}^{N-2} \\prod_{k=1}^{\\lfloor \\frac{N-i}{2} \\rfloor} R_{Y_i Z_{i+k}}\\big(4\\theta(n\\Delta t)h_i^x\\big) R_{Z_i Y_{i+k}}\\big(4\\theta(n\\Delta t)h_{i+k}^x\\big) \\right] \\\\\n", + "& \\times \\prod_{i=1}^{N-3} \\prod_{t=1}^{\\lfloor \\frac{N-i-1}{2} \\rfloor} \\prod_{k=t+1}^{N-i-t} \\bigg( R_{Y_i Z_{i+t} Z_{i+k} Z_{i+k+t}}\\big(8\\theta(n\\Delta t)h_i^x\\big) \\\\\n", + "& \\quad \\times R_{Z_i Y_{i+t} Z_{i+k} Z_{i+k+t}}\\big(8\\theta(n\\Delta t)h_{i+t}^x\\big) \\\\\n", + "& \\quad \\times R_{Z_i Z_{i+t} Y_{i+k} Z_{i+k+t}}\\big(8\\theta(n\\Delta t)h_{i+k}^x\\big) \\\\\n", + "& \\quad \\times R_{Z_i Z_{i+t} Z_{i+k} Y_{i+k+t}}\\big(8\\theta(n\\Delta t)h_{i+k+t}^x\\big) \\bigg)\n", + "\\end{aligned}\n", + "\\end{equation}\n", + "$$\n", + "\n", + "
\n",
+ "
\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 25,
+ "id": "8c2e6d7d-03b2-482f-bc36-d572e6d4855a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import cudaq\n",
+ "import numpy as np\n",
+ "from math import pi\n",
+ "\n",
+ "# RZZ gate implementation (needed for both blocks)\n",
+ "@cudaq.kernel\n",
+ "def rzz_gate(q0: cudaq.qubit, q1: cudaq.qubit, theta: float):\n",
+ " \"\"\"\n",
+ " Implements e^(-i*theta*Z⊗Z/2) using:\n",
+ " RZZ(θ) = CNOT(q0,q1) · RZ(θ,q1) · CNOT(q0,q1)\n",
+ " \"\"\"\n",
+ " x.ctrl(q0, q1)\n",
+ " rz(theta, q1)\n",
+ " x.ctrl(q0, q1)\n",
+ "\n",
+ "# Two-qubit rotation block R_YZ(θ) R_ZY(θ)\n",
+ "@cudaq.kernel\n",
+ "def apply_two_body_block(qubits: cudaq.qview, i: int, j: int, theta: float):\n",
+ " \"\"\"\n",
+ " Apply the 2-qubit block from Figure 3 in the paper.\n",
+ " This implements R_YiZj(θ) R_ZiYj(θ)\n",
+ " \n",
+ " Args:\n",
+ " qubits: quantum register\n",
+ " i, j: qubit indices (i < j)\n",
+ " theta: rotation angle (will be 4*θ(nΔt) in the main circuit)\n",
+ " \"\"\"\n",
+ " # First block: R_YZ(theta)\n",
+ " rx(pi/2, qubits[i]) # Basis change Y -> Z\n",
+ " h(qubits[j]) # Basis change Z -> X for CNOT\n",
+ " \n",
+ " rzz_gate(qubits[i], qubits[j], theta)\n",
+ " \n",
+ " rx(-pi/2, qubits[i]) # Restore basis\n",
+ " h(qubits[j])\n",
+ " \n",
+ " # Second block: R_ZY(theta)\n",
+ " h(qubits[i]) # Basis change Z -> X\n",
+ " rx(pi/2, qubits[j]) # Basis change Y -> Z\n",
+ " \n",
+ " rzz_gate(qubits[i], qubits[j], theta)\n",
+ " \n",
+ " h(qubits[i]) # Restore basis\n",
+ " rx(-pi/2, qubits[j])\n",
+ "\n",
+ "# Four-qubit rotation block\n",
+ "@cudaq.kernel\n",
+ "def apply_four_body_block(qubits: cudaq.qview, i: int, i_t: int, i_k: int, i_kt: int, theta: float):\n",
+ " \"\"\"\n",
+ " Apply the 4-qubit block from Figure 4 in the paper.\n",
+ " This implements the four rotations:\n",
+ " R_YZZZ, R_ZYZZ, R_ZZYZ, R_ZZZY\n",
+ " \n",
+ " Args:\n",
+ " qubits: quantum register\n",
+ " i, i_t, i_k, i_kt: four qubit indices for positions i, i+t, i+k, i+k+t\n",
+ " theta: rotation angle (will be 8*θ(nΔt) in the main circuit)\n",
+ " \"\"\"\n",
+ " indices = [i, i_t, i_k, i_kt]\n",
+ " \n",
+ " # The figure shows a complex decomposition with 10 RZZ gates and 28 single-qubit gates\n",
+ " # For each of the 4 rotations (Y on one qubit, Z on the others):\n",
+ " \n",
+ " # Rotation 1: Y on qubit i, Z on others\n",
+ " rx(pi/2, qubits[i]) # Y basis for qubit i\n",
+ " h(qubits[i_t])\n",
+ " h(qubits[i_k])\n",
+ " h(qubits[i_kt])\n",
+ " \n",
+ " # Multi-controlled Z rotation (simplified - actual circuit from paper is more complex)\n",
+ " # Apply pairwise RZZ gates\n",
+ " rzz_gate(qubits[i], qubits[i_t], theta/4)\n",
+ " rzz_gate(qubits[i], qubits[i_k], theta/4)\n",
+ " rzz_gate(qubits[i], qubits[i_kt], theta/4)\n",
+ " rzz_gate(qubits[i_t], qubits[i_k], theta/4)\n",
+ " rzz_gate(qubits[i_t], qubits[i_kt], theta/4)\n",
+ " rzz_gate(qubits[i_k], qubits[i_kt], theta/4)\n",
+ " \n",
+ " h(qubits[i_t])\n",
+ " h(qubits[i_k])\n",
+ " h(qubits[i_kt])\n",
+ " rx(-pi/2, qubits[i])\n",
+ " \n",
+ " # Rotation 2: Y on qubit i+t, Z on others\n",
+ " rx(pi/2, qubits[i_t])\n",
+ " h(qubits[i])\n",
+ " h(qubits[i_k])\n",
+ " h(qubits[i_kt])\n",
+ " \n",
+ " rzz_gate(qubits[i_t], qubits[i], theta/4)\n",
+ " rzz_gate(qubits[i_t], qubits[i_k], theta/4)\n",
+ " rzz_gate(qubits[i_t], qubits[i_kt], theta/4)\n",
+ " rzz_gate(qubits[i], qubits[i_k], theta/4)\n",
+ " rzz_gate(qubits[i], qubits[i_kt], theta/4)\n",
+ " rzz_gate(qubits[i_k], qubits[i_kt], theta/4)\n",
+ " \n",
+ " h(qubits[i])\n",
+ " h(qubits[i_k])\n",
+ " h(qubits[i_kt])\n",
+ " rx(-pi/2, qubits[i_t])\n",
+ " \n",
+ " # Rotation 3: Y on qubit i+k, Z on others\n",
+ " rx(pi/2, qubits[i_k])\n",
+ " h(qubits[i])\n",
+ " h(qubits[i_t])\n",
+ " h(qubits[i_kt])\n",
+ " \n",
+ " rzz_gate(qubits[i_k], qubits[i], theta/4)\n",
+ " rzz_gate(qubits[i_k], qubits[i_t], theta/4)\n",
+ " rzz_gate(qubits[i_k], qubits[i_kt], theta/4)\n",
+ " rzz_gate(qubits[i], qubits[i_t], theta/4)\n",
+ " rzz_gate(qubits[i], qubits[i_kt], theta/4)\n",
+ " rzz_gate(qubits[i_t], qubits[i_kt], theta/4)\n",
+ " \n",
+ " h(qubits[i])\n",
+ " h(qubits[i_t])\n",
+ " h(qubits[i_kt])\n",
+ " rx(-pi/2, qubits[i_k])\n",
+ " \n",
+ " # Rotation 4: Y on qubit i+k+t, Z on others\n",
+ " rx(pi/2, qubits[i_kt])\n",
+ " h(qubits[i])\n",
+ " h(qubits[i_t])\n",
+ " h(qubits[i_k])\n",
+ " \n",
+ " rzz_gate(qubits[i_kt], qubits[i], theta/4)\n",
+ " rzz_gate(qubits[i_kt], qubits[i_t], theta/4)\n",
+ " rzz_gate(qubits[i_kt], qubits[i_k], theta/4)\n",
+ " rzz_gate(qubits[i], qubits[i_t], theta/4)\n",
+ " rzz_gate(qubits[i], qubits[i_k], theta/4)\n",
+ " rzz_gate(qubits[i_t], qubits[i_k], theta/4)\n",
+ " \n",
+ " h(qubits[i])\n",
+ " h(qubits[i_t])\n",
+ " h(qubits[i_k])\n",
+ " rx(-pi/2, qubits[i_kt])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 27,
+ "id": "28648106-42d0-42d7-b52a-5383599e29d8",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "@cudaq.kernel\n",
+ "def dcqo_circuit(N: int, n_steps: int, G2: list[list[int]], G4: list[list[int]], \n",
+ " thetas_2body: list[float], thetas_4body: list[float]):\n",
+ " \"\"\"\n",
+ " Complete DCQO circuit for LABS problem.\n",
+ " \n",
+ " Args:\n",
+ " N: number of qubits (sequence length)\n",
+ " n_steps: number of Trotter steps\n",
+ " G2: list of 2-body interaction pairs [(i,j), ...]\n",
+ " G4: list of 4-body interaction quartets [(i,i+t,i+k,i+k+t), ...]\n",
+ " thetas_2body: rotation angles for 2-body terms [n_steps x len(G2)]\n",
+ " thetas_4body: rotation angles for 4-body terms [n_steps x len(G4)]\n",
+ " \"\"\"\n",
+ " # Allocate qubits\n",
+ " qubits = cudaq.qvector(N)\n",
+ " \n",
+ " # Step 1: Prepare initial state |+⟩^⊗N\n",
+ " # This is the ground state of H_i = Σ h_i^x σ_i^x\n",
+ " for i in range(N):\n",
+ " h(qubits[i]) # |0⟩ -> |+⟩\n",
+ " \n",
+ " # Step 2: Apply Trotter steps\n",
+ " for step in range(n_steps):\n",
+ " \n",
+ " # Apply all 2-body terms for this Trotter step\n",
+ " for idx in range(len(G2)):\n",
+ " i, j = G2[idx]\n",
+ " # Get theta for this step and interaction\n",
+ " theta_idx = step * len(G2) + idx\n",
+ " theta = thetas_2body[theta_idx]\n",
+ " \n",
+ " # Apply R_YZ and R_ZY rotations with angle 4*theta\n",
+ " apply_two_body_block(qubits, i, j, theta)\n",
+ " \n",
+ " # Apply all 4-body terms for this Trotter step\n",
+ " for idx in range(len(G4)):\n",
+ " i, i_t, i_k, i_kt = G4[idx]\n",
+ " # Get theta for this step and interaction\n",
+ " theta_idx = step * len(G4) + idx\n",
+ " theta = thetas_4body[theta_idx]\n",
+ " \n",
+ " # Apply the four rotations with angle 8*theta\n",
+ " apply_four_body_block(qubits, i, i_t, i_k, i_kt, theta)\n",
+ " \n",
+ " # Step 3: Measure all qubits\n",
+ " mz(qubits)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 28,
+ "id": "8b6673f7-832b-474e-b0e9-0182475b89d9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import auxiliary_files.labs_utils as utils\n",
+ "\n",
+ "def compute_theta(N, n_steps, T, G2, G4):\n",
+ " \"\"\"\n",
+ " Compute rotation angles for each Trotter step.\n",
+ " \n",
+ " Args:\n",
+ " N: number of qubits\n",
+ " n_steps: number of Trotter steps\n",
+ " T: total evolution time\n",
+ " G2: 2-body interaction indices\n",
+ " G4: 4-body interaction indices\n",
+ " \n",
+ " Returns:\n",
+ " thetas_2body: flat list of angles for 2-body terms\n",
+ " thetas_4body: flat list of angles for 4-body terms\n",
+ " \"\"\"\n",
+ " dt = T / n_steps # Time step size\n",
+ " \n",
+ " thetas_2body = []\n",
+ " thetas_4body = []\n",
+ " \n",
+ " for n in range(n_steps):\n",
+ " t = (n + 0.5) * dt # Mid-point of time step\n",
+ " \n",
+ " # Compute λ(t) and its derivative\n",
+ " lambda_t = utils.lambda_schedule(t, T)\n",
+ " lambda_dot = utils.lambda_derivative(t, T)\n",
+ " \n",
+ " # Compute α(t) from Γ₁ and Γ₂\n",
+ " gamma1 = utils.compute_gamma1(N, G2, G4, lambda_t)\n",
+ " gamma2 = utils.compute_gamma2(N, G2, G4, lambda_t)\n",
+ " alpha_t = -gamma1 / gamma2\n",
+ " \n",
+ " # θ(t) = Δt × α(t) × λ̇(t)\n",
+ " theta_t = dt * alpha_t * lambda_dot\n",
+ " \n",
+ " # For 2-body terms: angle is 4*θ(t) (but we pass θ and multiply in kernel)\n",
+ " for _ in G2:\n",
+ " thetas_2body.append(4 * theta_t)\n",
+ " \n",
+ " # For 4-body terms: angle is 8*θ(t)\n",
+ " for _ in G4:\n",
+ " thetas_4body.append(8 * theta_t)\n",
+ " \n",
+ " return thetas_2body, thetas_4body"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 29,
+ "id": "306eac9b-5965-4bbe-96c7-4671943c9617",
+ "metadata": {},
+ "outputs": [
+ {
+ "ename": "AttributeError",
+ "evalue": "module 'auxiliary_files.labs_utils' has no attribute 'lambda_schedule'",
+ "output_type": "error",
+ "traceback": [
+ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+ "\u001b[0;31mAttributeError\u001b[0m Traceback (most recent call last)",
+ "Cell \u001b[0;32mIn[29], line 11\u001b[0m\n\u001b[1;32m 8\u001b[0m G2, G4 \u001b[38;5;241m=\u001b[39m get_interactions(N)\n\u001b[1;32m 10\u001b[0m \u001b[38;5;66;03m# Compute rotation angles\u001b[39;00m\n\u001b[0;32m---> 11\u001b[0m thetas_2body, thetas_4body \u001b[38;5;241m=\u001b[39m \u001b[43mcompute_theta\u001b[49m\u001b[43m(\u001b[49m\u001b[43mN\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mn_steps\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mT\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mG2\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mG4\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 13\u001b[0m \u001b[38;5;66;03m# Run the quantum circuit\u001b[39;00m\n\u001b[1;32m 14\u001b[0m result \u001b[38;5;241m=\u001b[39m cudaq\u001b[38;5;241m.\u001b[39msample(dcqo_circuit, N, n_steps, G2, G4, thetas_2body, thetas_4body, \n\u001b[1;32m 15\u001b[0m shots_count\u001b[38;5;241m=\u001b[39mn_shots)\n",
+ "Cell \u001b[0;32mIn[28], line 27\u001b[0m, in \u001b[0;36mcompute_theta\u001b[0;34m(N, n_steps, T, G2, G4)\u001b[0m\n\u001b[1;32m 24\u001b[0m t \u001b[38;5;241m=\u001b[39m (n \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m0.5\u001b[39m) \u001b[38;5;241m*\u001b[39m dt \u001b[38;5;66;03m# Mid-point of time step\u001b[39;00m\n\u001b[1;32m 26\u001b[0m \u001b[38;5;66;03m# Compute λ(t) and its derivative\u001b[39;00m\n\u001b[0;32m---> 27\u001b[0m lambda_t \u001b[38;5;241m=\u001b[39m \u001b[43mutils\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mlambda_schedule\u001b[49m(t, T)\n\u001b[1;32m 28\u001b[0m lambda_dot \u001b[38;5;241m=\u001b[39m utils\u001b[38;5;241m.\u001b[39mlambda_derivative(t, T)\n\u001b[1;32m 30\u001b[0m \u001b[38;5;66;03m# Compute α(t) from Γ₁ and Γ₂\u001b[39;00m\n",
+ "\u001b[0;31mAttributeError\u001b[0m: module 'auxiliary_files.labs_utils' has no attribute 'lambda_schedule'"
+ ]
+ }
+ ],
+ "source": [
+ "# Problem parameters\n",
+ "N = 8 # Sequence length (start small for testing)\n",
+ "n_steps = 10 # Number of Trotter steps\n",
+ "T = 1.0 # Total evolution time\n",
+ "n_shots = 1000 # Number of measurements\n",
+ "\n",
+ "# Get interaction indices from Exercise 4\n",
+ "G2, G4 = get_interactions(N)\n",
+ "\n",
+ "# Compute rotation angles\n",
+ "thetas_2body, thetas_4body = compute_theta(N, n_steps, T, G2, G4)\n",
+ "\n",
+ "# Run the quantum circuit\n",
+ "result = cudaq.sample(dcqo_circuit, N, n_steps, G2, G4, thetas_2body, thetas_4body, \n",
+ " shots_count=n_shots)\n",
+ "\n",
+ "# Extract bitstrings and energies\n",
+ "bitstrings = []\n",
+ "energies = []\n",
+ "\n",
+ "for bits, count in result.items():\n",
+ " # Convert bitstring to list\n",
+ " bit_list = [int(b) for b in bits]\n",
+ " \n",
+ " # Compute LABS energy\n",
+ " energy = labs_energy_list(bit_list)\n",
+ " \n",
+ " # Store\n",
+ " for _ in range(count):\n",
+ " bitstrings.append(bit_list)\n",
+ " energies.append(energy)\n",
+ "\n",
+ "# Find best solution from quantum sampling\n",
+ "best_idx = np.argmin(energies)\n",
+ "best_quantum = bitstrings[best_idx]\n",
+ "best_energy = energies[best_idx]\n",
+ "\n",
+ "print(f\"Best quantum solution: {best_quantum}\")\n",
+ "print(f\"Energy: {best_energy}\")\n",
+ "\n",
+ "# Use this to seed MTS (Exercise 6)\n",
+ "quantum_population = [best_quantum] * 100 # Replicate best solution\n",
+ "# Or use diverse solutions:\n",
+ "# quantum_population = bitstrings[:100]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "fb89d90e-66e2-4700-85b9-40df9fca22c1",
+ "metadata": {},
+ "source": [
+ "## Generating Quantum Enhanced Results\n",
+ "\n",
+ "Recall that the point of this lab is to demonstrate the potential benefits of running a quantum subroutine as a preprocessing step for classical optimization of a challenging problem like LABS. you now have all of the tools you need to try this for yourself.\n",
+ "\n",
+ "Exercise 5:
\n", + "\n", + "You are now ready to construct the entire circuit and run the counteradiabatic optimization procedure. The final kernel needs to apply Equation 15 for a specified total evolution time $T$ and the `n_steps` number of Trotter steps. It must also take as input, the indices for the two and four body terms and the thetas to be applied each step, as these cannot be computed within a CUDA-Q kernel.\n", + "\n", + "The helper function `compute_theta` computes the theta parameters for you, using a few additional functions in accordance with the equations defined in the paper.\n", + "
\n",
+ "
\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "8e5f02e6-41bf-4634-9cb1-f0f543ef2e3f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# TODO - write code here to sample from your CUDA-Q kernel and used the results to seed your MTS population"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "0a5756e2-e4cb-42f4-a4b3-7f627095f4f4",
+ "metadata": {},
+ "source": [
+ "The results clearly show that a population sampled from CUDA-Q results in an improved distribution and a lower energy final result. This is exactly the goal of quantum enhanced optimization. To not necessarily solve the problem, but improve the effectiveness of state-of-the-art classical approaches. \n",
+ "\n",
+ "A few major caveats need to be mentioned here. First, We are comparing a quantum generated population to a random sample. It is quite likely that other classical or quantum heuristics could be used to produce an initial population that might even beat the counteradiabatic method you used, so we cannot make any claims that this is the best. \n",
+ "\n",
+ "Recall that the point of the counteradiabatic approach derived in the paper is that it is more efficient in terms of two-qubit gates relative to QAOA. The benefits of this regime would only truly come into play in a setting (e.g. larger problem instance) where it is too difficult to produce a good initial population with any know classical heuristic, and the counteradiabatic approach is more efficiently run on a QPU compared to alternatives.\n",
+ "\n",
+ "We should also note that we are comparing a single sample of each approach. Maybe the quantum sample got lucky or the randomly generated population was unlucky and a more rigorous comparison would need to repeat the analysis many times to draw any confidently conclusions. \n",
+ "\n",
+ "The authors of the paper discuss all of these considerations, but propose an analysis that is quite interesting related to the scaling of the technique. Rather than run large simulations ourselves, examine their results below. \n",
+ "\n",
+ "\n",
+ "Exercise 6:
\n", + "\n", + "Use your CUDA-Q code to prepare an initial population for your memetic search algorithm and see if you can improve the results relative to a random initial population. If you are running on a CPU, you will need to run smaller problem instances. The code below sets up the problem\n", + "\n", + "
![](\"images/quantum_enhanced_optimization_LABS/tabu_search_results.png\")
\n",
+ "\n",
+ "The authors computed replicate median (median of solving the problem repeated with same setup) time to solutions (excluding time to sample from QPU) for problem sizes $N=27$ to $N=37$. Two interesting conclusions can be drawn from this. First, standard memetic tabu search (MTS) is generally faster than quantum enhanced (QE) MTS. But there are two promising trends. For larger problems, the QE-MTS experiments occasionally have excellent performance with times to solution much smaller than all of the MTS data points. These outliers indicate there are certain instances where QE-MTS could provide much faster time-to-solution. \n",
+ "\n",
+ "More importantly, if a line of best fit is calculated using the median of each set of medians, the slope of the QE-MTS line is smaller than the MTS! This seems to indicate that QE solution of this problem scales $O(1.24^N)$ which is better than the best known classical heuristic ($O(1.34^N)$) and the best known quantum approach (QAOA - $O(1.46^N)$).\n",
+ "\n",
+ "For problems of size of $N=47$ or greater, the authors anticipate that QE-MTS could be a promising technique and produce good initial populations that are difficult to obtain classically. \n",
+ "\n",
+ "The study reinforces the potential of hybrid workflows enhanced by quantum data such that a classical routine is still the primary solver, but quantum computers make it much more effective. Future work can explore improvements to both the quantum and classical sides, such as including GPU accelerated memetic search on the classical side."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "7aab7af9",
+ "metadata": {
+ "jp-MarkdownHeadingCollapsed": true
+ },
+ "source": [
+ "## Self-validation To Be Completed for Phase 1\n",
+ "\n",
+ "In this section, explain how you verified your results. Did you calculate solutions by hand for small N? Did you create unit tests? Did you cross-reference your Quantum energy values against your Classical MTS results? Did you check known symmetries?"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 [cuda-q-v0.13.0]",
+ "language": "python",
+ "name": "python3_zogs"
+ },
+ "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.11.9"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/tutorial_notebook/01_quantum_enhanced_optimization_LABS.ipynb b/tutorial_notebook/01_quantum_enhanced_optimization_LABS.ipynb
index f96cb45c..fc76d7d4 100644
--- a/tutorial_notebook/01_quantum_enhanced_optimization_LABS.ipynb
+++ b/tutorial_notebook/01_quantum_enhanced_optimization_LABS.ipynb
@@ -145,7 +145,176 @@
"metadata": {},
"outputs": [],
"source": [
- "#TODO - Write code to perform MTS"
+ "#TODO - Write code to perform MTS\n",
+ "\n",
+ "import cudaq\n",
+ "import numpy as np\n",
+ "from math import floor\n",
+ "import auxiliary_files.labs_utils as utils\n",
+ "\n",
+ "\n",
+ "import random \n",
+ "\n",
+ "#creates bitstrings of length N\n",
+ "def bitstring(N):\n",
+ " string = []\n",
+ " for bit in range(N):\n",
+ " string.append(random.choice([0,1])) \n",
+ " return string\n",
+ "\n",
+ "\n",
+ "\n",
+ "#creates s_max bitstrings of lenth N\n",
+ "def populate(s_max,N):\n",
+ " pop = []\n",
+ " for i in range(s_max):\n",
+ " pop.append(bitstring(N))\n",
+ " return pop\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "#makes a array from 0 to N\n",
+ "def cut(N):\n",
+ " k_list = []\n",
+ " for bit in range(N):\n",
+ " k_list.append(bit) \n",
+ " return k_list\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "#cuts two parent bitstrings p1 and p2 of length N at index k and combines the halves into a new bitstring child of length N\n",
+ "def combine(p1,p2):\n",
+ " if len(p1)==len(p2):\n",
+ " k = random.choice(cut(len(p1)))\n",
+ " child = p1[0:k]\n",
+ " for element in p2[k:len(p2)]:\n",
+ " child.append(element)\n",
+ " return child\n",
+ "\n",
+ "\n",
+ "\n",
+ "#pm probability of flippinng each bit in a bitstring s \n",
+ "def mutate(s,pm):\n",
+ " for i in range(len(s)):\n",
+ " if random.random()